T-REX: the tile-based representation of lateral exchange processes in ICON-Land

The vast majority of land-surface models uses a tiling-approach to capture the effects of subgrid-scale spatial heterogeneity in the land-surface properties. In most cases, however, the tiles, which represent patches with homogeneous characteristics at and below the surface, are treated independently of each other and the lateral exchange between them is not being taken into consideration. The present manuscript describes an approach for a tile-based representation of lateral exchange processes in heterogeneous landscapes that was recently implemented into ICON-Land, the land surface component of the ICON framework. The scheme captures the horizontal fluxes on a broad range of spatial scales and represents 5 lateral exchange processes, namely gravity-driven moisture fluxes and the corresponding advective heat transport, diffusive- and conductive fluxes of water and heat as well as the redistribution of snow between small-scale topographic features. In the approach, the relationships between any two tiles are determined by a set of characteristic connectivities which are treated as inherent properties of the pair of tiles – invariant in time and independent of the location – and derive from the internal logic underlying the definition of the tiles. The characteristic connectivities are used to calculate the spatio-geometric relationships between two tiles, such as the (geometrical) contact length or the characteristic center-to-center distance between two adjacent surface clusters. These, in turn, define gradients as well as the time-lag factors that govern the lateral transport in the model. In addition to a description of the model development, we present two example applications of the new scheme, which address the effect of sub-grid scale fluxes on the model's ability to capture the spatial variability in the state of the surface and sub-surface and the overall terrestrial water storage. Here, our results suggest that lateral exchange processes, especially on small horizontal scales, are highly relevant for the spatial variability in the soil temperatures and for the simulated extent of surface water bodies, while the effects on the grid-cell mean state and the turbulent exchange with the atmosphere appear to be largely negligible.

1 Introduction

Most land surface models (LSM), especially those included in Earth system models, have been designed to run with horizontal grid spacings of tens to hundreds of kilometers. While such resolutions might cover an important part of the spatial variability in the atmospheric state variables, this is not necessarily the case for soil- and land-surface properties which can vary at the (sub-) meter scale. Here, the coarse resolution may entail aggregation errors and hinders the representation of distinct ecosystem states which shape the physical and biochemical processes at and below the surface (Rastetter et al., 1992; Beer, 2016). Even with the grid spacing of some applications approaching the kilometer-scale (Stevens et al., 2019; Bogenschutz et al., 2023; Lee and Hohenegger, 2024), many of the relevant processes can still not be resolved. To deal with the issue, most models employ a tiling-approach (Avissar and Pielke, 1989; Koster and Suarez, 1992). This approach is based on the assumption that the land surface within one grid cell can be split into sub-divisions, which each comprise patches with homogeneous characteristics, and that model processes can be determined for each of these sub-units, the tiles, separately. Most often, the individual tiles do not interact directly, but only through the exchange with the overlying atmosphere (Molod et al., 2003; Best et al., 2004; de Vrese et al., 2016b; Huang et al., 2022), allowing for a computationally inexpensive treatment of the subgrid-scale variability since the landscape heterogeneity does not need to be be resolved in a spatially explicit way.

Land cover constitutes a key determinant of the local energy and water cycles, with vegetation affecting virtually all terms of the surface energy balance and the moisture exchange with the atmosphere (Seneviratne et al., 2010; Duveiller et al., 2018). Topography is another important factor of heterogeneity and induces variability on all scales. The microtopography often determines the position of the local water table and soil temperatures, which, in turn, control soil respiration rates and their sensitive to fluctuations in the controlling variables (Sommerkorn, 2008). At the other end of the scale-range, topographic ridge-to-valley gradients are the main cause for a moisture convergence at the hillslope- and catchment-scale (Fan et al., 2019). And while topography and land cover arguably constitute the more prominent factors, soil properties may also have a strong impact on the local state of the surface and subsurface. Soil structural properties are closely linked to infiltration rates, the water holding capacity (Bordoloi et al., 2018; Basset et al., 2023) and the soil thermophysical properties (Zhang and Wang, 2017). Here, horizontal gradients in the soil organic matter content, in particular, have a large potential for inducing sub-grid-scale variability in soil moisture as well as in surface and below-ground temperatures (Rawls et al., 2004; Koven et al., 2009; Fekete et al., 2012; Chadburn et al., 2015; Langer et al., 2016; Zhu et al., 2019; Siewert et al., 2021; Zhang et al., 2021).

With land cover being one of the most important sources of heterogeneity, the tiles have traditionally been used to represent different vegetation classes or plant functional types. These are key determinants of the properties shaping the physical land-atmosphere interactions – albedo, roughness and surface resistance to evapotranspiration – as well as for the processes governing the terrestrial carbon cycle. However, studies have shown that, depending on the region and processes considered, other factors may be as important as the vegetation cover, with one of the most prominent examples being agricultural irrigation (de Rosnay et al., 2003; Boucher et al., 2004; Sacks et al., 2009; Guimberteau et al., 2011; Puma and Cook, 2010; de Vrese et al., 2016a; de Vrese and Hagemann, 2017; Chou et al., 2018; Singh et al., 2018; Hauser et al., 2019; Cook et al., 2020; Al‐Yaari et al., 2022; McDermid et al., 2023). Here, the soil moisture in well-confined areas is maintained at high levels often in (semi-) arid regions and several models employ tiles to distinguish between managed and unmanaged agricultural areas (Yao et al., 2025).

Since the tiles possess distinct properties, they also develop distinct hydrological and thermophysical states, e.g. with average surface temperature differences between irrigated and non-irrigated grid-cell fractions and between grass- and tree tiles being in the order of degrees (Schultz et al., 2016; Thiery et al., 2017). In reality, numerous processes act to both reduce and increase horizontal gradients in the state of the surface and the subsurface, with most models (implicitly) making extreme assumptions on the effectiveness of these processes. As stated above, the most prominent assumption is that there is no direct interaction between the tiles, other than through the exchange with the atmosphere. The other extreme assumption, which is used e.g. in JSBACH3 the land surface component of the Max Planck Institute for Meteorology's MPI-ESM and the Alfred Wegener Institute's AWI-ESM, is that the lateral exchange is highly effective so that the physical state variables, most importantly soil moisture as well as soil and surface temperatures, are continuously averaged between tiles (Reick et al., 2021). It is obvious that these assumptions are only valid for either extremely large or extremely small horizontal scales but misrepresent a large fraction of the real-world landscape heterogeneity. A broadly applicable representation of subgrid-scale heterogeneity requires an explicit representation of the lateral exchange between tiles and many advances have been made to incorporate such interactions, especially with respect to the hydrological fluxes (Aas et al., 2019; Fan et al., 2019; Swenson et al., 2019; Blyth et al., 2021; Chaney et al., 2021; Smith et al., 2022; de Vrese et al., 2024; Li et al., 2024).

In the following, we describe a recently developed scheme for including horizontal subgrid-scale fluxes of water and heat in ICON-Land, the land surface component used within the ICON modelling framework (Jungclaus et al., 2022); namely a tile-based representation of lateral exchange processes in heterogeneous landscapes (T-REX). First, we give a concise overview of the assumptions and the governing equations of the scheme (Sect. 2), before investigating the effects that the lateral exchange processes have on the simulated state of the surface and subsurface, using two example applications (Sect. 3). Here, the first application targets the effect of lateral subgrid-scale fluxes on the terrestrial water storage as well as the resulting impact on the turbulent land-atmosphere exchange and the state of the surface (Sect. 3.1). With the second example, we investigate the effect of the lateral heat transport on the spatial subgrid-scale variability in the soil temperatures, focusing on typical patterned-ground structures commonly found in periglacial regions (Set. 3.2).

2 T-REX

T-REX is designed to enable simulations in which the prevalent spatial variability in the state variables and the exchange fluxes with the atmosphere can be captured by the tiling scheme. Thus, the model needs to be capable of representing the interactions between all those tiles that are required to resolve the subgrid-scale heterogeneity in the determining factors of the local hydrological conditions and the thermophysical state at and below the land surface. In reality, the spatial variability comprises a wide range of length-scales (Blöschl and Sivapalan, 1995; Seyfried and Wilcox, 1995), requiring the model to determine the lateral transport processes for a similarly broad range of scales (Fig. 1): On the micro scale, moisture conditions can vary dramatically across distances of a few meters, which can also entail pronounced temperature gradients, especially if there is a strong effect on the soil thermophysical properties – that is heat conductivity and capacity – and on the energy released by or required for phase changes in the soil (Langer et al., 2011). Macro-scale moisture convergence – that is on length scales of hundreds of meters to a few kilometers – may result in pronounced vegetation gradients and landscape types such as gallery forests and raised bogs (Li et al., 2024), with a large potential of macro-scale soil moisture differences to notably affect temperatures via their effect on the evaporative cooling (Thiery et al., 2017).

Figure 1 Spatial subgrid-scale heterogeneity. Shown is an idealized representation of the spatial heterogeneity that is not resolved by the horizontal grid and how this is treated by the tiling scheme of the model. Here, the land surface may exhibit heterogeneity with respect to a number of characteristics, most importantly differences in topography, land cover or soil properties. Small scale features (here, represented by the model by the tiles 2 and 4), that is patches or clusters with horizontal length scale of a few meters to tenths of meters, interact with the larger features they are encompassed in (tiles 1 and 3). Depending on the kind of heterogeneity, non-negligible interactions can also persist across distances of hundreds of meters to a few kilometers, that is between the surface clusters represented by tile 1 and tile 3. The model does not treat the clusters or patches individually but aggregates them into tiles, that are represented by a cover fraction and by a set of characteristic properties. The state of the clusters and the lateral exchange between them is then determined based on the characteristic properties that describe the tiles themselves and also the connections between each pair of tiles.

T-REX includes 5 lateral exchange processes that can be categorized with respect to the transported quantity (water or heat), the scales the processes are relevant for – that is the micro-scale (horizontal distances of roughly 1 to 100 m) and the macro-scale (roughly 100 to 10 000 m) – and the mechanism driving the transport (gravity-driven-, wind-driven- and conductive/diffusive fluxes) (Fig. 2). On the micro scale the scheme represents the diffusive- and conductive fluxes of water and heat (Sects. 2.3.1 and 2.4.1) as well as a wind-driven redistribution of snow (Sect. 2.3.3). Furthermore, the model accounts for gravity-driven fluxes of water, which we simply refer to as surface- and subsurface runoff (Sect. 2.3.2), and the corresponding advective heat transport (Sect. 2.4.2). A diffusive and conductive macro-scale exchange is not included in the model, with gradients being negligible due to the large distances considered. On this scale, the lateral transport is limited to the gravity-driven fluxes of water and the resulting advective heat fluxes. In general, the representation of the transport processes is based on a set of factors that describe the degree to which two tiles are connected (Sect. 2.1) in combination with a number of explicit or implicit assumptions on the spatial relationship between the tiles (Sect. 2.2). In our technical implementation, the interaction between two tiles is performed sequentially. This involves computing fluxes and gradients between an individual tile, t, and an equivalent connected tile, which we refer to as a sibling tile, s.

Figure 2 Lateral transport processes. Overview of the five lateral transport processes that are included in T-REX, grouped by the transported quantity, the underlying mechanisms and the scales at which the exchange is assumed to have a notable effect on the physical state of the land surface. Concerning the redistribution of snow, it should be noted that, in reality, wind is the dominant force responsible for the movement across the land surface. However, due to structural limitations of our land-surface model, wind speed and direction cannot be taken into consideration, and the scheme merely parameterizes the bulk effect of primarily wind-driven snow redistribution between small-scale topographic features (see Sect. 2.3.3).

Finally, we would like to highlight that numerous approaches exist to describe any of the processes that the model represents. Some are more simplistic than others and they may also vary substantially with respect to their computational costs and data requirements. For the present implementation, one of our main goals was to design the parametrizations as consistent as possible with the assumptions that ICON-Land's soil-hydrology- and thermophyiscs routines are based on, since the lateral exchange processes are tightly coupled to the vertical movement of water and energy. Here, however, T-REX should be thought of less as a final suit of parametrizations of lateral exchange processes and more of a framework that allows for a consistent treatment of the inter-tile exchange, with enough flexibility to facilitate the implementation of improved parametrizations and new processes.

2.1 Connectivities

Representing subgrid-scale heterogeneity in land-surface models requires balancing realism with computational complexity and -demand as well as with the availability of data needed to describe interactions between surface elements. In this first T-REX implementation, this trade-off is addressed by adopting a tiling approach that focuses on a small number of generalized landscape elements whose interactions follow a broadly applicable physical logic. Consequently, it is most naturally suited to simplified and idealized tile configurations and not necessarily designed to capture every aspect of spatial heterogeneity that may be important for a specific region. Instead, it targets landscape elements that are expected to play a key role across a wide range of environments, and whose interactions follow a consistent structural logic that can be applied throughout the model domain. This approach allows T-REX to start from simple configurations that are straightforward to set up for idealized cases (see Sect. 3.1 and 3.2) while providing a framework that can be extended to increase the degree of realism in future simulations (see Sect. 4).

The main requirement of T-REX is that the relationships between any two tiles must be describable by a set of characteristic connectivities that reflect their interactions reasonably well. For example, surface runoff generated on elevated terrain will generally flow toward surrounding lowland areas via slopes, suggesting that any grid-cell can be subdivided into elevations, slopes and lowlands, with elevations being hydrologically connected to the slopes which in turn are connected to the lowlands. These characteristic connectivities are used to calculate the spatial relationships between two tiles (see Sect. 2.2) as well as the time-lag factors (see Sect. 2.3.2) that govern the lateral transport (note that all key variables are summarized in Appendix A). In the present implementation, they are treated as inherent properties of the pair of tiles, invariant in time and independent of the location and composition of a grid cell. Here, the characteristic connectivities merely describe the general connection between two types of surface clusters, based on the internal logic underlying the tile definition of a given setup. However, the characteristic connectivities do not necessarily constitute the degree to which two tiles are connected within a specific grid cell. Since in most cases the inter-tile exchange ultimately depends on the grid-cell composition, the actual connectivities also account for the fractional cover of the tiles within a given grid cell (see below). The characteristic connectivities ${\mathrm{\Lambda }}_{\mathrm{p}}^{\mathrm{t},\mathrm{s}}$ between a tile t and its siblings s need to be provided to the model in the form of a n×n connection matrix, where n is the number of tiles considered in the setup. T-REX uses several matrices to describe the lateral interactions between tiles, since the connectives depend on the process p that is being considered (Table 1).

Table 1 Types of connectivities considered by T-REX (note that superscripts referring to the connected tile pair are dropped).

The most basic of these matrices provides the characteristic connectivities due to direct contact between two tiles (Λ_ctc). Here, the matrix elements (${\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}$) constitute the relative (geometrical) contact lengths, corresponding to the fraction of the boundary of a surface patch of tile t that interfaces with a patch of given sibling tile s. These connectivities are used to determine the spatial relationships between a tile and its siblings and all those lateral fluxes that are directly proportional to the inter-tile gradient in the corresponding state variables – that is the conductive and diffusive fluxes of heat and water – as well as the redistribution of snow. The matrix elements need to be provided in such a way that (i) if a connectivity ${\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}>\mathrm{0.0}$ is assumed between t and s, a connectivity ${\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{s},\mathrm{t}}>\mathrm{0.0}$ may not be given for the relation between s and t, since this would result in certain fluxes being computed twice per time-step for the same connection. Here, the order should be in such a way that (ii) the connectivity is provided for t, with t being the tile that is encompassed by sibling s or in a way that t is the tile with the smaller characteristic area (see Sect. 2.2.1). (iii) The sum of all connectivities that are defined for tile t ($\sum _{s=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}$, with nct being the number of connected tiles) may not exceed 1., while the sum of all its connectivities, including those that have been defined for the tile's siblings ($\sum _{s=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}+{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{s},\mathrm{t}}$) may be smaller than 1, implying that a given fraction of the patch's boundary interfaces with non-interacting surface patches.

As stated above, the characteristic connectivities ${\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}$ alone may not be sufficient to accurately describe the effective connectivity between a tile and a sibling within a specific grid cell, particulalrly for tiles with multiple connections. Here, an important factor that is being neglected is the relative abundance of the connected siblings within the cell which, for many factors, can be regarded as a reasonable proxy for the likelihood that a larger fraction the tile is connected to a specific sibling rather than to the other siblings. Thus, neglecting the grid-cell composition can produce highly implausible fluxes – for example, a sibling occupying 70 % of the area may receive only half of the runoff generated on the tile, while a second sibling also receives half of the runoff but occupies only 1 % of the grid cell. To avoid such situations, the effective connectivity is assumed to be proportional to the relative abundance of the connected sibling within the cell. In the standard setup, this is implemented by weighting the characteristic connectivities by the cover fraction of each connected sibling (c^s), while preserving the overall geometric assumptions of the original connectivity matrix:

$$\begin{array}{}\text{(1)}& {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{t},\mathrm{s}}={\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}\cdot c^{\mathrm{s}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},k}\cdot c^k}}\cdot \sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},k}.\end{array}$$

This weighting is optional and only relevant when a tile interacts with multiple siblings and, if the connectivity between tiles is inherently independent of the grid cell composition, ${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{t},\mathrm{s}}$ can simply be set equal to ${\mathrm{\Lambda }}_{\mathrm{ctc}}^{\mathrm{t},\mathrm{s}}$.

With respect to surface- and subsurface runoff and the corresponding advective heat fluxes, it is assumed, that the movement of water does not happen uniformly across the surface or through the soil, but through dominant hydrological flow paths. Consequently, the runoff fluxes from a tile to its connected siblings may not be directly proportional to the contact lengths. A good example for this are tiles that are connected to two siblings with different surface elevations. Here, the relative contact length may be the same with respect to the two siblings but water at the surface would only run off to the sibling that has a lower surface elevation than the tile. The connectivities are distinct for the micro- and macro-scale fluxes as well as for the surface and subsurface flow. Consequently, the present T-REX setup considers four additional matrices, representing micro-scale connectivities at the surface (${\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{micro},\mathrm{srf}}^{\mathrm{t},\mathrm{s}}$) and below ground (${\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{micro},\mathrm{blg}}^{\mathrm{t},\mathrm{s}}$) as well as macro-scale connectivities at the surface (${\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{macro},\mathrm{srf}}^{\mathrm{t},\mathrm{s}}$) and below ground (${\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{macro},\mathrm{blg}}^{\mathrm{t},\mathrm{s}}$). On the macro-scale, connections may even exist for two tiles that do not interface directly, but where the hydrological flow paths are assumed to pass the area of another tile without interacting with the latter. This simplification was required to be able to represent runoff across hillslopes - where water may initially pool in micro-scale depressions – with only a small number of tiles. Hillslope runoff may best be described by a depressional fill–spill cascade, in which a fraction of the runoff initially pools in small scale depressions until it overflows and moves down-slope across the even surface of the slope to again pool in a lower lying depression (McDonnell et al., 2021). For many applications, it may, however, not be feasible to resolve the lateral fluxes representing individual height bands, as proposed by e.g. Chaney et al. (2021), and the entire hillslope is subdivided into two tiles, representing the relatively even slope surface – that is the areas where water may not pool – and small scale depressions. In such a setup, the model can not represent hill-slope flow cascades and the inclusion of a flux from the depression tile back to the tile representing the even surface fraction results in a circular connection that has the potential to retain water on the slope indefinitely. To avoid such closed loops, the runoff from the depression tile can be routed to a connected lowland tile directly, bypassing the tile representing the relatively even surface areas of the slope.

Table 2 Key variables used to describe the spatial relationship between connected tiles.

2.2 Spatial relationships

On the macro scale, all spatial relationships between connected tiles are merely implicitly included in the parametrizations of the model, via the slopes and the time-lag factors used to determine the macro-scale runoff (see Sect. 2.3.2). On the micro scale, however, the lateral exchange between a tile and its connected sibling tiles is in large parts determined by spatial relationships that the model determines explicitly, more specifically by differences in surface elevation, the distance between the centers of the patches that the tiles represent and the characteristic interface length of the connection. These are determined based on the characteristic surface elevation, the characteristic area covered by a patch of a given tile, the fraction of the area of the patch that is typically connected to a patch of a specific sibling-tile and by the cover fractions of the tiles, all of which need to be provided as input parameters. A concise overview of the key variables used to describe the spatial relationship between a given pair of tiles is provided in Table 2.

2.2.1 Horizontal relationships

The scheme assumes a constellation in which one of two (micro-scale-) connected tiles, the outer tile (t^o), encompasses the second, the inner tile (tⁱ), at least partially (Fig. 3). Furthermore, it is assumed that tⁱ represents well distributed subgrid patches, the shape of which can be approximated by a circle (oⁱ) with the characteristic area Aⁱ. This allows the calculation of the characteristic interface length between tiles tⁱ and t^o, l^i,o, as the circular arc whose fraction of oⁱ's circumference is equal to the relative contact length (${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{i},\mathrm{o}}$).

$$\begin{array}{}\text{(2)}& l^{\mathrm{i},\mathrm{o}}={\left(A^{\mathrm{i}}\cdot \mathit{\pi }\cdot \mathrm{4}\right)}^{\mathrm{0.5}}\cdot {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{i},\mathrm{o}}.\end{array}$$

oⁱ forms the inner region of a second circle (o^a) which contains the characteristic areas of tⁱ and of t^o, and whose radius (r^a) is used to approximate the center-to-center distance between the two tiles (d^i,o). The latter is assumed to be equal to half the distance between the circumferences of oⁱ and o^a added to rⁱ, the radius of oⁱ. The area of o^a is determined so that a vector with a central angle α (in radians) equal to the relative contact length contains an area equal to A^o, the characteristic area of the encompassing tile t^o, and the area of a vector of oⁱ with the central angle $\mathit{\alpha }={\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{i},\mathrm{o}}$. To ensure the consistency between the assumed circular shape of the patches that compose tⁱ and the cover fractions of the two tiles, the scheme does not use the characteristic area of the outer tile directly, but an approximation ($A_{*}^{\mathrm{o}}$) that is based on the characteristic area of the inner tile (Aⁱ) and on the cover fractions of the outer (c^o) and the inner tile (cⁱ).

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& d^{\mathrm{i},\mathrm{o}}=r^{\mathrm{i}}+{\displaystyle \frac{r^{\mathrm{a}}}{\mathrm{2}}}\mathrm{with}\\ & r^{\mathrm{i}}={\left({\displaystyle \frac{A^{\mathrm{i}}}{\mathit{\pi }}}\right)}^{\mathrm{0.5}}\mathrm{and}\\ & r^{\mathrm{a}}={\left({\displaystyle \frac{\frac{A_{*}^{\mathrm{o}}}{{\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{i},\mathrm{o}}}+A^{\mathrm{i}}}{\mathit{\pi }}}\right)}^{\mathrm{0.5}}\mathrm{and}\\ & A_{*}^{\mathrm{o}}=A^{\mathrm{i}}\cdot \left({\displaystyle \frac{c^{\mathrm{i}}+c^{\mathrm{o}}}{c^{\mathrm{i}}}}-\mathrm{1}\right).\end{array}\end{array}$$

In case of a multilevel, nested connectivity – that is a constellation with the encompassed tile surrounding a third tile – the calculation of l^m,o and d^m,o between the middle (t^m) and the outer tile (t^o), are based on the characteristic area (Aⁱ) of the innermost tile (tⁱ). Here, the current implementation of the scheme only allows for a two-level nesting – i.e. t^m may only contain one tile tⁱ, which itself does not encompass another tile – in which tⁱ is fully surrounded by t^m (${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{i},\mathrm{m}}=\mathrm{1}$). In such cases, l^m,o is calculated by:

$$\begin{array}{}\text{(4)}& l^{\mathrm{m},\mathrm{o}}={\left(A_{*}^{\mathrm{c}}\cdot \mathit{\pi }\cdot \mathrm{4}\right)}^{\mathrm{0.5}}\cdot {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{m},\mathrm{o}},\end{array}$$

where, $A_{*}^{\mathrm{c}}$ constitutes the combined areas of the inner and of the middle tile and is calculated based on the same assumptions used to determine $A_{*}^{\mathrm{o}}$ in Eq. (3).

$$\begin{array}{}\text{(5)}& A_{*}^{\mathrm{c}}=A^{\mathrm{i}}\cdot \left({\displaystyle \frac{c^{\mathrm{i}}+c^{\mathrm{m}}}{c^{\mathrm{i}}}}\right).\end{array}$$

d^m,o is determined analogously to Eq. (3), with the difference being that $A_{*}^{\mathrm{c}}$ is used instead of Aⁱ and that $A_{*}^{\mathrm{o}}$ now accounts for the cover fractions of all three tiles:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& d^{\mathrm{m},\mathrm{o}}={\displaystyle \frac{r^{\mathrm{i}}+r^{\mathrm{a}}}{\mathrm{2}}}\mathrm{with}\\ & r^i={\left({\displaystyle \frac{A^{\mathrm{i}}}{\mathit{\pi }}}\right)}^{\mathrm{0.5}}\mathrm{and}\\ & r^{\mathrm{a}}={\left({\displaystyle \frac{\frac{A_{*}^{\mathrm{o}}}{{\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{m},\mathrm{o}}}+A_{*}^{\mathrm{c}}}{\mathit{\pi }}}\right)}^{\mathrm{0.5}}\mathrm{and}\\ & A_{*}^{\mathrm{o}}=A^{\mathrm{i}}\cdot \left({\displaystyle \frac{c^{\mathrm{i}}+c^{\mathrm{m}}+c^{\mathrm{o}}}{c^{\mathrm{i}}}}-\mathrm{1}\right).\end{array}\end{array}$$

It is possible to use the model for higher levels of nesting, however, in these cases the horizontal relation between two tiles is determined for each pair of tiles individually using Eqs. (2) and (3), which means that a different shape is assumed for each tile depending on whether it constitutes the inner tile or the outer tile in the connection to a given sibling-tile. With the exception of the innermost tile in the hierarchy, this, however results in an underestimation of the interface-length and an overestimation of the center-to-center distance since it assumes the inner tile to have the shape of a circle rather than that of a ring.

Figure 3 Assumed horizontal relationships between connected tiles. The top row shows the most basic case of land-surface heterogeneity, in which the patches of an inner tile are randomly distributed and fully contained within the patches of an outer tile. The second row shows a configuration in which the inner tile is connected to two outer tiles, with different cover fractions but approximately the same connectivity to the inner tile (resulting in the same interface length but different distances between the centers of the tiles). The third row shows a 2-level nesting, where a middle tile contains an inner tile and is connected to two outer tiles with different cover fractions as well as a different connectivity with the middle tile (resulting in roughly the same center-to-center distance with the middle tile). With the patches of the inner tile fully contained by the patches of the middle tile, the assumed relation between these two tiles is the same as between those shown in the top row. The bottom row shows a n-level nesting (n>2), in which the horizontal relationships for each pair are determined based exclusively on the characteristics of the two adjacent tiles.

2.2.2 Vertical relationship

The horizontal moisture fluxes and the associated heat transport between two adjacent tiles also depend on the assumed vertical relation between the two. Here, the most important parameter is the difference in surface elevation (Δh^t,s) between a tile (t) and its connected sibling-tile (s). There are several possibilities to account for difference in surface elevation in the parametrization of the lateral movement of water, with the most straight-forward being the consideration of the geodetic head gradient in the Richards equation. However, since gravity-driven and suction-driven fluxes are being treated separately, T-REX follows a different approach.

Concerning the gravity-driven fluxes, the difference in surface elevation are used to determine the micro-scale slope ($m_{\mathrm{mic}}^{\mathrm{t},\mathrm{s}}$) between two tiles:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& m_{\mathrm{mic}}^{\mathrm{t},\mathrm{s}}={\displaystyle \frac{\mathrm{\Delta }{h}^{\mathrm{t},\mathrm{s}}}{d^{\mathrm{t},\mathrm{s}}}}\mathrm{with}\\ & \mathrm{\Delta }{h}^{\mathrm{t},\mathrm{s}}=h^{\mathrm{t}}-h^{\mathrm{s}}.\end{array}\end{array}$$

The micro-scale slopes, in turn, are used to derive the effective slope of a given tile (m^t) which determines the amount of runoff generated under a given level of saturation of the soil (see Sect. 2.3.3). The effective slope is approximated by the sum of the macro-scale slope ($m_{\mathrm{mac}}^{\mathrm{t}}$), which is an input parameter for the model, and the weighted sum of all sibling-specific micro-scale slopes:

$$\begin{array}{}\text{(8)}& m^{\mathrm{t}}=m_{\mathrm{mac}}^{\mathrm{t}}+\sum _{s=\mathrm{1}}^{\mathrm{nct}}{m}_{\mathrm{mic}}^{\mathrm{t},\mathrm{s}}\cdot {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{t},\mathrm{s}}+m_{\mathrm{mic}}^{\mathrm{s},\mathrm{t}}\cdot {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{s},\mathrm{t}}.\end{array}$$

Here, the sibling-specific slopes are weighted using a modified connectivity ${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{t},\mathrm{s}}$ (${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{s},\mathrm{t}}$) which is based on the interface lengths (l^t,s) and only accounts for those connections with a positive (negative) offset in surface elevation:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{t},\mathrm{s}}=\left\{\begin{array}{ll}\frac{l^{\mathrm{t},\mathrm{s}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{l}^{\mathrm{t},k}+l^{k,\mathrm{t}}}& \mathrm{if}\mathrm{\Delta }{h}^{\mathrm{t},\mathrm{s}}>\mathrm{0}.\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & {\mathrm{\Lambda }}_{{\mathrm{ctc}}^{+}}^{\mathrm{s},\mathrm{t}}=\left\{\begin{array}{ll}\frac{l^{\mathrm{s},\mathrm{t}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{l}^{\mathrm{t},k}+l^{k,\mathrm{t}}}& \mathrm{if}\mathrm{\Delta }{h}^{\mathrm{t},\mathrm{s}}<\mathrm{0}.\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

With respect to the diffusive and conductive fluxes, the lateral exchange is determined for directly interfacing sections of the ground – i.e. sections that have the same absolute vertical position (z^*). Here, the absolute vertical position of any point (p) in the below-ground column of a tile (t) is simply given by its distance to the surface ($z_{\mathrm{p}}^{\mathrm{t}}$), while, for the connected sibling tile (s), the absolute vertical position of a point with the same distance to the surface additionally accounts for the height offset between the two, i.e. $z_{\mathrm{p}}^{*}=z_{\mathrm{p}}^{\mathrm{s}}+\mathrm{\Delta }{h}^{\mathrm{t},\mathrm{s}}$. The model determines the fluxes using a refined grid which encompasses all layer boundaries of the native grid of the tile, all levels of the connected sibling shifted by the height offset and, depending on the process considered, the bedrock boundaries and the depth of standing water at the surface (Fig. 4). As the latter varies with time, the refined grid is determined at the beginning of each time step. This refined grid is also used to map the below-ground runoff on the micro-scale between two connected tiles, while the below-ground runoff on the macro-scale is being transported between two interfacing layers on the tiles' native vertical grid (see Sect. 2.3.4).

Figure 4 Refined vertical grid. Shown is the vertical grid used by the model to determine the diffusive and advective fluxes between two (micro-) connected tiles. This refined grid z^* encompasses all the layers of the vertical grid on the tile z^t and the layers in the vertical grid on the sibling tile z^s shifted by Δh^t,s the height offset between the two. Colors indicate corresponding layers on the tile and its sibling, both on the native and on the refined grid.

2.3 Hydrological fluxes

The lateral hydrological fluxes accounted for in T-REX encompass the diffusive (Sect. 2.3.1) and gravity-driven fluxes of water (Sect. 2.3.2) and the micro-scale redistribution of snow (Sect. 2.3.3), with a concise overview of the key variables used to describe the inter-tile exchange given in Table 3.

Table 3 Key variables in Sect. 2.3 (note that superscripts referring to tiles and subscripts referring to layers have been dropped).

2.3.1 Diffusive exchange of water

The parametrization of suction-driven fluxes (as well as the conductive inter-tile heat exchange) follows the approach outlined in Aas et al. (2019), Nitzbon et al. (2019) and Smith et al. (2022). Here, the fluxes between a tile and a given sibling ($q_{\mathrm{D},i^{*}}^{\mathrm{t},\mathrm{s}}$) are calculated for each layer (i^*) of the refined grid that is located above the bedrock boundary within both tiles (Fig. 3). Since the interfacing soil layers have the same absolute vertical position, the fluxes can simply be calculated using the (geometric) mean hydraulic conductivity ($k_{i^{*}}^{\mathrm{t},\mathrm{s}}$), the difference in matric potential between the two tiles ($\mathrm{\Delta }{\mathit{\psi }}_{i^{*}}^{\mathrm{t},\mathrm{s}}$), the horizontal distance along which the difference in ${\mathit{\psi }}_{i^{*}}^{\mathrm{t},\mathrm{s}}$ occurs (d^t,s) and the interface area ($l^{\mathrm{t},\mathrm{s}}\cdot \mathrm{\Delta }{z}_{i^{*}}$, with $\mathrm{\Delta }{z}_{i^{*}}$ being the thickness of layer i^*):

$$\begin{array}{}\text{(10)}& q_{\mathrm{D},i^{*}}^{\mathrm{t},\mathrm{s}}=l^{\mathrm{t},\mathrm{s}}\cdot \mathrm{\Delta }{z}_{i^{*}}\cdot -k_{i^{*}}^{\mathrm{t},\mathrm{s}}{\displaystyle \frac{\mathrm{\Delta }{\mathit{\psi }}_{i^{*}}^{\mathrm{t},\mathrm{s}}}{d^{\mathrm{t},\mathrm{s}}}}.\end{array}$$

To prevent oscillating fluxes at large time-steps, the lateral exchange is additionally limited to the fluxes that lead to the same degree of saturation in the interfacing soil layers.

The ingress of water ($q_{\mathrm{D},i^{*}}^{\mathrm{d},\mathrm{e}}$) into a section of the soil of a relative elevation (e) that is located within the depth of the surface water body of a relative depression (d) is calculated in the same way, using a matric potential of zero for standing water at the surface. Since the model does not contain any information on the relative location of surface water bodies, it is assumed that these are distributed evenly within the surface area of a given tile. Consequently, the relative contact length l^d,e is reduced according to the maximum inundated fraction of the depression $c_{\mathrm{pnd},\max }^{\mathrm{d}}$.

$$\begin{array}{}\text{(11)}& q_{\mathrm{D},i^{*}}^{\mathrm{d},\mathrm{e}}=c_{\mathrm{pnd},\max }^{\mathrm{d}}\cdot l^{\mathrm{d},\mathrm{e}}\cdot \mathrm{\Delta }{z}_{i^{*}}\cdot k_{i^{*}}^{\mathrm{e}}{\displaystyle \frac{{\mathit{\psi }}_{i^{*}}^{\mathrm{e}}}{d^{\mathrm{d},\mathrm{e}}}}.\end{array}$$

It should be noted that the model does not simulate the lateral egress of water from layers that are located above the surface of the water body as these fluxes are accounted for in the formulation of subsurface runoff (see Sect. 2.3.2).

2.3.2 Surface and subsurface runoff

Runoff generation

In the standard configuration, ICON-Land uses the ARNO-rainfall–runoff model to determine the partitioning of precipitation (q_P), into infiltration (q_I) and runoff (q_R) (Todini, 1996; Dümenil and Todini, 1992; Reick et al., 2021). The scheme has been designed for applications at the watershed-scale – that is horizontal grid spacings > 10 km (Blöschl and Sivapalan, 1995) – and assumes variations in topography to induce a sub-grid-scale soil moisture distribution which, leads to spatially non-uniform infiltration and runoff rates. The use of the ARNO-model becomes problematic if the resolution of the model increases to a point when grid cells no longer represent entire catchments or drainage basins, since the assumptions on the sub-grid-scale soil moisture distribution may no longer be valid. Here, an increase in resolution may not only refer to the actual grid-spacing of the model but also to setups in which the land surface model is applied to represent specific sites, implying a plot-scale resolution of the model. In case that T-REX is applied to tiles representing certain subsystems of a catchment, the ARNO-scheme can also not be used to determine the tile-specific generation of surface and subsurface runoff. For such cases a set of alternative formulations has recently been implemented in the model, which are based on a point-scale approach and make no assumptions about any sub-grid scale variability or spatial extent of the area represented by the model (for more details see Appendix B).

In the point scale approach, surface runoff $q_{\mathrm{R},\mathrm{srf}}^{\mathrm{t}}$ from a given tile (t) occurs if the infiltratable water – that is the sum of rain (or throughfall in case of a vegetation cover), snowmelt, water stored in the surface reservoir and runon from connected tiles – exceeds the soil's infiltration capacity ($q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}$) and the depression storage capacity ($q_{\mathrm{U}}^{\mathrm{t}}$). Here, the model distinguishes between the water that may infiltrate across the entirety of the surface area ($q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}^{\mathrm{t}}$) and the water that may additionally infiltrate below inundated areas ($q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}^{\mathrm{t}}$, with $c_{\mathrm{pnd}}^{\mathrm{t}}$ being the inundated fraction). With respect to the infiltration capacity, the model assumes the latter to be equal to the saturated hydraulic conductivity of the uppermost soil layer, reduced according to the layer's ice content. Runoff also occurs if infiltration results in super saturated soil layers, i.e. the infiltrated water can not percolate downwards fast enough or drain from the soil either at the soil bedrock interface or laterally from the soil layers above the bedrock boundary:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& q_{\mathrm{R},\mathrm{srf}}^{\mathrm{t}}=q_{\mathrm{R},\mathrm{srf},\mathrm{h}}^{\mathrm{t}}+q_{\mathrm{R},\mathrm{srf},\mathrm{d}}^{\mathrm{t}}\mathrm{with}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{h}}^{\mathrm{t}}=q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}^{\mathrm{t}}+q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}^{\mathrm{t}}\mathrm{and}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}^{\mathrm{t}}=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{sl}}^{\mathrm{t}}-q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot}}^{\mathrm{t}}>q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}^{\mathrm{t}}=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}^{\mathrm{t}}-q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}\cdot c_{\mathrm{pnd}}^{\mathrm{t}}-q_{\mathrm{U}}^{\mathrm{t}}& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot}}^{\mathrm{t}}>q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}\cdot c_{\mathrm{pnd}}^{\mathrm{t}}+q_{\mathrm{U}}^{\mathrm{t}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{d}}^{\mathrm{t}}=\sum _{i=\mathrm{1}}^{\mathrm{nsoil}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{upd},i}^{\mathrm{t}}-{\mathit{\theta }}_{\mathrm{sat},i}^{\mathrm{t}}}{\mathrm{\Delta }t}}\mathrm{for}\mathrm{all}{\mathit{\theta }}_{\mathrm{upd},i}^{\mathrm{t}}>{\mathit{\theta }}_{\mathrm{sat},i}^{\mathrm{t}}.\end{array}\end{array}$$

Here, $q_{\mathrm{R},\mathrm{srf},\mathrm{h}}^{\mathrm{t}}$ constitutes the Horton overland flow, with $q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}^{\mathrm{t}}$ referring to the runoff that is generated across the entirety of the tile area and $q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}^{\mathrm{t}}$ to the additional overflow of the surface depression storage. $q_{\mathrm{R},\mathrm{srf},\mathrm{d}}^{\mathrm{t}}$ constitutes the Dunne-type overland flow, with ${\mathit{\theta }}_{\mathrm{upd},i}^{\mathrm{t}}$ being the soil moisture on layer i after the state of the soil has been updated (with infiltration, as the upper boundary condition, limited by $q_{\mathrm{I},\mathrm{cap}}^{\mathrm{t}}$, $q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}^{\mathrm{t}}$ and $q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}^{\mathrm{t}}$) and ${\mathit{\theta }}_{\mathrm{sat},i}^{\mathrm{t}}$ the soil porosity.

The (lateral) subsurface runoff ($q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t}}$) generation on a given layer i of a tile t, is determined as a function of the hydraulic conductivity, $k_i^{\mathrm{t}}$, and the effective slope m^t (see Sect. 2.2):

$$\begin{array}{}\text{(13)}& q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t}}=k_i^{\mathrm{t}}\cdot \sin \left(m^{\mathrm{t}}\cdot {\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right).\end{array}$$

Since high-resolution grid-cells, tiles and even the plots in site-level simulations still cover a given area, the use of the above point-scale formulations is not without problems. Most importantly it can not be assumed that the generated runoff reaches any tributary or connected, lower situated area instantaneously. Consequently, the fluxes are initially stored in intermediary reservoirs which, conceptually, represent the flow paths along which the water moves downslope. With respect to the overland flow we assume that the fluxes are comparatively fast allowing us to neglect any interactions with the underlying ground or the overlying atmosphere. Thus, the outflow ($q_{\mathrm{R},\mathrm{srf}}^{\mathrm{t},*}$), that is the rate at which water reaches the connected downslope areas or the stream system, is exclusively a function of the water content of the intermediary reservoir ($V_{\mathrm{L},\mathrm{srf}}^{\mathrm{t}}$) and the assumed retention time (${\mathit{\tau }}_{\mathrm{srf}}^{\mathrm{t}}$; see below):

$$\begin{array}{}\text{(14)}& q_{\mathrm{R},\mathrm{srf}}^{\mathrm{t},*}={\displaystyle \frac{V_{\mathrm{L},\mathrm{srf}}^{\mathrm{t}}}{{\mathit{\tau }}_{\mathrm{srf}}^{\mathrm{t}}}}.\end{array}$$

The outflow ($q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t},*}$) from a given layer (i) of the below-ground reservoir ($V_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}$) is determined analogously:

$$\begin{array}{}\text{(15)}& q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t},*}={\displaystyle \frac{V_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}}{{\mathit{\tau }}_{\mathrm{blg}}^{\mathrm{t}}}}.\end{array}$$

Here, however, interactions with the surrounding soil can not be neglected. On the one hand, the lateral movement is comparatively slow, allowing the water to percolate downwards and drain beyond the bedrock boundary. On the other hand, the lateral movement is not independent of the level of saturation of the soil, even though a large fraction of the flow may occur through preferential flow networks, partially disconnecting it from the soil matrix. To account for these two effects, the water content of the below-ground reservoir is reduced by a bottom drainage flux ($q_{\mathrm{B},i}^{\mathrm{t}}$) and by a (soil moisture) restoration flow ($q_{X,i}^{\mathrm{t}}$) which moves water from the intermediary reservoir back into the pore spaces, once the water content of a given layer drops below a certain threshold. For the bottom drainage flux we assume that the water drains freely from the reservoir, with the drainage rate being limited to the saturation hydraulic conductivity of (fractured) bedrock,   (k_rock; assumed to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹) and the hydraulic conductivity (k_b) of the soil layer containing the bedrock interface (b). Here, the water is not removed from the soil starting at the bedrock boundary – i.e. the lowest layer of the intermediary reservoir – but from the top of the soil column. This is done to account for the vertical movement within the reservoir, assuming that water percolates downward with the rate it drains into the bedrock. With respect to the restoration flow, we assume that these fluxes occur once the excess water has drained from the soil, setting the moisture threshold to the field capacity (${\mathit{\theta }}_{\mathrm{fc},i}^{\mathrm{t}}$) of a given soil layer.

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}}{\mathrm{\Delta }t}}=q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t}}-q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t},*}-q_{X,i}^{\mathrm{t}}-q_{\mathrm{B},i}^{\mathrm{t}},\mathrm{with}\\ & q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t},*}={\displaystyle \frac{V_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}}{{\mathit{\tau }}_{\mathrm{blg},i}^{\mathrm{t}}}},\mathrm{and}\\ & q_{X,i}^{\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{0},& \mathrm{if}{\mathit{\theta }}_{\mathrm{fc},i}^{\mathrm{t}}>{\mathit{\theta }}_i^{\mathrm{t}}\\ \frac{{\mathit{\theta }}_{\mathrm{fc},i}^{\mathrm{t}}-{\mathit{\theta }}_i^{\mathrm{t}}}{\mathrm{\Delta }t},& \mathrm{if}{\mathit{\theta }}_{\mathrm{fc},i}^{\mathrm{t}}<{\mathit{\theta }}_i^{\mathrm{t}}\mathrm{and}{\mathit{\theta }}_{\mathrm{fc},i}^{\mathrm{t}}-{\mathit{\theta }}_i^{\mathrm{t}}<V_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}\\ \frac{V_{\mathrm{L},\mathrm{blg},i}^{\mathrm{t}}}{\mathrm{\Delta }t}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & q_{\mathrm{B},i}^{\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{0}.& \mathrm{if}\frac{\sum _{l=\mathrm{1}}^{i-\mathrm{1}}{V}_{\mathrm{L},\mathrm{blg},l}^{\mathrm{t}}}{\mathrm{\Delta }t}>k_{\mathrm{bot}}\\ \frac{V_{\mathrm{L},\mathrm{blg},i}^t}{\mathrm{\Delta }t}& \mathrm{if}\frac{\sum _{l=\mathrm{1}}^i{V}_{\mathrm{L},\mathrm{blg},l}^{\mathrm{t}}}{\mathrm{\Delta }t}<k_{\mathrm{bot}}\\ k_{\mathrm{bot}}-\frac{\sum _{l=\mathrm{1}}^{l-\mathrm{1}}{V}_{\mathrm{L},\mathrm{blg},l}^{\mathrm{t}}}{\mathrm{\Delta }t}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & k_{\mathrm{bot}}=\left\{\begin{array}{ll}{k}_{\mathrm{rock}}& \mathrm{if}{k}_{\mathrm{rock}}<k_{\mathrm{b}}\\ k_{\mathrm{b}}& \mathrm{otherwise}\end{array}\right.\end{array}\end{array}$$

Retention time and partitioning of runoff fluxes between sibling tiles

While the generation of gravity-driven fluxes, that is surface- and subsurface runoff, depends mainly on the characteristics of a given tile (t), i.e. the slope and soil characteristics, the run-on fluxes are largely determined by the relationships with the connected sibling-tiles (s). Here, the model employs lag factors to account for the fact that some time is required for the generated runoff – which is determined based on point-scale parametrizations that assume no spatial extent of a tile – to reach the connected sibling-tiles. With respect to micro-scale connections, the lag factors ($f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}$, where l either stands for surface [srf] or below-ground [blg]) are determined based on the time-step length (Δt), the lateral flow velocities (v_l) and the center-to-center distance of two connected tiles (d^t,s):

$$\begin{array}{}\text{(17)}& f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}={\displaystyle \frac{\mathrm{\Delta }t\cdot d^{\mathrm{t},\mathrm{s}}}{v_{\mathrm{l}}}}.\end{array}$$

It is difficult to determine appropriate lateral flow velocities since this involves a number of highly uncertain assumptions. With respect to the lateral velocities at the surface (v_srf), it is conceivable to parameterize the fluxes assuming an open channel flow with a given slope and surface roughness which, however, requires an assumption on the distribution and geometry of the emerging flow paths at the surface, such as rivulets. With respect to the below-ground fluxes (v_blg) it appears obvious to employ the saturation hydraulic conductivity also in combination with the local slope, assuming that the respective flow is adequately described by Darcy's law. This however, neglects the possibility that, even in areas with generally less permeable soils, preferential flow, such as funnel-, pipe- and bypass flow, allows for high transmissivties – i.e. through regions with more permeable soils, cracks or crevices and macro-pores (Barcelo and Nieber, 1982; McDonnell, 1990; Uchida et al., 2001; Di Prima et al., 2018). The latter may be accounted for in the estimate of the transmissivty via a macropore enhancement of the hydraulic conductivity (Milly et al., 2014). However, given the uncertainty of the factors involved and the simple nature of our approach, we instead use globally uniform values for the flow velocities, assuming a value of 1 m h⁻¹ for v_blg, with the order of magnitude corresponding to the conductivity of soils including the effects of macropores and pipes or the saturated hydraulic conductivity of sand. Here, the values are independent of the degree of saturation of the soil, meaning there is no distinction between groundwater- and interflow. This was done because the model strongly inhibits the lateral fluxes in the vadose zone, since the water in the intermediary reservoir migrates back into the soil column of the runoff-generating tile once the latter starts drying and soil moisture levels drop below the field capacity. For the velocity of the overland flow, we simply assume that it is an order of magnitude larger than the below-ground values, setting v_srf to 10 m h⁻¹, which is at the lower end of the range of previous estimates (McCaig, 1983).

An accurate description of the lateral fluxes on the macro-scale is even more difficult as it involves additional uncertainty. On the catchment- or hill-slope-scale, the distances are usually large enough to allow water to percolate downwards to the bedrock boundary, where horizontal fluxes occur as sheet flow or through preferential flow networks, possibly returning to the surface at the bottom of the slope (Kirkby, 1988). Approaches that use Darcy's law to describe the long-distance lateral transport across sloped terrain risk underestimating the subsurface runoff, with observed lateral flow velocities sometimes being orders of magnitude larger than what can be expected from a flow through a porous medium (Graham et al., 2010). Since our model does not include a representation of the respective dynamics, these fluxes need to be parameterized. At least for the surface runoff, it is possible to estimate specific lag factors – analogously to the determination of unit-hydrographs – based on the geomorphological characteristics of the grid-cell, while the velocity of sub-surface runoff may be approximated based on observed rainfall and base flow rates (Lohmann et al., 1996; Kumar et al., 2007; Lehner et al., 2008; Singh et al., 2014; Mizukami et al., 2016). For the present implementation, however, we use a simpler approach and calculate the macro-scale lag factors $f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}$ using the same assumptions as the ARNO-scheme, namely that the drainage of excess water from any catchment-sized area takes place over a period of a few days (Todini, 1996; Dümenil and Todini, 1992). In the scheme's implementation in ICON-Land, the minimum and maximum drainage rates are globally uniform and independent of the model resolution, which further assumes that the overall retention time is largely determined by the time it takes water to pass through the ground and that the distances to the nearest tributary are not only similar in all drainage basins but also small relative to the horizontal grid spacing of the model. Accordingly, we determine the lag factors ($f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}$) based on the time-step length and globally uniform retention times (${\mathit{\tau }}_{\mathrm{ret},\mathrm{mac},\mathrm{l}}$):

$$\begin{array}{}\text{(18)}& f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}={\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\tau }}_{\mathrm{ret},\mathrm{mac},\mathrm{l}}}}.\end{array}$$

For standard applications, in which the macro-scale transport moves water across typical drainage-basin scales, ${\mathit{\tau }}_{\mathrm{mac},\mathrm{blg}}^{\mathrm{ret}}$ is set to 120 h and ${\mathit{\tau }}_{\mathrm{mac},\mathrm{srf}}^{\mathrm{ret}}$ is set to 10 h – again assuming that the flow velocities at the surface are a magnitude larger than those below ground. In case that the distances represented by the macro-scale connections are below the catchment scale, ${\mathit{\tau }}_{\mathrm{mac},\mathrm{blg}}^{\mathrm{ret}}$ and ${\mathit{\tau }}_{\mathrm{mac},\mathrm{srf}}^{\mathrm{ret}}$ can be specified via a namelist parameter. However, for T-REX applications that implicitly resolve slopes by a number of tiles (e.g. corresponding to different height bands), it may be more appropriate not do so using macro-scale connections. In such cases, the relations between tiles should rather be based on micro-scale connectivities where distances are explicitly accounted for in the estimation of the lag factors.

The weighted sum of the sibling-specific lag factors, constitutes the average outflow lag factors of each tile ($f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t}}$ and $f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t}}$). For macro-scale connections, the additional assumption is being made that for those fractions of the tile that feature micro-scale connections, macro-scale fluxes may not occur. Here, the average lag factors are scaled according to the cover fraction of the tile and the cover fractions of the (micro-) connected sibling-tiles:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t}}=\sum _{s=\mathrm{1}}^{\mathrm{nct}}{f}_{\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}\cdot {\mathrm{\Lambda }}_{{\mathrm{hfp}}^{*},\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}\mathrm{and}\\ & f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t}}=\sum _{s=\mathrm{1}}^{\mathrm{nct}}{f}_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}\cdot {\mathrm{\Lambda }}_{{\mathrm{hfp}}^{*},\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}\mathrm{with}\\ & {\mathrm{\Lambda }}_{{\mathrm{hfp}}^{*},\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}={\displaystyle \frac{c^{\mathrm{s}}\cdot {\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{c}^k\cdot {\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}}}\mathrm{and}\\ & {\mathrm{\Lambda }}_{{\mathrm{hfp}}^{*},\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}=\left\{\begin{array}{ll}\frac{c^{\mathrm{s}}\cdot {\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{c}^k\cdot {\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mac},\mathrm{l}}^{\mathrm{t},k}}\cdot \frac{c^{\mathrm{t}}-\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}\cdot c^k}{c^{\mathrm{t}}}& \mathrm{if}{c}^{\mathrm{t}}-\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}\cdot c^k>\mathrm{0}.\\ \mathrm{0}.& \mathrm{otherwise}.\end{array}\right.\end{array}\end{array}$$

The average retention times for each tile (${\mathit{\tau }}_{\mathrm{l}}^{\mathrm{t}}$, where l either stands for surface [srf] or below-ground [blg]), are calculated based on the average lag factors for surface and below-ground fluxes ($f_{\mathrm{l}}^{\mathrm{t}}$), which can be obtained simply by adding the respective lag factors for fluxes on the micro- and macro scale:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathit{\tau }}_{\mathrm{l}}^{\mathrm{t}}={\displaystyle \frac{\mathrm{\Delta }t}{f_{\mathrm{l}}^{\mathrm{t}}}}\mathrm{with}\\ & f_{\mathrm{l}}^{\mathrm{t}}=f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t}}+f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t}}.\end{array}\end{array}$$

Finally, the fraction of the surface- and subsurface runoff ($p_{\mathrm{l}}^{\mathrm{t},\mathrm{s}}$) from tile (t) received by any of its siblings (s) is determined by the ratio of the sibling-specific lag factor and the average lag factor. For tiles that are not fully connected to downstream tiles – i.e. $\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}+{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mac},\mathrm{l}}^{\mathrm{t},k}<\mathrm{1}$ – the fractions of outflow received by the tile's connected siblings are reduced accordingly, with the residual outflow contributing to the stream-flow directly.

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& p_{\mathrm{l}}^{\mathrm{t},\mathrm{s}}=\\ & \left\{\begin{array}{ll}\frac{f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}+f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}}{f_{\mathrm{l}}^{\mathrm{t}}}\cdot \sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}+{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mac},\mathrm{l}}^{\mathrm{t},k}& \mathrm{if}\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mic},\mathrm{l}}^{\mathrm{t},k}+{\mathrm{\Lambda }}_{\mathrm{hfp},\mathrm{mac},\mathrm{l}}^{\mathrm{t},k}<\mathrm{1}\\ \frac{f_{\mathrm{mic},\mathrm{l}}^{\mathrm{t},\mathrm{s}}+f_{\mathrm{mac},\mathrm{l}}^{\mathrm{t},\mathrm{s}}}{f_{\mathrm{l}}^{\mathrm{t}}}& \mathrm{otherwise}\end{array}\right.\end{array}\end{array}$$

Applying this partitioning to the fluxes from the intermediary reservoirs ($q_{\mathrm{R},\mathrm{l}}^{\mathrm{t},*}$), the lateral runoff flux ($q_{\mathrm{R},\mathrm{l}}^{\mathrm{t},\mathrm{s},*}$) from tile t to any of its siblings s is calculated according to:

$$\begin{array}{}\text{(22)}& q_{\mathrm{R},\mathrm{l}}^{\mathrm{t},\mathrm{s},*}=q_{\mathrm{R},\mathrm{l}}^{\mathrm{t},*}\cdot p_{\mathrm{l}}^{\mathrm{t},\mathrm{s}}.\end{array}$$

In case of overland flow, the runoff contributes to the infiltratable water of the receiving sibling tile ($q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}^{\mathrm{s}}$ and $q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}^{\mathrm{s}}$), with the inflow primarily increasing the water content of the surface depression storage ($\mathrm{\Delta }{V}_{\mathrm{sfc}}^{\mathrm{s}}$):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}^{\mathrm{s}}=q_{\mathrm{P}}\cdot \left(\mathrm{1}-c_{\mathrm{pnd},\max }^{\mathrm{s}}{}^{\frac{\mathrm{1}}{\mathrm{3}}}\right)+q_{\mathrm{R},\mathrm{srf}}^{\mathrm{t},\mathrm{s},*}-{\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{sfc}}^{\mathrm{s}}}{\mathrm{\Delta }t}}\mathrm{and}\\ & q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}^{\mathrm{s}}=q_{\mathrm{P}}\cdot c_{\mathrm{pnd},\max }^{\mathrm{s}}{}^{\frac{\mathrm{1}}{\mathrm{3}}}+{\displaystyle \frac{V_{\mathrm{sfc}}^{\mathrm{s}}+\mathrm{\Delta }{V}_{\mathrm{sfc}}^{\mathrm{s}}}{\mathrm{\Delta }t}}\mathrm{with}\\ & \mathrm{\Delta }{V}_{\mathrm{sfc}}^{\mathrm{s}}=\left\{\begin{array}{ll}{q}_{\mathrm{R},\mathrm{srf}}^{\mathrm{t},\mathrm{s},*}\cdot \mathrm{\Delta }t& \mathrm{if}{q}_{\mathrm{R},\mathrm{srf}}^{\mathrm{t},\mathrm{s},*}\cdot \mathrm{\Delta }t<V_{\mathrm{sfc},\max }^{\mathrm{s}}-V_{\mathrm{sfc}}^{\mathrm{s}}\\ V_{\mathrm{sfc},\max }^{\mathrm{s}}-V_{\mathrm{sfc}}^{\mathrm{s}}& \mathrm{otherwise}.\end{array}\right.\end{array}\end{array}$$

Here, q_P constitutes precipitation and $c_{\mathrm{pnd},\max }^{\mathrm{s}}$ the maximum pond fraction of the sibling tile. The subsurface runoff from a given layer on a tile contributes to the soil moisture on the corresponding layer of the sibling, with the height offset between the two accounted for in case of micro-scale connections:

$$\begin{array}{}\text{(24)}& {\mathit{\theta }}_{\mathrm{upd},i^{*}}^{\mathrm{s}}={\mathit{\theta }}_{i^{*}}^{\mathrm{s}}+q_{\mathrm{R},\mathrm{blg},i^{*}}^{\mathrm{t},\mathrm{s},*}\cdot \mathrm{\Delta }t.\end{array}$$

In case of macro-scale connections, the height offset is being ignored to retain consistency with our assumptions on the vertical movement of water with the soil column, in particular allowing the water to flow along the bedrock boundary within the grid cell:

$$\begin{array}{}\text{(25)}& {\mathit{\theta }}_{\mathrm{upd},i}^{\mathrm{s}}={\mathit{\theta }}_i^{\mathrm{s}}+q_{\mathrm{R},\mathrm{blg},i}^{\mathrm{t},\mathrm{s}*}\cdot \mathrm{\Delta }t.\end{array}$$

2.3.3 Snow redistribution

The inter-tile redistribution of snow also follows Aas et al. (2019), Nitzbon et al. (2019) and Smith et al. (2022). The approach is based on the assumption that wind-blown snow predominantly acts to level the surface of the snowpack between two adjacent areas with different surface elevations. The scheme further assumes that drifting mainly affects the freshly deposited, powdery snow, rather than the settled, more granular snow that has begun metamorphosis. Thus, while the approach neglects effects of abrasion over longer periods, it allows to determine the redistribution of snow simply as a change in the snow deposition flux ($\mathrm{\Delta }{q}_{\mathrm{S}}^{\mathrm{h},\mathrm{l}}$), rather than a modification of the existing snowpack. Here, a redistribution may only occur if the thickness of the existing snowpack exceeds the height for which it can be assumed that drifting would be inhibited by the existing vegetation. For this threshold a snow water equivalent of 0.01 m was assumed, which corresponds to a snow depth of between roughly 3 and 20 cm depending on the simulated snow densities. If the snow water equivalent of the snowpack is above this value, the hift in snow deposition flux from the tile with the higher (h) to the tile with the lower snow surface (l) is determined as a function of the relative contact lengths (${\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{h},\mathrm{l}}$), the initial snow deposition flux (q_S) and the initial difference in surface elevation of the adjacent snow surfaces ($\mathrm{\Delta }{h}_{\mathrm{snw}}^{\mathrm{h},\mathrm{l}}$):

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& \mathrm{\Delta }{q}_{\mathrm{S}}^{\mathrm{h},\mathrm{l}}={\displaystyle \frac{{\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{h},\mathrm{l}}}{\sum _{k=\mathrm{1}}^{\mathrm{nct}}{\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{\mathrm{h},k}+{\mathrm{\Lambda }}_{{\mathrm{ctc}}^{*}}^{k,\mathrm{h}}}}\cdot \left\{\begin{array}{ll}{q}_{\mathrm{S}}\left(\frac{\mathrm{1}}{c_{\mathrm{rel}}^{\mathrm{l}}}-\mathrm{1}\right)& \mathrm{if}{q}_{\mathrm{S}}\le \frac{\mathrm{\Delta }{h}_{\mathrm{snw}}^{\mathrm{h},\mathrm{l}}}{\mathrm{\Delta }t}\frac{{\mathit{\rho }}_{\mathrm{snw}}}{{\mathit{\rho }}_{\mathrm{w}}}{c}_{\mathrm{rel}}^{\mathrm{l}}\\ \frac{\mathrm{\Delta }{h}_{\mathrm{snw}}^{\mathrm{h},\mathrm{l}}}{\mathrm{\Delta }t}\frac{{\mathit{\rho }}_{\mathrm{snw}}}{{\mathit{\rho }}_{\mathrm{w}}}\left(\mathrm{1}-c_{\mathrm{rel}}^{\mathrm{l}}\right)& \mathrm{otherwise}\end{array}\right.\mathrm{with}\\ & c_{\mathrm{rel}}^{\mathrm{l}}={\displaystyle \frac{c_{\mathrm{l}}}{c_{\mathrm{l}}+c_{\mathrm{h}}}}\end{array}\end{array}$$

where ρ_w the density of liquid water, ρ_snw the density of snow and c^l and c^h the cover fractions of the adjacent tiles.

This is a first-order simplification that neglects key factors such as wind direction and speed, which are critical drivers of real-world snow transport. This is mainly because of the structural design of the land-surface component, where subgrid tiles coexist statistically within a grid cell but are not explicitly spatially arranged. Consequently, the model cannot determine whether wind preferentially connects specific tile pairs or how fluxes would vary with wind direction. As a result, the present implementation cannot be used for the snow transport on the macro-scale but only for micro-scale heterogeneity where the relative arrangement of tiles is sufficiently random for snow transport to become approximately independent of wind direction (see Discussion). For such conditions however, this elevation-based redistribution can at least account for the preferential deposition of snow in microscale depressions (Clark et al., 2011).

2.4 Heat exchange

Concerning the lateral heat fluxes, T-REX distinguishes between conductive heat transport (Sect. 2.4.1) – the heat transfer through direct molecular diffusion – and advective heat transport (Sect. 2.4.2) – the heat transfer via bulk water movement – with a concise overview of the key variables used to describe the inter-tile exchange given in Table 4.

Table 4 Key variables in Sect. 2.4 (note that superscripts referring to tiles and subscripts referring to layers have been dropped).

2.4.1 Conductive heat exchange

The conductive heat transport (${\mathit{\varphi }}_{\mathrm{C},i^{*}}^{\mathrm{t},\mathrm{s}}$) follows the same approach as the diffusive transport of water:

$$\begin{array}{}\text{(27)}& {\mathit{\varphi }}_{\mathrm{C},i^{*}}^{\mathrm{t},\mathrm{s}}=\mathrm{\Delta }{z}_{i^{*}}\cdot l^{\mathrm{t},\mathrm{s}}\cdot -{\mathit{\lambda }}_{i^{*}}^{\mathrm{t},\mathrm{s}}{\displaystyle \frac{\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}}{d^{\mathrm{t},\mathrm{s}}}},\end{array}$$

with ${\mathit{\lambda }}_{i^{*}}^{\mathrm{t},\mathrm{s}}$ being the average soil heat conductivity and $\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}$ the difference in soil temperature between the adjacent tiles. To prevent oscillating fluxes at large time-steps, the lateral exchange is additionally limited to the fluxes that lead to the same soil temperature in the two interfacing soil layers. The calculations are performed on, essentially, the same refined vertical grid with the difference being that the conductive heat exchange is also calculated for layers (i^*) that are located below the bedrock boundary.

2.4.2 Advective heat fluxes

ICON-Land does not simulate the temperatures of water in the soil explicitly. Instead the below-ground temperatures are determined for the soil-water compound – by employing a soil heat capacity and conductivity that is determined based on the water and ice content of a given layer – which implies that all water has the same temperature as the surrounding soil matrix. Based on this assumption, the advective heat transport on a layer of the refined grid (i^*) between a tile (t) and its connected siblings (s) can be determined simply as a function of the volume flux of water ($q_{\mathrm{W},i^{*}}^{\mathrm{t},\mathrm{s}}$) from the tile to a given sibling and the difference in soil temperature ($\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}$) between the two:

$$\begin{array}{}\text{(28)}& {\mathit{\varphi }}_{\mathrm{A},i^{*}}^{\mathrm{t},\mathrm{s}}=\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}\cdot {\mathrm{\Theta }}_{\mathrm{w}}\cdot q_{\mathrm{W},i^{*}}^{\mathrm{t},\mathrm{s}},\end{array}$$

with Θ_w being the volumetric heat capacity of water and $q_{\mathrm{W},i^{*}}^{\mathrm{t},\mathrm{s}}$ the sum of the subsurface runoff and diffusive fluxes from the tile to a given sibling.

Similarly, standing water at the surface does not have an explicit temperature, with the exception being large lakes, but these are represented by separate tiles that do not interact with the surrounding land surface. Furthermore, the heat storage of (all other) surface water bodies, is not accounted for in the surface energy balance of the model. This makes it impossible to simulate the advective heat fluxes due to surface runoff, with the latter initially being given to the surface water reservoir of the receiving tile before it infiltrates into the ground or runs off. Here an optional formulation was introduced into the model which modifies the surface heat capacity and is based on the assumptions that, (i) any precipitation has the same temperature as the surface (which is the same assumption as in the standard model), (ii) that the surface water bodies have no thermal stratification – which is a permissible simplification as long as the surface water bodies are limited to puddles, ponds and politic lakes – and (iii) that the heat exchange between surface water and the ground is sufficiently fast for the uppermost soil layer to have approximately the same temperature as the overlying surface water body. This modified volumetric heat capacity of the top soil layer is calculated as:

$$\begin{array}{}\text{(29)}& {\mathrm{\Theta }}_{\mathrm{soil},\mathrm{1}}^{*}={\displaystyle \frac{{\mathrm{\Theta }}_{\mathrm{soil},\mathrm{1}}\cdot \mathrm{\Delta }{z}_{\mathrm{1}}+{\mathrm{\Theta }}_{\mathrm{w}}\cdot h_{\mathrm{wtr}}+{\mathrm{\Theta }}_{\mathrm{ice}}\cdot h_{\mathrm{ice}}}{\mathrm{\Delta }{z}_{\mathrm{1}}}}.\end{array}$$

Using the modified heat capacity, any change in temperature due to a surface advective heat flux calculated by Eq. (28) (with $\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}$ referring to the uppermost soil layer of both tiles and $q_{\mathrm{W},i^{*}}^{\mathrm{t},\mathrm{s}}$ being the surface runoff) is (implicitly) being applied to the surface water body along with the uppermost soil layer.

With large time steps and lateral fluxes of water, it is possible that the lateral heat fluxes calculated above could increase (lower) the temperature on the receiving tile above (below) those of the tile the water originated on. To prevent this, the lateral heat fluxes are limited by the temperature differences and the heat capacity of the receiving tile:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\mathit{\varphi }}_{\mathrm{A},i^{*}}^{\mathrm{t},\mathrm{s}}=& \mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}\cdot {\displaystyle \frac{{\mathrm{\Theta }}_{\mathrm{soil},i^{*}}^{\mathrm{s}}\cdot \mathrm{\Delta }{z}_{i^{*}}}{\mathrm{\Delta }t}}\mathrm{if}|{\mathit{\varphi }}_{\mathrm{A},i^{*}}^{\mathrm{t},\mathrm{s}}\\ & |>|\mathrm{\Delta }{T}_{i^{*}}^{\mathrm{t},\mathrm{s}}\cdot {\displaystyle \frac{{\mathrm{\Theta }}_{\mathrm{soil},i^{*}}^{\mathrm{s}}\cdot \mathrm{\Delta }{z}_{i^{*}}}{\mathrm{\Delta }t}}|.\end{array}\end{array}$$

3 Example applications

3.1 Lateral water transport

For the first example application, we aim to investigate the effect of micro- and macro-scale lateral transport processes on the land's capacity to retain water and the resulting impacts on the turbulent land-atmosphere exchange and surface temperatures. On the macro scale, we separate the grid cell into those parts that generate runoff, the uplands, and those that receive the generated runoff, i.e. the lowlands. With respect to the micro scale, the setup distinguishes between those grid-cell areas where water (including snow) may pool on top of the surface, i.e. local depressions, and those where excess water runs of immediately, i.e. relative elevations. Here, the logic of the tile definition provides the following characteristic connectivities, assuming that the depressions, both on the uplands and lowlands (DU, DL), are fully surrounded by and receive the runoff from relative elevations (EU, EL) and that the surface runoff from the uplands occours via surface channels that are primarily connected to lowland depressions, while the below-ground runoff similarly constitutes inflow to the lowland elevations and depressions (with the actual connectivity being determined by the cover fractions of the two lowland tiles):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Lambda }}_{\mathrm{ctc}}=\begin{array}{ccccc}& \mathrm{EU}& \mathrm{DU}& \mathrm{EL}& \mathrm{DL}\\ \mathrm{EU}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{DU}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{EL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{DL}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\\ {\displaystyle }& {\displaystyle }{\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{micro},\mathrm{srf}}=\begin{array}{ccccc}& \mathrm{EU}& \mathrm{DU}& \mathrm{EL}& \mathrm{DL}\\ \mathrm{EU}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{DU}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{EL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{DL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ {\displaystyle }& {\displaystyle }{\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{micro},\mathrm{blg}}=\begin{array}{ccccc}& \mathrm{EU}& \mathrm{DU}& \mathrm{EL}& \mathrm{DL}\\ \mathrm{EU}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{DU}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{EL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{DL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ {\displaystyle }& {\displaystyle }{\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{macro},\mathrm{srf}}=\begin{array}{ccccc}& \mathrm{EU}& \mathrm{DU}& \mathrm{EL}& \mathrm{DL}\\ \mathrm{EU}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{DU}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{EL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{DL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ {\displaystyle }& {\displaystyle }{\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{macro},\mathrm{blg}}=\begin{array}{ccccc}& \mathrm{EU}& \mathrm{DU}& \mathrm{EL}& \mathrm{DL}\\ \mathrm{EU}& \mathrm{0}& \mathrm{0}& \mathrm{0.5}& \mathrm{0.5}\\ \mathrm{DU}& \mathrm{0}& \mathrm{0}& \mathrm{0.5}& \mathrm{0.5}\\ \mathrm{EL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{DL}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}.\end{array}$$

The respective fractions are determined by combining the 30 m resolution version of the Copernicus DEM (European Space Agency, 2024) together with the CEH CTI dataset (Marthews et al., 2015a) at 15^′′ resolution as follows: In a first step, we use the Whitebox Tools Open Core (Lindsay, 2014) to fill all local depressions in the DEM data. The difference between this and the original data provides us with the depression depth. Next, we aggregate this data to a 15^′′ resolution thereby also computing the depression fraction as grid-cell area belonging to local depressions relative to the total cell area. Furthermore, we use the Connected-Components implementation of the OpenCV Toolbox (Bradski, 2000) to label all grid cells belonging to the same depression. Thus, we can divide the overall depression area by the number of different depressions to gain an estimate of the mean depression size for each 15^′′ cell. These depression characteristics are then cross-referenced with the CTI data. Wherever the CTI ≤ 5.7 we consider the cell to be an upland and therefore assign the depression fraction, depth and average size to the upland region and vice versa for lowlands with a CTI > 5.7. Here, we chose a threshold value of 5.7, because it is the median of the global dataset and because it roughly corresponds to the value commonly used separate lowlands from uplands and slopes (Marthews et al., 2015b). Finally, the data is further aggregated onto the R2B4 grid (≈160 km resolution) which is used for our global simulations. Grid cells with missing values are set to the average of their respective quantities. We conduct ICON-Land-standalone simulations driven by climate forcing from the Global Soil Wetness Project Phase 3 (GSWP3, Dirmeyer et al., 2006; Kim, 2017). The simulations were performed for the period 1979–2018, with the analysis using the model output from the period 2009–2018. In the following, we compare simulations in which all subgrid-scale lateral transport of water is disabled to simulation in which either micro-scale or macro-scale transport processes are accounted for as well as a simulation in which transport processes on both scales are enabled in the model.

Water is transported laterally mainly in those regions of the world where the water availability at the land surface exceeds the atmospheric moisture demand, at least temporarily. Thus, the model simulates only small lateral subgrid-scale fluxes in the arid and semi-arid regions, while they can exceed 500 mm yr⁻¹ in the humid areas (Fig. 5a). However, also the humid regions in the tropics feature extensive areas with comparatively little lateral water movement, such as the Amazon basin. Here, it is most often the topography that does not favor large quantities of water being moved laterally through the soil or at the surface. Comparatively small slopes in the contributing areas generate only little surface- and lateral subsurface runoff and, in the model, most of the water drains at the bottom of the soil column, contributing to the channel flow as base-flow below the bedrock boundary (not shown). In contrast, the flat areas in the cold regions often feature larger lateral fluxes, since here the ground is predominantly frozen, at least during the snowmelt season, and water can not percolate deep into the soil but is forced to run off laterally within a comparatively shallow layer close to the surface. In case of the micro-scale fluxes another minor effect stems from the redistribution of snow, which already shifts water from relative elevations to the depressions during the cold season (see Sect. 4).

Figure 5 Effects of the lateral water transport on terrestrial water storage, latent heat flux and subgrid-scale temperature variability. Shown is the total subgrid-scale lateral flux as the sum of surface runoff, lateral drainage, diffusive fluxes, and snow redistribution for a simulation in which micro- and macro-scale transport processes are enabled (a), a simulation in which only the micro-scale transport is active (b) and a simulation that only accounts for the macro-scale lateral fluxes (c; one these scales diffusive fluxes and snow redistribution are not simulated). Panels (d)–(f) show the resulting impacts on the total terrestrial water storage. These were calculated as the differences between a given simulation with (parts of) the lateral transport activated and a reference simulation in which the lateral transport is switched off. Panels (g)–(i) show the impact on the grid-cell mean latent heat flux, while panels (j)–(l) show the change in the maximum spatial (subgrid-scale) temperature variability. The later was calculated as the temperature difference between the warmest and coldest tile within a given grid cell.

Lateral transport processes enable low-lying parts of the gridcell to store the runoff from upslope areas and allow a redistribution of water between those areas with a large depression storage and those that may not retain water at the surface. Thus, they almost exclusively increases the terrestrial water storage (Fig. 5d), with marked increases in the total soil water content and in the wetland area (not shown). Additionally, snow redistribution further increases total water storage, as the preferential accumulation of snow in small-scale depressions delays the moment when the grid cell becomes completely snow-free. There are, however, a few grid cells, mainly in the (sub-) polar regions, where the total terrestrial water, in particular the soil water content, is reduced. Here, the lateral transport initially increases the moisture content in the near-surface layers after snowmelt, which increases the heat conductivity of the ground (not shown). This, in turn, raises the ground heat flux and allows a higher fraction of the soil ice to melt. And with more water being liquid, the resulting increase in evapotranspiration and drainage during summer and fall causes an overall reduction in the annual mean soil water content. Substantial increases in the terrestrial water storage are mainly limited to the temperate zone and the (sub) polar region, where the additional retention capacity allows to partially compensate the differences in the seasonality of water availability and (atmospheric) moisture demand. Depressions, in particular, can store parts of the spring snowmelt, with the water subsequently diffusing into the surrounding areas increasing the water availability during spring and summer when potential evaporation rates are high. Here, the effect on the simulated wetland extent is especially pronounced and increases in the inundated area – relative to the potentially inundated area – of up to 70 % across large parts of the temperate and most of the polar and subpolar regions show that lateral exchange processes are essential for maintaining standing water at the surface (not shown). Pronounced effects on the total soil water storage are, however, limited to those regions that feature a deep bedrock boundary, where increases in total soil water of up to 0.2 m often correspond to a relative increase in the soil water content of about 20 % (not shown). In regions where the bedrock boundary is comparatively shallow, such as most of eastern Siberia, there are only small increases in the soil water content, despite large lateral fluxes. In contrast, some large effects are observable in regions where the lateral fluxes are comparatively small, such as the West Siberian Lowlands, where a large depression capacity, deep soils and low potential evaporation rates facilitate the water storage over longer periods allowing even comparatively small lateral fluxes to have a notable impact on the soil moisture and depression storage.

Within the subtropics and the tropics, the effects on the terrestrial water storage are largely negligible (Fig. 5d). Here, most precipitation transpires and evaporates locally and the model simulates only minor lateral fluxes. Additionally, the temporal offset between water availability and demand is much smaller than in the temperate and polar region and most of the water that is redistributed laterally evaporates or is transpired swiftly rather than being stored in the soil or in depressions over longer periods. Consequently, there are areas in the tropics and subtropics with only a minor impact on the terrestrial water storage but with a notable effect on the latent heat flux (Fig. 5g). In general, the inclusion of the lateral exchange fluxes increases the latent heat flux notably, mainly in those areas that also feature large changes in the terrestrial water storage where respective effects can correspond to relative increases of up to 10 %. And since the increase in latent heat flux is mainly balanced by a decrease in the sensible heat flux, the Bowen ratio changes by as much as 20 %. Yet even a 20-percent change in the Bowen ratio does not drastically alter the overall surface energy balance on the grid-cell level. Consequently, the effects on the (grid cell) average surface temperatures are comparatively small and the cooling due to the inclusion of the lateral water transport rarely exceeds −0.2 K (not shown). However, the spatial subgrid-scale temperature variability is extremely sensitive to the inclusion of the lateral transport processes, since the local effects on the Bowen ratio in the lowlands and depressions can far exceed the grid-cell mean value. Here, the temperature differences between the warmest and the coldest grid-cell fraction can increase by as much as 3 K (Fig. 5j), corresponding to an order of magnitude increase across most of North America and Eastern Siberia.

A comparison of the micro- (Fig. 5b) and the macro-scale fluxes (Fig. 5c) reveals a large similarity with respect to the spatial patterns and there are only few areas in which the model simulates relatively large micro-scale fluxes but only relatively small macro-scale fluxes (and vice versa). This suggests that the global patterns are mainly controlled by scale-independent factors, such as the state of the atmosphere – which determines precipitation and potential evaporation but also the general thermal state of the soil (frozen or unfrozen) – and the soil properties, and less by the specific topography within a given grid cell. In contrast the overall magnitude is quite different, with the lateral fluxes being predominantly larger for the micro-scale exchange – on average by about 40 %. With respect to the total terrestrial water storage (Fig. 5e and f), the effect of the micro-scale transport is about 25 % larger than that of the macro-scale transport, with the impact on the total inundated area being almost 50 % larger. In contrast, the increase in latent heat flux is actually about 15 % smaller for the micro- than for the macro-scale lateral fluxes, due to the different ways the water is redistributed within the grid-cell system (Fig. 5h and i). Here, macro-scale runoff increases the soil moisture in a comparatively large area, i.e. the lowland region, while the micro-scale transport processes predominantly lead to a moisture convergence in a comparatively small part of the grid cell, i.e. the depressions. And since evapotranspiration rates have a non-linear dependency on the water availability, most importantly they are limited by the available energy, a small increase in the soil moisture within a large region has a stronger effect on the latent heat flux than a large increase in soil moisture in a small region. The above, however, also means, that the effect on the maximum (spatial) temperature gradient is notably larger – on average by about 35 % – for the micro-scale fluxes (Fig. 5k and l). Due to the moisture convergence in a comparatively small part of the grid cell, the local temperatures within the depressions predominantly show a stronger reduction due to micro-scale fluxes, than the temperatures in the lowland fraction due to the macro-scale fluxes. Finally, it should be noted that the effects in the simulation that enable both micro- and macro-scale transport processes are not necessarily the same as combining the effects from the micro-scale-flux-enabling and the macro-scale-flux-enabling simulations. The lateral-movement of water in the fully-enabled simulation is very close to the sum of the effects in the two simulations that only enable one flux component (Fig. 5a–c). However, in case of the terrestrial water storage the former is about 25 % smaller than the sum of the effects (Fig. 5d–f) and for the effect on the latent heat flux it only amounts to about 60 % (Fig. 5g–i). This discrepancy arises because the terrestrial water storage and the latent heat flux are not only limited by the available water, with former additionally being limited by the pore volume and the potential depression storage and the latter by the available energy. Thus, an increase in the laterally transported water does not necessarily entail an increase in the water storage or in evapotranspiration.

3.2 Lateral heat transport

In the second application, we aim to investigate the effect of the lateral heat transport on the subgrid-scale temperature variability. Here, we set up the model to represent typical patterned-ground structures – e.g. non-sorted circles, ice-wedge polygons – that are often found in the arctic tundra. The presence of such patterns is limited to periglacial regions and, while the coverage of environmentally suitable spaces has been estimated at around 3 million km² for ice-wedge polygons (Karjalainen et al., 2020), other studies conclude the actual extent of the polygonal tundra to be merely around 0.3 million km² (Höfle et al., 2013). However, these structures present an ideal test case for our model, since soil organic matter concentrations, and with that hydrological and thermal soil properties, can vary dramatically on the sub-meter scale (Ping et al., 2015). For this application, we assume circular structures, that consists of an organic-rich center (C), with a volumetric organic matter fraction of 95 % at the surface – i.e. the uppermost 0.1 m of the soil, encompassed by the rim section (R), which has a near-surface soil organic matter fraction of 35 %. These two are in turn surrounded by an outer section (O) with mainly mineral soil throughout the below ground column, i.e. an organic matter fraction of 5 % at the surface, resulting in the following connectivity matrix:

$${\mathbf{\Lambda }}_{\mathrm{ctc}}=\begin{array}{cccc}& C& R& O\\ C& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ R& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ O& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}.$$

We assume the radius of these circles to be evenly subdivided into the center, rim and outer section, with water draining from the circle either at the soil bedrock interface or laterally through the outer section, with the former only being possible if the active layer is sufficiently deep. For simplicity reasons, we assume these circles to be omnipresent in the permafrost region, allowing us to conduct our analysis for the pan-Arctic average rather than focusing on a specific grid cell. We conduct ICON-Land-standalone simulations driven by climate forcing from the Global Soil Wetness Project Phase 3 (GSWP3, Dirmeyer et al., 2006; Kim, 2017). The simulations were performed for the period 1979–2018, with the analysis using the model output from the period 2009–2018. For different setups we compare the spatial, below-ground temperature variability, i.e. the temperature differences between center, rim and outer section, between simulations in which the lateral heat transport is disabled and simulations with the heat transport active.

In a first step, we assume a radius of one meter, with center, rim and outer section each covering 0.33 m. Additionally, we assume the circle to be level so that there is no vertical offset between the surface of the sections. The latter may not be the most common configuration, since the circular structures often have a distinct micro-topography, e.g. an elevated center in case of a hummock-type structure, but it makes the analysis much more straight-forward, since effects due to the lateral water transport and a lateral coupling with a vertical surface offset are (largely) absent. This provides the following connectivity matrices for connections due to preferential hydrological flow paths, assuming that they are the same for above- and below-ground fluxes (note that macro-scale fluxes are not present in this particular setup):

$${\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{micro}}=\begin{array}{cccc}& C& R& O\\ C& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ R& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ O& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}.$$

More importantly, this setup allows us to compare our simulations with those of the Dynamic Soil Model (Thurner et al., 2026, DynSoM), a pedon-scale soil model, which was developed to explore the movement of energy specifically in permafrost-affected soils. DynSoM simulates the key physical variables, such as soil temperature, and soil water and ice content, as well as the surface and below-ground heat fluxes both in the vertical and horizontal direction (note that the lateral movement of water is not represented in the present DynSom setup). Here, DynSoM, in contrast to ICON-Land, is a true 2D model which explicitly resolves the soil in one horizontal direction. For the present setup, the model uses a 10 cm horizontal grid-spacing and a variable spacing in the vertical direction, with the uppermost meter of the soil column being resolved in 10 cm steps. The simulations with DynSoM are based a 1 m soil transect of a non-sorted circle in the Chersky region (Gentsch et al., 2015; Supplements profile CH-E [note that in the original transect the inner circle – i.e. the circle center – features the low organic matter concentration and the inter circle area – i.e. the outer section – the high concentrations, but the order can simply be reversed without affecting the results of the DynSoM simulation]) and are forced with CRUNCEP data (Viovy, 2018) for the respective region. The setup of the DynSoM soil transect differs to the synthetic ICON-Land setup in that the organic rich region (surface organic matter fraction of 95 %) and the organic poor region (organic matter fraction of 5 %) both only cover 0.2 m, while the intermediary region (surface organic matter fraction of 35 %) covers 0.6 m. Furthermore, the mineral soil below the top layer also shows differences along the transect, while the present ICON-Land setup assumes spatially homogeneous mineral soil properties within one grid cell. However, since we do not attempt a site-level validation of our model but merely aim to compare its general behavior to a model that explicitly resolves the horizontal exchange on the pedon-scale, we assume that respective setups of the models are sufficiently similar.

Figure 6 Effects of the lateral heat transport on spatial temperature differences in patterned ground. Shown are the monthly temperature differences until a depth of 5 m between the organic rich circle center and the rim (first and third column) and between the center and the outer section (second and fourth column), for simulations that do not account for lateral transport processes (first and second column) and simulations where these processes are included in the model (third and fourth column). The first row shows the differences between simulations with the DynSoM model, for a flat-centered circle with a total radius of 1 m. The second row shows the corresponding ICON-Land simulations. The third row shows the temperature differences in a ICON-Land simulations for a circle radius of 10 m and the fourth row for a circle radius of 100 m. The fifth row shows the temperature differences for a ICON-Land setup of a low-centered circle with a 1 m radius and the last row for a setup of a high-centered circle with a 1 m radius.

In DynSoM, soil organic matter predominantly lowers the soil heat conductivity, which impedes the ground heat fluxes during the snow-free season (not shown). When neglecting the lateral exchange, this entails notably lower soil temperatures in the organic rich center than in the rim and the outer section, with differences at the surface starting to emerge after the snowmelt in early summer (Fig. 6a and b). The signal subsequently propagates downwards and the peak temperature differences in depth > 2 m are being reached the following spring. At the same time, the lower heat conductivity in the center also results in a more pronounced heating of the surface during late summer and a less pronounced cooling during early fall, with the near-surface temperatures exceeding those in the rim and the outer section between August and October. However, as soon as the soil is insulted by an increasingly thick snow cover in fall, also the near surface temperatures in the center fall below those of the rim and the outer section. During the cold season, the lower heat conductivity of the organic matter plays a subordinate role mainly because the soil is insulted by the snow cover. Additionally, ice has a much higher heat conductivity than liquid water and the relative differences in the near-surface heat conductivity between the sections are much smaller when the soils are frozen. Consequently, the cooling of the center, relative to the rim and the outer section, during summer is not fully balanced by a relative warming during the cold season. A temperature convergence during the cold season is observable in some soil layers, however, the temperatures in the center remain below those in the rim and the outer section even in these layers. Here, the temperature differences between center and rim and between center and outer section are largely similar. However, the signal is slightly more pronounced for the center-rim comparison, despite the near-surface soil organic matter concentrations in the center actually being closer to those in the rim than those in the outer section. The reason for this lays in the fact that, at depth > 5 m, the rim has a higher heat conductivity than the outer section, which has a notable effect on the magnitude of the bottom heat flux (not shown). Here, DySoM does not employ a zero-flux condition, but a prescribed temperature. With 0 °C, this temperature is higher than the simulated soil temperatures in the superjacent layers resulting in a positive bottom heat flux. Owing to the higher conductivities, the rim exhibits a larger warming at the bottom of the soil column, hence, higher temperatures than the outer section. Thus, the temperature differences to the cooler center are also more pronounced for the rim.

In the ICON-Land simulation, the temperature differences between the sections are equally driven by differences in the top- and f). Consequently, the overall effects are quite similar between the two models, especially given the fact that we are comparing simulations of a specific site (DynSoM) to the pan-Arctic average based on a synthetic setup (ICON-Land). However, there are also important differences, most notably in that large temperature differences between the sections are confined to the uppermost meter of the column and that the near-surface temperatures during the cold season are actually higher in the circle center. Both of these differences are related to the hydraulic properties of the soil organic matter – which can differ quite strongly to those of the mineral soil fractions – and to the assumed drainage pathways of the circle. Contrary to the DySoM simulation, the soil in the ICON-Land simulation is not necessarily water saturated and the degree of saturation of the soil in a given circle section depends on the overall water holding capacity, the infiltration rates and the position on the drainage path. Here, organic matter has a higher porosity than the mineral soil fractions – hence can hold more water – and also a higher hydraulic conductivity – hence allows for higher infiltration rates. Furthermore, the circle center drains laterally through the rim and the outer section and without vertical offset between the three sections, lateral fluxes mainly occur if the rim and the outer section are less saturated than the circle center. As a result, the soil column is more water saturated in the organic rich center than in the rim and the outer section, which entails higher heat conductivities below the organic layer (not shown). During the snow-free season, when the soil warms from the top, the combination of lower heat conductivities at the surface and higher heat conductivities below does not only reduce the overall heat uptake, but also enhances the downward transport, resulting in a reduced vertical temperature gradient. This difference in heat distribution causes the cooling of the center (relative to the rim and the outer section) to be more pronounced closer to the surface and less pronounced deeper within the soil. Similarly, the higher heat conductivities also lead to a reduced vertical temperature gradient during the cold season. However, during this period, the column cools from the top and the reduced temperature gradient corresponds to an enhanced upwards heat transport in the center and, consequently, higher near-surface temperatures than in the rim and the outer section. Here, the respective effects are more pronounced in the center-outer-section than in the center-rim comparison, i.e. larger positive and negative temperature differences close to the surface and a slightly less pronounced relative cooling below depths > 1 m. This is because the near-surface soil organic matter concentrations in the center are closer to those in the rim than those in the outer section and because the soil column in the latter, as the outermost circle section, is predominantly less water saturated than in the rim.

As discussed above, there are some differences in the general effect of a heterogeneous soil organic matter distribution between the models. Nonetheless, they agree well on the impact of the lateral heat fluxes on the spatial temperature variability. In both simulations, the temperature differences between the circle sections are strongly reduced and any notable spatial variability is largely limited to the uppermost meter of the column as well as to the snow-free season (Fig. 6c, d, g and h). Here, the fact that the lateral heat fluxes have a stronger influence on the temperature differences at greater depths has two main causes. On the one hand, the temperature differences deeper within the soil are not only balanced by the lateral heat fluxes on the respective layers, but they are also reduced since the vertical heat exchange with the superjacent layers is spatially more homogeneous due to the lateral heat exchange in the layers above. On the other hand, the temporal temperature variability is strongly dampened with increasing depth. Thus, while large near-surface temperature differences between the sections may not be fully balanced on the short term, the signal is attenuated while propagated downwards to the extent that temperature differences in the deeper soil layers only emerge if the longer-term mean temperature near the surface is notably different.

Overall, the above results suggest that, at the pedon-scale, the lateral heat fluxes are sufficiently strong to largely balance the below-ground temperature differences resulting from the spatial variability in the surface heat fluxes. Furthermore, it appears that these effects can be captured reasonably well by LSMs employing a tiling-approach. The latter allows us to use ICON-Land simulations to investigate the extent to which the ability of the lateral heat exchange to the equilibrate soil temperatures depends on the horizontal distances considered. For a radius of 10 m, with the circle sections each covering 3.33 m, we find that the temperature differences between center, rim and outer section above 1 m, appear to be hardly affected by the lateral heat exchange (Fig.6i–l). This strongly suggest that at these length scales the lateral heat exchange is too weak to balance the short-term effects resulting from the differences in the seasonal variability of the surface heat fluxes. However, the near-surface temperatures in the center do in fact show a slight longer-term increase due to the lateral fluxes (e.g. for the center-outer-section comparison this is visible in the [positive] temperature differences during the cold season being slightly larger) and the resulting annual mean temperature is very close to that in the rim and the outer section. Consequently, in those soil layers that are less sensitive to the seasonal variability in the surface heat flux, the spatial temperature variability is reduced substantially. At depths > 1 m, there are no marked temperature differences neither between the center and rim nor between center and the outer section. Here, it should be noted that the absence of large temperature differences comparatively close to the surface is partly caused by the rather coarse vertical resolution of the model, since the damping of the temporal temperature variability depends on the latter. With a thickness of almost 3 m, the first soil layer below a depth of 1 m possesses a substantial thermal inertia and the temporal variability of the vertical heat fluxes only has a small effect on the temperatures of this layer. Thus it is highly plausible that, in reality or with a higher vertical resolution, the spatial temperature variability extends much further downward than in our simulation. Increasing the radius to 100 m, with the sections each covering 33.33 m, further reduces the lateral heat fluxes. As a result, the latter are not only too weak to balance the short-term effects due to the variability of the surface heat fluxes, but also the longer-term mean state of the near surface layers is substantially different. Consequently, notable differences in the below-ground temperatures persist even at depths > 4 m and, for these length scales, the impact of the lateral heat fluxes appears to become negligible (Fig. 6m–p).

As stated above, a level surface of the circle may not be the most common configuration as in reality we often find high-centered hummock-type structures, or low-centered hollow-type structures. On the one hand, the vertical offset between the sections of the circle affects the soil thermophysical properties, since it determines the flow of water, hence, the saturation of the soil and, thusly, the heat conductivity. On the other hand, the ground features a strong vertical temperature gradient during most of the year, with the vertical offset between two sections determining the distance to the surface of any two interfacing points (Fig. 3). Thus, in the final example we investigate the extend to which the micro-scale topography of the circle effects the spatial temperature variability. First we focus on a low-centered circle with a 1 m radius, in which the surface of the center and outer section are both located 0.1 m below the surface of the rim, with the connectivities for hydrological flow paths being set in the following way:

$${\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{micro}}=\begin{array}{cccc}& C& R& O\\ C& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ R& \mathrm{0.5}& \mathrm{0}& \mathrm{0.5}\\ O& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}.$$

For, this configuration we find notable differences between the temperatures of the sections despite the small horizontal extent of the circle (Fig. 6q–t). The near surface soil layers in the center are cooler during summer and warmer during winter, with the differences being larger in the center-rim than in the center-outer-section comparison. Here, the temperature differences between center and rim are less determined by the differences in heat conductivities and mainly reflect the vertical temperature gradient across the uppermost 0.1 m of the soil. The lateral heat fluxes are driven by the temperature differences of two points with the same absolute vertical position, hence the surface temperatures of the center are essentially adjusting to the temperatures at 0.1 m depth in the rim. During the summer the (downwards negative) vertical temperature gradient near the surface reaches a maximum of about 20 K m⁻¹ and during winter roughly half that value (downwards positive), meaning that temperature differences of up to +2 and −1 K, are fully explainable by the vertical offset between the two tiles. Furthermore, since the vertical temperature gradient averaged over longer-periods is close to zero, there are no notable differences in the temperatures in the deeper soil layers. The temperature differences between center and the outer section (which have the same vertical offset to the rim) show, however, that there is also a notable variability in the temperatures of points with the same absolute vertical position, i.e. that the surface temperatures in the rim and the outer section do not fully adjust to the temperatures at 0.1 m depth in the rim. This variability is much larger for the low-centered than for the level circle (Fig. 6h and t), indicating that the lateral heat fluxes are smaller in case of the former. The reason for the reduced lateral exchange lies in the different levels of saturation of the rim. Due to its raised position – in case of the low-centered circle – the rim drains more readily into the other sections, resulting in overall drying of the soil column. This in turn leads to a lower heat conductivity and, consequently, reduced lateral fluxes. For a high-centered structure, in which the surface of the circle center is located 0.1 m above the surface of the rim, which in turn is located 0.1 m above the surface of the outer section – resulting in the following connectivity:

$${\mathbf{\Lambda }}_{\mathrm{hfp},\mathrm{micro}}=\begin{array}{cccc}& C& R& O\\ C& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ R& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ O& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}$$

we find that the lateral heat fluxes have essentially the same effect as for a low-centered circle. However, since the center is elevated relative to the rim and the outer section, it features higher near surface temperatures during summer and colder temperatures during winter (Fig. 6u–x), with the magnitude of the effects increasing with the vertical offset. As a result, the temperature differences, resulting from the lateral heat fluxes actually oppose the ones in the simulation in which the lateral fluxes are being neglected.

4 Summary and discussion

Rapid increases in feasible model resolutions raise the question of whether further development of statistical representations of land-surface heterogeneity and the associated lateral fluxes is still warranted. But while much of what we consider macro-scale spatial heterogeneity is already, or can soon be, explicitly resolved even in multi-decadal simulations, many micro-scale structures – such as the surface depressions or patterned ground considered in our example applications – have characteristic length scales that are orders of magnitude smaller than feasible grid resolutions. Moreover, the ability to resolve landscapes at increasingly fine spatial scales does not necessarily imply that interactions between connected landscape elements are well represented in the models. In fact, the standard configurations of most current-generation land surface models lack an explicit representation of lateral exchanges between grid cells (beyond river routing), and only a few research extensions and coupled frameworks implement such lateral interactions (Qiu et al., 2024). Importantly, increases in horizontal resolution also raise computational costs exponentially: for example, refining resolution from a 100 to 1 km grid spacing increases the number of grid cells by a factor of 10 000, whereas the computational demand associated with a more complex statistical representation of surface heterogeneity scales approximately linearly with the number of additional tiles. Consequently, increasing spatial resolutions and more sophisticated tiling approaches may each be more suitable to address different research questions and will likely remain complementary strategies for representing land-surface heterogeneity in the foreseeable future.

In the present work, we described a new scheme that provides a lateral coupling between the tiles capturing the spatial subgrid-scale heterogeneity of the land surface in ICON-Land, the land component of the ICON framework. The scheme represents 5 lateral exchange processes, namely gravity-driven moisture fluxes and the corresponding advective heat transport, diffusive- and conductive fluxes of water and heat as well as a parameterized snow-redistribution component that approximates the net effect of wind-driven snow transport. The main assumption the scheme is based on is that the relationships between any two homogeneous surface clusters, the tiles, can be described by a set of connectivities. The latter result from the internal logic underlying the tile-definition and reflect the relative (geometrical) contact lengths between two tiles and the dominant hydrological flow paths connecting the two. Taking into account the composition of a given grid cell, they are used to calculate the spatial relationships between two tiles and the time-lag factors that govern the lateral transport. Here, the actual parametrizations of the scheme were chosen to integrate well with the model's existing description of the vertical transport processes, e.g. the scheme merely transports the runoff that is calculated within ICON-Land's point-scale hydrology module, while the time-lag factors for the macro-scale fluxes of water are based on the assumptions underlying the ARNO-rainfall–runoff model, which is used to describe the runoff from catchment-sized areas in standard ICON-Land-setups.

In it is present implementation, the scheme does not provide a final suit of parametrizations that necessarily work for all types of land-surface heterogeneity. And, given the vast range of land processes and their dependencies on different soil and surface properties – all with potentially different spatial distributions, it is highly questionable whether a setup accounting for the entirety of land-surface heterogeneity can even be run with reasonable computational costs. Instead it appears more likely that the model setup, hence, the required lateral-flux parametrizations, will depend on the scientific question that a specific simulation is used to help answer. Thus, our goal was to provide a framework that allows for a consistent treatment of the inter-tile exchange, with enough flexibility to facilitate the implementation of improved or adjusted parametrizations and new processes.

The showcases discussed above represent highly idealized applications that focus on a small number of surface and soil characteristics. The scheme is, however, not necessarily limited to such configurations, and more realistic cases can be constructed utilizing its infrastructure. One way to extend the approach to increasingly complex landscape representations is to use a more extensive tiling configuration that captures a broader range of surface heterogeneity while still providing consistent descriptions of the relationships between tiles. For example, the topography-based tiling used here could readily be expanded to include additional (quasi-) topographic landscape elements, such as floodplain tiles that predominantly interact with lowland areas. Tiles could also be further subdivided according to additional properties such as soil type or vegetation. However, defining globally applicable connectivity rules describing tile interactions becomes increasingly difficult as the number of tile characteristics increases. In the absence of such rules, one might assume a random distribution for certain characteristics within a grid cell, implying that the respective connectivities could be treated as being proportional to the cover fractions of the tiles. In many cases, however, this would be a poor approximation. Achieving higher degrees of realism with such a strategy would therefore require additional flexibility in the treatment of tile connectivities, and incorporating grid-cell-specific connectivity matrices derived from external data describing landscape structure represents a natural direction for future model development.

Another strategy, which is already feasible with the current model capabilities, is to define tiles primarily by their characteristic connectivities rather than by specific surface or soil properties. In this case, the connectivities between tiles are globally uniform, while the land-surface characteristics associated with each tile may vary between grid cells. All tile properties, including vegetation, soil parameters, and cover fractions, can be provided to the model individually for each grid cell. High-resolution data could therefore be used to identify areas within a grid cell that, due to their relative position and topographic setting, are likely to interact with their surroundings in a way that matches the connectivity pattern prescribed for a given tile type. For example, a tile A may always interact with tiles B and C in the same prescribed manner, while the surface properties associated with these tiles differ between grid cells. In one grid cell, tile A might represent shrub-covered areas on loamy soils covering only a small fraction of the grid cell, while tiles B and C represent forest and grassland areas on predominantly sandy soils. In another grid cell, tile A could instead occupy a much larger fraction of the landscape and consist mainly of bare ground, while tiles B and C represent wetlands and shrublands. In this way, the connectivity structure remains fixed, while the land-surface characteristics associated with each tile can vary across the domain.

There are, however, certain aspects which limit the degree of realism that can be achieved with the present implementation of the scheme. One of these is the assumed circular shape of the subgrid scale patches. The representation of the patch geometry in the spatial relationship module is necessarily a simplification, which has important implications for the simulated lateral exchanges. The lateral fluxes depend on the characteristic dimensions of the adjacent patches, and any assumption about their shapes will influence the magnitude of the fluxes: The diffusive fluxes scale linearly with the contact length between two connected patches and are inversely proportional to the center-to-center distance, while subsurface runoff is proportional to the slope, which depends on the height difference between adjacent patches relative to the center-to-center distance. In the present implementation, the circular geometry was chosen as a pragmatic compromise, as it provides a reasonably general approximation for many landscape features, including small lakes, ponds, and patterned terrain, while allowing the patch dimensions to be described by a single parameter, the characteristic area, which can easily be determined from high-resolution observational data. For highly elongated or irregular features, however, these assumptions may lead to over- or underestimation of the lateral fluxes by up to an order of magnitude, potentially affecting the ability of the fluxes to equilibrate the states of neighboring tiles (see Sect. 3.2). To overcome these limitations, future developments could improve the flexibility of the scheme by introducing (tile- and grid-cell specific) shape factors (S₁, S₂) in the key ratios of perimeter (C) and radius (r):

$$\begin{array}{}\text{(31)}& {\displaystyle }& {\displaystyle }C=\mathrm{2}\mathit{\pi }r{S}_{\mathrm{1}}\text{(32)}& {\displaystyle }& {\displaystyle }r=\sqrt{A/\mathit{\pi }}{S}_{\mathrm{2}},\end{array}$$

where A is the area of the tile.

Another problematic aspect is the highly simplified approach to incorporate snow redistribution. In natural landscapes the latter is a complex process, driven by a combination of topography, snowpack properties, vegetation, and especially near-surface wind direction and speed (Clark et al., 2011). Land surface models typically do not represent this process explicitly but some do account for topographic variability when relating (grid-cell mean) snow water equivalent to the snow cover fraction, though not in a way that takes into consideration wind speed or direction (Roesch and Roeckner, 2006; Ma et al., 2019). Similarly, simple topography-based schemes have been employed successfully in previous studies, to account for the effect of snow redistribution between microtopographic features in permafrost landscapes (Aas et al., 2019; Nitzbon et al., 2019; Smith et al., 2022). However, while these approaches may work sufficiently well on the micro-scale, it can be argued that for the hillslope scale the wind direction may no longer be ignored where drifting leads to preferential deposition in the lee of elevations and on sheltered aspects, making them less suitable for the macro-scale transport.

While the current snow redistribution scheme is based solely on differences in surface elevation and neglects near-surface wind effects, it is conceptually possible to extend the approach to account for wind direction and magnitude. In principle, direction-specific connectivity matrices could be pre-computed during the model setup, reflecting the preferential pathways for snow transport under prevailing wind directions. Furthermore, the connectivity between two tiles could be scaled according to the fraction of the boundary aligned with the wind in combination wind-speed-dependent redistribution functions which could, for example, be derived by running dedicated snow pack and snow drift models (Quéno et al., 2024). However, while a wind-aware formulation is feasible in principle, it is beyond the scope of the current study and is left as a potential avenue for future model development. It is important to emphasize that while the representation of snow redistribution is a key process for capturing subgrid-scale hydrological and energy dynamics in certain regions, its omission or inclusion has only minor impacts on the overall patterns and magnitudes of the integrated fluxes and state variables presented in Fig. 5. This is because, regardless of whether snow redistribution is active or not, the majority of snowmelt water ultimately converges toward small-scale depressions – either through drifting snow or surface and subsurface runoff. Consequently, enabling or disabling snow redistribution primarily affects the timing of lateral fluxes, rather than their overall magnitude or system-scale impact. A key exception occurs during the snowmelt season, when the subgrid-scale distribution of the snow cover influences the spatial pattern of surface energy fluxes and temperatures, as energy is required to melt snow in areas where it accumulates – either on small-scale elevations or in depressions. However, this spatial variability does not change the grid-cell mean fluxes shown in Fig. 5 substantially. Importantly, is also does not significantly alter the maximum subgrid-scale temperature variability, since this is primarily governed by water availability during the summer months.

Despite these shortcomings of the model and its incompleteness, we hope that the two example applications demonstrate the potential usefulness of our scheme. Here, our simulations showed that lateral-transport processes are a key factor determining the subgrid-scale temperature and moisture variability. Among other things they showed that the model's ability to simulate standing water at the surface, was mainly determined by the inclusion of the lateral water movement, with the subrid-scale temperature variability increasing by an order of magnitude across most of North America and Eastern Siberia. Furthermore, the small-scale (spatial) variability in soil temperatures was highly sensitive to the lateral heat fluxes, with the magnitude and even the direction of effects depending on the assumed spatial-scales and the topography.

It should be noted, however, that in our simulations the grid-cell mean state and exchange with the atmosphere did not show substantial changes due to the inclusion of the lateral transport processes. Even more fundamentally, there was very little to suggest that capturing the general physical land-atmosphere interactions necessarily requires a tiling of the land surface and even the wetland extent can most likely be represented sufficiently well with an appropriate parametrization of the depression storage and the retention due to the flow along slopes. In large parts this may be attributable to the nature of the experiments we conducted, since we ran the simulations with prescribed atmospheric conditions, excluding all land-atmosphere feedback effects. Here, fully coupled simulations with a similar setup showed that even comparatively small changes in Bowen ratio can have substantial climate effects in regions where the state of the atmosphere, in particular the low-altitude cloud cover, is sensitive to changes in terrestrial evapotranspiration (de Vrese et al., 2024). Thus, while our offline results show negligible changes in grid-cell mean fluxes and states, coupled simulations would be needed to assess any feedback on climate fully. Furthermore, a tiling may still be required to represent those processes that have a highly non-linear dependency on the state of the (sub) surface and which may not be well be described based on the grid-cell mean state. Here, one of the most prominent example are the terrestrial methane fluxes, which are dominated by the emissions from soils and wetlands. The latter are the result of anaerobic decomposition processes and, therefore, require anoxic conditions. Since fully water saturated soils can mostly only be found in a fraction of any grid-cell-sized area, a tiling approach may still be the most valid strategy to represent the state and fluxes in the respective areas as it allows for a consistent treatment of the local hydrological, thermophysical and biochemical processes. And if tiles are being used to represent different surface clusters, the respective lateral coupling may not be ignored. In case of methane emissions, for example, the lateral exchange fluxes appear to not only have a substantial impact on the extent of areas with fully water saturated soils but also on the subgrid-scale temperature distribution, another key determinant of the local decomposition rates.

Appendix A: Overview of variables

Table A1 Comprehensive overview of T-REX variables in alphabetical order.

Appendix B: Soil hydrology in ICON-Land

ICON-Land includes two approaches for determining infiltration, surface- and subsurface runoff, the first of which is suitable for coarse-resolution simulations and is based on the ARNO-rainfall–runoff model (Dümenil and Todini, 1992; Todini, 1996; Reick et al., 2021). The ARNO model determines the partitioning of the moisture fluxes at the surface based on an assumed subgrid-scale soil-moisture distribution. Conceptually, the soil moisture in any grid cell is subdivided into a number of local water storages, the water content (V_w) of which is assumed to follow the cumulated distribution function:

$$\begin{array}{}\text{(B1)}& f_{\left(\mathrm{w}\right)}=\mathrm{1}\mathrm{-}{\left(\mathrm{1}\mathrm{-}{\displaystyle \frac{V_{\mathrm{w}}}{V_{\mathrm{w},\max }}}\right)}^b,\end{array}$$

where V_w,max is the maximum water-holding capacity of the local storages and b a steepness-parameter, which accounts for the subgrid-scale topography and is estimated based on the standard deviation of the topographic height (σ_z) and two fitting parameters (σ₀ and σ_mx), whose standard values are 100 and $\mathrm{1000}\frac{\mathrm{64}}{n_{\mathrm{lat}}}$ (where n_lat is the number of grid cells in latitudinal direction considered in a setup) respectively.

$$\begin{array}{}\text{(B2)}& b=\left\{\begin{array}{ll}\mathrm{0.5}& \mathrm{if}{\mathit{\sigma }}_{mx}<{\mathit{\sigma }}_z\\ \frac{{\mathit{\sigma }}_z-{\mathit{\sigma }}_{\mathrm{0}}}{{\mathit{\sigma }}_z-{\mathit{\sigma }}_{mx}}& \mathrm{if}{\mathit{\sigma }}_{\mathrm{0}}<{\mathit{\sigma }}_z<{\mathit{\sigma }}_{mx}\\ \mathrm{0.01}& \mathrm{otherwise}\end{array}\right..\end{array}$$

Further assuming that the surface fraction on which runoff occurs $\left(\frac{A_{\mathrm{Rsrf}}}{A_{\mathrm{cell}}}\right)$ for a given grid-cell-mean water content is described by the above distribution function –that is $\frac{A_{\mathrm{Rsrf}}}{A_{\mathrm{cell}}}=f_{\left(\mathrm{w}\right)}$, infiltration (q_I) can be determined as the residual of precipitation (q_P) – more specifically that fraction of precipitation that is not intercepted by vegetation, but including snow melt –and surface runoff (q_R,srf):

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& q_{\mathrm{I}}=q_{\mathrm{P}}-q_{\mathrm{R},\mathrm{srf}}\mathrm{with}\\ & q_{\mathrm{R},\mathrm{srf}}=q_{\mathrm{P}}-{\displaystyle \frac{\left(V_{\mathrm{rz},\max }-V_{\mathrm{rz}}\right)}{\mathrm{\Delta }t}}+{\displaystyle \frac{\mathrm{1}}{(b+\mathrm{1}{)}^{b+\mathrm{1}}\cdot V_{\mathrm{rz},\max }^{\mathrm{b}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\\ & \left\{\begin{array}{ll}{\left[(b+\mathrm{1})\cdot V_{\mathrm{rz},\max }{\left(\mathrm{1}-\frac{V_{\mathrm{rz}}}{V_{\mathrm{rz},\max }}\right)}^{\frac{\mathrm{1}}{\mathrm{1}+b}}-q_{\mathrm{P}}\cdot \mathrm{\Delta }t\right]}^{b+\mathrm{1}}& \mathrm{if}(b+\mathrm{1})\cdot V_{\mathrm{rz},\max }{\left(\mathrm{1}-\frac{V_{\mathrm{rz}}}{V_{\mathrm{rz},\max }}\right)}^{\frac{\mathrm{1}}{\mathrm{1}+b}}>q_{\mathrm{P}}\cdot \mathrm{\Delta }t\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right.\end{array}\end{array}$$

with V_rz being the water content of the root zone and V_rz,max the maximum root-zone soil moisture.

Lateral drainage or subsurface runoff q_R,blg is calculated using two fixed drainage rates that are scaled by the saturation of the soil:

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}& q_{\mathrm{R},\mathrm{blg}}=\\ & \left\{\begin{array}{ll}{d}_{mn}\cdot \frac{V_{\mathrm{rz}}}{V_{\mathrm{rz},\max }}+\left(d_{mx}-d_{mn}\right)\cdot {\left[\frac{V_{\mathrm{rz}}-V\mathrm{rz},\mathrm{crt}}{V_{\mathrm{rz},\max }-V_{\mathrm{rz},\mathrm{crt}}}\right]}^{d_{xp}}& \mathrm{if}{V}_{\mathrm{rz}}>V_{\mathrm{rz},\mathrm{crt}}\\ d_{mn}\cdot \frac{V_{\mathrm{rz}}}{V_{\mathrm{rz},\max }}& \text{otherwise}\end{array}\right.\end{array}\end{array}$$

with

$$V_{\mathrm{rz},\mathrm{crt}}=V_{\mathrm{rz},\max }\cdot f_{\mathrm{dr},\mathrm{crt}}$$

where d_mx ($\mathrm{2.81}\times {\mathrm{10}}^{-\mathrm{8}}$) corresponds to the maximum drainage that occurs under conditions close to saturation – that is for the soil-moisture values about a critical threshold ($V_{\mathrm{rz},\max }\cdot f_{\mathrm{dr},\mathrm{crt}}$, where $f_{\mathrm{dr},\mathrm{crt}}=\mathrm{0.9}$) –and d_mn ($\mathrm{2.81}\times {\mathrm{10}}^{-\mathrm{10}}$) to the maximum drainage rate in unsaturated soils. d_xp (1.5) is a drainage parameter used for model tuning. Assuming fully water-saturated soils – i.e. $q_{\mathrm{R},\mathrm{blg}}=d_{mx}$ – the above parameter values allow for a volume of water corresponding to a typical excess-water volume (i.e. $V_{\mathrm{rz},\max }-V_{\mathrm{rz},\mathrm{crt}}$ between 0.005 and 0.03 m (Todini, 1996; Dümenil and Todini, 1992)) to drain from a catchment-sized area within 2.1–12.5 d.

In the second, the point-scale, approach soil moisture is assumed to be distributed evenly across the domain and the formation of surface runoff is exclusively based on the speed with which water can infiltrate into and move vertically through the soil. Here, infiltration (q_I) and surface runoff (q_R,srf) are determined in two steps, calculating the contribution of infiltration-excess (Horton overland flow; $q_{\mathrm{R},\mathrm{srf},\mathrm{h}}$) and saturation-excess (Dunne overland flow; $q_{\mathrm{R},\mathrm{srf},\mathrm{d}}$) separately. In the first step, the infiltration capacity (q_I,cap) is calculated as the saturated hydraulic conductivity of the uppermost soil layer, k_sat, scaled by an ice-impedance factor (Swenson et al., 2012), f_imp:

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& q_{\mathrm{I},\mathrm{cap}}=k_{\mathrm{sat}}{f}_{\mathrm{imp}}\mathrm{with}\\ & f_{\mathrm{imp}}={\mathrm{10}}^{-\mathrm{6}\frac{{\mathit{\theta }}_{\mathrm{ice},\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{sat},\mathrm{1}}}},\end{array}\end{array}$$

where θ_ice,1 is the ice content of the uppermost soil layer and θ_sat,1 the layer's porosity.

Horton overland flow, $q_{\mathrm{R},\mathrm{srf},\mathrm{h}}$, occurs if the available water exceeds the infiltration capacity and the surface-storage capacity, q_U. Here, the model distinguishes between the water that can potentially infiltrate across the entire surface area ($q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}$) and the water, $q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}$, that may additionally infiltrate within the inundated fraction c_pnd. Consequently, the Horton overland flow consists of two elements, the infiltration-excess that is generated across the entirety of the surface area, $q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}^{\mathrm{t}}$, and the additional overflow of the surface depression storage $q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}^{\mathrm{t}}$:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}& q_{\mathrm{R},\mathrm{srf},\mathrm{h}}=q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}+q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}\mathrm{with}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{sl}}=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{sl}}-q_{\mathrm{I},\mathrm{cap}}& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot},\mathrm{sl}}>q_{\mathrm{I},\mathrm{cap}}\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{pnd}}=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}-q_{\mathrm{I},\mathrm{cap}}\cdot c_{\mathrm{pnd}}-q_{\mathrm{U}},& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}>q_{\mathrm{I},\mathrm{cap}}\cdot c_{\mathrm{pnd}}+q_{\mathrm{U}}\\ \mathrm{0},& \mathrm{otherwise}.\end{array}\right.\end{array}\end{array}$$

Potential infiltration rates depend on precipitation q_P, or throughfall in vegetated areas, and on the water that is stored within surface depressions. Here, the assumption is being made that the precipitated water does not infiltrate uniformly across the surface area but predominantly pools within the depressions:

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}& q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}=q_{\mathrm{P}}\cdot \left(\mathrm{1}-c_{\mathrm{pnd},\max }^{\frac{\mathrm{1}}{\mathrm{3}}}\right)\mathrm{and}\\ & q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}=q_{\mathrm{P}}\cdot c_{\mathrm{pnd},\max }^{\frac{\mathrm{1}}{\mathrm{3}}}+{\displaystyle \frac{V_{\mathrm{sfc}}}{\mathrm{\Delta }t}},\end{array}\end{array}$$

where c_pnd,max constitutes the maximum inundated fraction and V_sfc is the current water content of the surface reservoir. The surface storage capacity is given by the maximum water content of the surface reservoir V_sfc,max:

$$\begin{array}{}\text{(B8)}& q_{\mathrm{U}}={\displaystyle \frac{V_{\mathrm{sfc},\max }}{\mathrm{\Delta }t}}.\end{array}$$

In the second step an upper limit for infiltration, q_I,max, is determined based on $q_{\mathrm{I},\mathrm{pot},\mathrm{sl}}$, $q_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}$, q_I,cap and the inundated fraction c_pnd:

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}& q_{\mathrm{I},\max }=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{sl}}+q_{\mathrm{I},\max ,\mathrm{pnd}}& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot},\mathrm{sl}}+q_{\mathrm{I},\max ,\mathrm{pnd}}<q_{\mathrm{I},\mathrm{cap}}\\ q_{\mathrm{I},\mathrm{cap}}& \mathrm{otherwise}\end{array}\right.\mathrm{with}\\ & q_{\mathrm{I},\max ,\mathrm{pnd}}=\left\{\begin{array}{ll}{q}_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}& \mathrm{if}{q}_{\mathrm{I},\mathrm{pot},\mathrm{pnd}}<q_{\mathrm{I},\mathrm{cap}}\cdot c_{\mathrm{pnd}}\\ q_{\mathrm{I},\mathrm{cap}}\cdot c_{\mathrm{pnd}}& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

This rate is used as the upper boundary condition in the computation of the vertical soil hydrology. After the state of the soil column and the drainage fluxes have been updated any water that exceeds the pore volume is removed from the soil forming the saturation-excess $q_{\mathrm{R},\mathrm{srf},\mathrm{d}}$. The latter is then used to estimate the actual infiltration and the total surface runoff:

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}& q_{\mathrm{I}}=q\mathrm{I},\max -q_{\mathrm{R},\mathrm{srf},\mathrm{d}}\mathrm{and}\\ & q_{\mathrm{R},\mathrm{srf}}=q_{\mathrm{R},\mathrm{srf},\mathrm{h}}+q_{\mathrm{R},\mathrm{srf},\mathrm{d}}\mathrm{with}\\ & q_{\mathrm{R},\mathrm{srf},\mathrm{d}}=\sum _{i=\mathrm{1}}^{n_{\mathrm{soil}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{upd},i}-{\mathit{\theta }}_{\mathrm{sat},i}}{\mathrm{\Delta }t}}\mathrm{for}\mathrm{all}{\mathit{\theta }}_{\mathrm{upd},i}>{\mathit{\theta }}_{\mathrm{sat},i}.\end{array}\end{array}$$

With respect to the subsurface lateral drainage, the ARNO-scheme determines the fluxes implicitly assuming a given retention time, to account for the temporal lag between runoff generation and the water reaching the closest tributary. For the point-scale approach the generation of subsurface runoff was decoupled from these implicit assumptions and q_R,lat is instead determined as a function of the hydraulic conductivity, k_i, on a given layer i and the local slope m:

$$\begin{array}{}\text{(B11)}& q_{\mathrm{R},\mathrm{blg},\mathrm{i}}=k_i\cdot \sin \left(m\cdot {\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right).\end{array}$$

In addition to this lateral-drainage component, we included the (vertical) drainage from the lowest soil layer b into the underlying bedrock, q_B,b, with the drainage rates being limited by the hydraulic conductivity of the lowest layer, k_b, and the hydraulic conductivity of (fractured) bedrock, k_rock, assumed to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹:

$$\begin{array}{}\text{(B12)}& q_{\mathrm{B},\mathrm{b}}=\left\{\begin{array}{ll}{k}_{\mathrm{b}}& \mathrm{if}{k}_{\mathrm{b}}<k_{\mathrm{rock}}\\ k_{\mathrm{rock}}& \mathrm{otherwise}\end{array}\right..\end{array}$$

In the vast majority of cases, the model is being used to represent areas that have a given spatial extent. This poses the problem, that any runoff determined by the above point-scale parametrizations cannot be assumed to contribute to the stream flow instantaneously, since the water requires some time to move across a given distance either at the surface or through the ground before reaching a tributary. Here, the model makes use of intermediary reservoirs, V_L,srf and V_L,blg – analogues to those used in the retention parametrization in the hydrological discharge model – in which q_R,srf and q_R,blg are initially stored. The outflow of these reservoirs, $q_{\mathrm{R},\mathrm{srf}}^{*}$ and $q_{\mathrm{R},\mathrm{blg}}^{*}$, is estimated based on assumed retention times τ_srf and τ_blg and constitutes the contribution to the streamflow.

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}& q_{\mathrm{R},\mathrm{srf}}^{*}={\displaystyle \frac{V_{\mathrm{L},\mathrm{srf}}}{{\mathit{\tau }}_{\mathrm{srf}}}}\mathrm{and}\\ & q_{\mathrm{R},\mathrm{blg},i}^{*}={\displaystyle \frac{V_{\mathrm{L},\mathrm{blg},i}}{{\mathit{\tau }}_{\mathrm{blg}}}}.\end{array}\end{array}$$

For the surface runoff, infiltration or evaporation from the emerging flow channels are being neglected and the change in water content of the respective reservoirs, ΔV_L,srf, is simply given by:

$$\begin{array}{}\text{(B14)}& {\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{L},\mathrm{srf}}}{\mathrm{\Delta }t}}=q_{\mathrm{R},\mathrm{srf}}-q_{\mathrm{R},\mathrm{srf}}^{*}.\end{array}$$

The change in the water content on a given soil layer i of the reservoir retaining the lateral subsurface fluxes, $V_{\mathrm{L},\mathrm{blg},i}$, is calculated analogously, with the difference being that the interactions between the intermediary reservoir and the surrounding ground need to be accounted for. Here, it is assumed that the lateral drainage ceases if the ground begins to dry out and all excess water has been removed from the soil. This behaviour is represented by a restoration flux q_X,i, which migrates water back into the soil matrix as soon as the soil water content θ_i decreases below the soil field capacity θ_fc,i. Furthermore, the water within the intermediary reservoir may percolate downwards until it drains across the soil-bedrock interface, constituting additional bottom drainage q_B,i. We assume this flux to be limited by the hydraulic conductivity k_rock of (fractured) bedrock ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹) and the hydraulic conductivity k_b of the lowest soil layer above the bedrock boundary b. Water is taken from a given layer of the intermediary reservoir until the volume of water that has been removed is equal to the bottom-drainage flux multiplied by the time-step length or until the intermediary reservoir is completely empty. This iterative process is started at the top of the soil column, and not at the soil-bedrock interface to account for the vertical movement within the reservoir. Here, it is assumed that the water percolates downward with the same rate as it drains into the bedrock:

$$\begin{array}{}\text{(B15)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{L},\mathrm{blg},i}}{\mathrm{\Delta }t}}=q_{\mathrm{R},\mathrm{blg},i}-q_{\mathrm{R},\mathrm{blg},i}^{*}-q_{X,i}-q_{\mathrm{B},i}\mathrm{with}\\ & q_{X,i}=\left\{\begin{array}{ll}\mathrm{0},& \mathrm{if}{\mathit{\theta }}_{\mathrm{fc},i}>{\mathit{\theta }}_i\\ \frac{{\mathit{\theta }}_{\mathrm{fc},i}-{\mathit{\theta }}_i}{\mathrm{\Delta }t},& \mathrm{if}{\mathit{\theta }}_{\mathrm{fc},i}<{\mathit{\theta }}_i\mathrm{and}{\mathit{\theta }}_{\mathrm{fc},i}-{\mathit{\theta }}_i<V_{\mathrm{L},\mathrm{blg},i}\\ \frac{V_{\mathrm{L},\mathrm{blg},i}}{\mathrm{\Delta }t}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & q_{\mathrm{B},i}=\left\{\begin{array}{ll}\mathrm{0}.,& \mathrm{if}\frac{\sum _{l=\mathrm{1}}^{i-\mathrm{1}}{V}_{\mathrm{L},\mathrm{blg},l}}{\mathrm{\Delta }t}>k_{\mathrm{bot}}\\ \frac{V_{\mathrm{L},\mathrm{blg},i}}{\mathrm{\Delta }t}& \mathrm{if}\frac{\sum _{l=\mathrm{1}}^i{V}_{\mathrm{L},\mathrm{blg},l}}{\mathrm{\Delta }t}<k_{\mathrm{bot}}\\ k_{\mathrm{bot}}-\frac{\sum _{l=\mathrm{1}}^{i-\mathrm{1}}{V}_{\mathrm{L},\mathrm{blg},l}}{\mathrm{\Delta }t}& \mathrm{otherwise}\end{array}\right.\mathrm{and}\\ & k_{\mathrm{bot}}=\left\{\begin{array}{ll}{k}_{\mathrm{rock}}& \mathrm{if}{k}_{\mathrm{rock}}<k_{\mathrm{b}}\\ k_{\mathrm{b}}& \mathrm{otherwise}.\end{array}\right.\end{array}\end{array}$$

Finally, the retention times τ_srf and τ_blg can be prescribed via namelist parameters since they are dependent on the setup of a given simulation. In case that the model is used in an actual point-scale simulation, the values for the retention times should be set to the time-step length bypassing the intermediary reservoirs altogether. If the spatial resolution of the model is coarse enough that grid cells encompass entire watersheds, the retention times should be similar to the time-scales assumed in the ARNO-scheme, with τ_blg being in the range of 2 to 13 d. For simulations in which grid cells do not encompass entire drainage basins – that is site-level simulations and high- and ultra-high resolution simulations – the choice of retention times is a lot more difficult, but as long as the distances that these retention times represent are in the order of hundreds of metres a retention time of several days still appears very appropriate.