Modeling the Link Between Air Convection and the Occurrence of Short‐Term Permafrost in a Low‐Altitude Cold Talus Slope

We extend a numerical modeling approach developed to explicitly model convective heat transfer in periglacial landforms to represent the ground thermal regime of low‐altitude talus slopes. Our model solves for heat conduction and accounts explicitly for air convection adopting a Darcy term with a Boussinesq approximation for air circulation in the porous ground. Numerical model experiments for the low‐altitude talus slope Dreveneuse, Switzerland, confirm that air convection is the key to forming and maintaining ground ice. In the model, the porous talus slope is underlain by a layer of water‐bearing morainic material. In years, where the gradient between air and talus temperature is sufficiently large to result in increased convection and therefore cooling, ground ice forms due to air convection within the porous material and lasts for more than a year. It is only by considering convection that the model is able to represent the occurrences of ground ice, in accordance with temperature observations on‐site. These findings are important, as they confirm that ground ice can be formed and maintained in landforms with a mean annual air temperature > 0°C if ground air convection is present combined with the presence of water.


| Introduction
Due to their exceptional climatic behavior, low-altitude talus slopes have a long research history.Ground temperatures at the base of such talus slopes are lower than in the surroundings creating a local ground microclimate that is characterized by colder mean temperatures and lower temperature variability.This has an effect on the local flora and fauna, as these low-altitude sites serve as a refuge for different species [1][2][3].Moreover, the existence of ground ice has been proven for low-elevation talus slopes by geophysical surveys [4], which is another factor that favors the occurrence of azonal permafrost in coarse blocky slopes [5].Although their thermal behavior has been extensively studied through field analysis with temperature measurements, geophysical investigations, and geomorphic research [6][7][8][9][10][11][12], attempts of a representation in a numerical modeling framework are sparse [13].In a recent modeling approach, Wicky and Hauck [14,15] developed a model that explicitly accounts for air convection in the ground in the context of coarse blocky alpine permafrost such as rock glaciers and talus slopes.The results showed that ground air convection is a crucial process governing the ground thermal regime and that the model is able to represent it well.Compared to a talus slope within the altitudinal belt of mountain permafrost, convection in low-altitude talus slopes not only influences the distribution of permafrost [16,17] but also is key to the existence of permafrost itself.The mean annual air temperature (MAAT) can be several degrees above the freezing point, and still year-round ground temperatures at or below 0°C can be observed [5,8].Moreover, the temporal dynamics of the ground ice occurrence and its distribution can be quite high in a low-altitude talus slope [8] and cannot exclusively and in general be explained by the incorporation of avalanche snow deposits into the subsurface, as postulated by Kenner et al. [18] for a particular alpine talus slope at significantly higher altitude.The question remains whether the complex processes involved in creating, sustaining, and modifying the thermal regime in undercooled low-elevation talus slopes can be simulated in a physically based model with explicit representation of air convection, that is, that permafrost conditions can be generated and modified by the convective processes themselves.To represent the specific thermal behavior of a low-altitude talus slope, we extend our two-dimensional model approach [14,15] by incorporating a realistic (i) geometry, (ii) subsurface heterogeneity, (iii) boundary conditions, and (iv) additional processes.We constructed a realistic two-dimensional geometry that (i) represents a specific slope profile, accounting for changes in slope angle and micro-topography.In contrast to previous publications that assumed a homogeneous porous medium, we now approximate the subsurface heterogeneities (ii) through the integration of geophysical data.Temperature measurements from different locations along the surface profile are now integrated as spatially heterogeneous boundary condition (iii), which provides a more realistic representation of the model forcing.Furthermore, the model is able to account for additional processes (iv) than in its previous version [14,15], such as the aggregation and melting of ice in porous sediments below the talus (respecting the latent heat exchange), and an observed seasonal asymmetry in the airflow intensity.
We use data from the Dreveneuse talus slope site, which is located below the timberline in the Swiss Prealps and showed the occurrence of permafrost only over some years in the recent past [8,19].The research site offers a broad and high-quality data set that can be efficiently used to constrain numerical modeling.We aim to answer the following research questions: To what extent is the model able to explain the exceptionally cold thermal behavior observed in borehole data through the process of convection?How can ground ice be formed and how can the temporal variability of occurrence and extent be characterized?Finally, the study aims to introduce a new model approach that is capable of combining two-dimensional geophysical data with air convection modeling to simulate the spatial heterogeneity within coarse blocky sediment slopes.

| Description
The Dreveneuse study site is located in the Valais Prealps (Switzerland) above the Rhone Valley on its western flank (coordinates WGS 84 lat/long 46.27344/6.88959,Figure 1A).The studied talus slope lies within a small hanging valley of roughly 2-km length and 800-m width composed of sedimentary rocks, mainly malm limestone.The slope stretches from 1560 to 1650 m.a.s.l, and the slope angle rises from about 30° at the base of the talus to 40° at the top and is oriented to E/ENE (see Figure 2).The MAAT in the years 2008/2009 has been 4.7°C [8].The lowest part is covered by dwarf spruces (mainly Picea abies, Figure 1B) and mosses, and the middle part is covered with mixed forest.The upper part is almost free of vegetation even though it is still well below the treeline (roughly at 2000 m.a.s.l. in the region); limestone blocks (Figure 1C) free of moss and lichens appear at the surface.The Dreveneuse low-altitude talus slope has been extensively studied on a geomorphic field-based approach combined with temperature, anemometer and geophysical measurements [8, [21][22][23], geophysical ERT monitoring [24], and soil moisture characteristics [25].Dreveneuse was part of the Swiss permafrost monitoring network (PERMOS) to represent low-elevation permafrost sites but was later withdrawn from the network because of the absence of permafrost conditions since 2012.A comprehensive overview of the data and site history is given in Permos [19].

| Temperature Data
Ground surface temperature (GST) monitoring data are available from GST loggers distributed along an altitudinal profile over the site (Figure 2), and the loggers with the ID DrB-10, 11, 13, and 14 are used in this study as thermal boundary data (Section 3.4).Logger DrB-12 is dismissed because the data gaps are of too long duration.Logger DrB-15 is located beneath bare blocks, where the energy balance is strongly influenced by the incoming radiation.Therefore, we use logger DrB-14 also for the boundary temperature in the uppermost section of the profile.The ground temperature data used for validation (Section 3.7) is monitored in two boreholes drilled in 2004 [19].Borehole DRE_0104 is reaching a depth of 14.5 m in the central part of the slope, and the underlying morainic material is reached at a depth of roughly 10 m [8].Borehole DRE_0204 in the lower part is reaching a depth of 4.8 m (Figure 2), where the accumulation of porous material is 4 m thick [20].

| Geophysical Surveys
Geophysical surveying at the Dreveneuse site was started in 2007 to characterize the talus slope, but also to monitor potential changes in the subsurface composition related to thermal changes [20].Initial electrical resistivity tomography (ERT) surveying revealed a clearly visible resistive anomaly in the central part of the slope (>50 kΩ m, Figure 3A), which is typical for coarse blocky material and often found in low-altitude and alpine talus slopes [4,6,24,26].The thickness of this resistive layer is around 10-15 m with a more conductive layer (~3 kΩ m) on top and decreasing resistivities below.
Collocated seismic refraction tomography (SRT) measurements from autumn 2020 (Figure 3B) reveal overall low (<1500 m/s) P-wave velocities in the region of the resistive anomaly, which indicate a low probability of substantial ground ice occurrences as the P-wave velocity of ice is markedly higher (~3500 m/s) than either water (1500 m/s) or air (300 m/s).Comparable SRT surveys in alpine talus slopes at higher altitudes usually show much higher P-wave velocities (3000-4000 m/s), which can be interpreted as ice-rich permafrost occurrences (e.g., at the Lapires talus slope; [26]).The absence of P-wave velocities > 2000 m s −1 in the region of the resistive anomaly (cf. Figure 3A) would therefore indicate that no larger ice content volumes have been present in the Dreveneuse talus slope at the time of the measurement in 2020.
A geophysical argument in favor of the former presence of small ground ice volumes until 2015 has been shown and discussed in PERMOS [19] and Mollaret et al. [24], that is, the continuous decrease of mean specific electrical resistivity at around 10-m depth between 2007 and 2015 (Figure 3A).Borehole temperatures in the two boreholes suggest warming at larger depths during this period (Figure 3C), but continuously negative temperatures were only present from 2005 to 2007 and for a short period in 2011/2012 (see Section 4.1).However, the continuous decrease in resistivity between 2007 and 2013 [24] and the geometrical change of the shape of the resistive anomaly between 2007 and 2015, which is still present afterwards (Figure 3A), could also be interpreted as evidence for still ongoing permafrost degradation or at least warming ground conditions.This is, however, in contrast with the measured borehole temperatures, which indicate the absence of ground ice from 2015 onwards (Section 4.1).The interpretation of the geophysical measurements thus leaves open the question of whether small ground ice occurrences (and hence permafrost) still existed in 2015 and whether the (point-scale) borehole temperature data can be considered representative for the whole slope in this respect.

| Seasonal Heat Transfer Cycle
Heat transfer by conduction, as described by Fourier [27], is the main heat transfer mechanism in frozen soils.However, further non-conductive heat transfer processes need also to be taken into consideration [28].Heat can further be transferred through radiation and convection as well as advective water flow.Especially in substrates with high permeability, heat transfer through air convection plays an important role.Convection can be described as heat transfer through a moving fluid.Natural convection in a flat porous ground occurs vertically and results in convection cells (Rayleigh-Bénard circulation), which is also often referred to as the Horton-Rogers-Lapwood problem for a porous medium heated from below [29,30].In mountain slopes, convection develops also a lateral component due to gravity resulting in heat transfer along the slope [31].This is also often referred to as advection.Advection is very pronounced in talus slopes and is referred to as the chimney effect, described by, for example, Wakonigg [32].Heat transfer through advective water flow in talus slopes has also been investigated (e.g., [33]) but is believed to be less important than air convection in coarse blocky terrain, where water drains quickly within the blocky layer [14,15].
The convective air circulation is driven by the temperature gradient between the talus (interior) temperature and the air (exterior) temperature, leading to a seasonally reversing internal air circulation.During winter, the temperature within the talus is warmer than the exterior temperature.Due to the density difference between warm and cold air, a buoyancydriven ascending air flow develops and is released in the upper part of the talus (cf. Figure 1).Snowmelt holes with hoar formations can sometimes be observed at the exits of the warm air [10].The upper part of the talus often shows round melt patterns of air expulsion, especially during very cold periods, and becomes earlier snow-free.As the release of warm air in the upper part of the talus leads to an underpressure in the lower part of the talus, the cold outside air is even aspirated through a thick snow cover, either through snow funnels or through the snow as a porous medium [20].During summer, the circulation reverses.When the air temperature exceeds the talus temperature, cold air sinks down within the talus and is expulsed at its base.The cold air expulsion can easily be observed in the field on warm summer days (e.g., [6]).Figure 4 shows ground temperatures at depth in borehole DRE_0104 during an intensive cooling period from mid-January to mid-February 2011.The temperature drop at the surface is propagated to lower depths not only through the comparatively slow process of conduction, but also through convection.In particular, at a depth of 8.5 and 9.5 m, the fast advection of cold air is visible and demonstrates the strong influence of convective heat transfer on the ground thermal regime of the talus slope.

| Governing Equations
We use a two-dimensional COMSOL-based numerical modeling approach with the same governing equations as described in Wicky and Hauck [15].Heat transfer by convection and conduction is described by Equations ( 1) and ( 2), respectively, where q is the heat flux, Q the storage term, q 0 the boundary flux, k the (effective) thermal conductivity, C p the (effective) volumetric heat capacity, ρ the density, v the velocity, T the temperature, and t the time.
The airflow in the porous domain is described by Darcy's law expressed in Equations ( 3) and ( 4), where Φ is the porosity, Q m the boundary flux, the permeability, μ the dynamic viscosity, g the gravitational acceleration, and p the pressure.Note the g term in Equation ( 4) denotes the Boussinesq approximation accounting for density changes in the air and thus for natural convection effects.

| Geometry
Based on (i) geophysical results showing the different entities of the subsurface along the two-dimensional profile and (ii) borehole data indicating consolidated, probably morainic deposits at a depth of roughly 10 m in the middle part of the talus [8], we (1) | Ground temperatures at different depths in borehole DRE_0104 during an intensive cooling period in January/February 2011.Surface temperature drops strongly around 18th of January and the non-conductive heat transfer is clearly visible through the stronger cooling in the following days, especially at a depth of 8.5 m.Note that almost no cooling is seen below a depth of 10 m, the vertical limit of the coarse blocky material.assume a geometry consisting of three regions.Figure 5 shows the extent of the numerical domain with the three subdomains representing the (1) porous talus, the (2) underlying morainic material, and the (3) bedrock (Table 1).

| Thermal Properties
The thermal properties of the solid part, corresponding to the talus debris and the underlying bedrock, are constant and set to values for Alpine limestone [34] with thermal conductivity k of 2.09 W m −1 K −1 , specific heat capacity C p of 800 J kg −1 K −1 , and density of 2700 kg m −3 .The thermal properties of the air (C p , k, ) in the porous part of the talus (Domain 1, Figure 5) are set to standard literature values as a function of temperature [35].The same applies to the thermal properties of water and ice (C p , k, ) in the morainic material (Domain 2, Figure 5).The change of phase between water and ice accounts for latent heat.The latent energy of freezing (or melting) is represented by a corresponding increase in heat capacity in a transition interval around the phase change temperature of 0°C of [−0.5°C, 0.5°C], where also the thermal properties change linearly.Both the talus and the underlying morainic material are porous media.The porosity of the morainic material is homogeneous and set to 0.2, whereas the talus porosity is linked to the measured normalized resistivity nR (cf.Section 3.3.2),with higher resistivity indicating higher porosity.The thermal properties of the porous domains are obtained through volume averaging.Note that water can only occur in the morainic material, and no lateral or vertical water flow is included.

| Permeability
Airflow in the porous domain is described by Darcy's law where the most important parameter describing the flow is the ground permeability (Equation 4).The permeability is often obtained through empirical relationships, whereby the Kozeny-Carman (K-C) equation is probably the most used one [36].It defines the permeability as where S p is the specific surface area, which is the surface area per unit volume and can also be described with the typical particle diameter d p [37]: However, the K-C relation has been developed with perfect spheres of equal sizes packed in a bed for d p values in the range of mm.For sand, Johansen [37] finds that d p is either the harmonic mean or d 10 , the diameter which 90% of the particles exceed corresponding to, that is, the smaller particle sizes present in the material.For bigger particle sizes, laboratory measurements have been conducted [38,39], as the permeability cannot be measured on-site.These laboratory measurements revealed that the K-C relation is likely to overestimate the permeability in crushed rocks.Furthermore, they derive permeability values of around 0.5-3.9* 10 −6 m 2 for their crushed rocks laboratory experiments, which is consistent with numerical modeling studies on the passive cooling capacities of road embankments [40].The probable ranges of permeability in an alpine permafrost context have been extensively discussed in Wicky and Hauck [15], and model simulations herein were consequently conducted with a range of different permeability values.

| Subsurface Heterogeneity
To account for the spatial heterogeneity of the subsurface, we integrate data from the geophysical measurements conducted along the slope, that is, the electrical resistivities obtained from the ERT survey.We use a single high-quality measurement from 2015-10-21 (Figure 3A) to account for spatial patterns of the heterogeneity in the subsurface.High permeability soils (in a comparatively dry state) tend to show lower conductivity values (= higher resistivity) due to their higher air content, air being an electrical isolator.Under the assumption that resistivity values in the dry, porous talus material are mainly controlled by the porosity (and only to a lesser extent by humidity), this allows to link the permeability to the electrical resistivity measurements as follows [41]: where a is a dimensionless shape factor (ranging from 1.7 to 3) and F is the electrical formation factor, described by Archie's law as follows [41]: where R w is the electrical resistivity of the pore fluid and R t is the total (or bulk) resistivity of the subsurface material, which is measured in ERT surveys and shown in Figure 3A.However, Lesmes and Friedman [41] note that S p is far more sensitive than F. As it is not possible to perform any on-site, we choose to apply different overall permeability values of 1 * 10 −7 , 1 * 10 −6 , and 3 * 10 −6 m 2 that were proven to be adequate for modeling of air convection in coarse blocky alpine permafrost in Wicky and Hauck [15] and that are within the range of proposed values from the literature, whether from laboratory values [38,39] or from modeling studies [31,40].This allows to assess the sensitivity of the results to permeability.On the other hand, Equation (7) shows that permeability is scaled by F, which in turn is defined by the measured two-dimensional resistivity distribution (Equation 8), assuming that R w is constant.The spatial resistivity distribution obtained from the ERT surveys allows us therefore to include a spatial heterogeneity of permeability into the model set-up, which spatially modify the otherwise constant overall permeability values of 1 * 10 −7 , 1 * 10 −6 , and 3 * 10 −6 m 2 .We thus normalize our measured two-dimensional resistivity values of R t from Figure 3A on a logarithmic scale to a scale from 0 to 1 as follows (Equation 9): With nR the resulting resistivity scaling factor (Figure 6), which is then multiplied with a porosity of 0.

| Seasonal Asymmetry
Ground temperature mapping showed an asymmetry between the summer air expulsion area and the winter air aspiration area [21].Further, roughly two times higher air velocities were observed during the winter aspiration phase [8,22] than during the summer expulsion phase, indicating that air convection is more effective during the winter.The exact physical process explaining this spatio-temporal asymmetry is not fully known, but an explanation could be that the summer and winter circulations do not use the same flow paths, a threedimensional effect that cannot be simulated explicitly in our two-dimensional model.
In order to simulate this seasonal asymmetry in the observed air flow velocity, numerous different model experiments were tested.The best results were obtained with the seasonal switch implemented in this study: reducing the permeability during summer to switch off the convective circulation when the outside temperature is warmer than the talus temperature.At each time step, the spatially averaged temperature of the porous talus (Domain 1) is calculated.If the mean talus temperature is lower than the mean air temperature measured on-site at the meteorological station at the corresponding time step, the permeability is set to 1 * 10 −7 m 2 .When the mean talus temperature exceeds the outside temperature and thus the gradient inverses, the permeability increases to either 1 * 10 −6 or 3 * 10 −6 m 2 before being scaled by nR (see Section 3.3.3).The sensitivity of this switch is further discussed in Section 5.2.
Other approaches, such as the cut-off of high summer temperatures, which could potentially be overestimated as a result of radiation effects at the surface, were dismissed.They did not lead to satisfying modeling results in comparison with the validation data.

| Boundary Conditions
The upper thermal boundary condition is of Dirichlet type and is described by measured GST.The loggers with the ID DrB-10, 11, 13, and 14 are used for the simulation of the hydrological years 2006-2018 (cf. Figure 5).The GST forcing data set is almost complete, except for logger DrB-11 for the hydrological year 2014, where data had to be replaced by measurements from the year 2013, which shows a quite similar meteorological behavior.Gap filling algorithms over such an extended period | Section from the model geometry (Figure 2) with nR values taken from the ERT measurements and Equation (9).Parts of the geometry without ERT results are filled with the closest value available.It is important to note that permeability and porosity are only linked to nR in the porous talus (Domain 1) and not in the morainic material (Domain 2).
are not very trustworthy, and therefore, we choose this option.
The measured GST data are prescribed as boundary condition at their respective location at the upper boundary (Figure 5).Daily altitudinal temperature gradients are calculated to yield temperature values for each surface node at the upper boundary between the four GST loggers using a linear interpolation.Below logger DrB-10 and above logger DrB-14, the temperature of the corresponding logger is prescribed.DrB-15 in the upper part was not used because its temperatures are strongly influenced by solar radiation, and therefore, short-term peaks of over 40°C occur in summer.
The lower thermal boundary is of Neumann type and is described by a constant geothermal heat flux of 0.03 W m −2 .The sidewalls of the model are isolated.Pressure at the upper boundary is set to atmospheric using an altitude-dependent pressure head boundary condition (p = g(−y)).Sidewalls and bottom do not allow for any flow.

| Initial Conditions
The initial conditions of the transient model run for the hydrological years 2006-2018 are obtained through a two-step spin-up procedure.The first step consists of a steady-state model run, where the lower boundary is prescribed by the constant geothermal heat flux and the upper boundary condition consists of a temperature of 0°C at the lowermost node with a constant gradient of 0.1°C m −1 .The gradient is chosen to assure that this steady-state run provides stable initial conditions to a transient spin-up.
The lower boundary condition of the transient spin-up is yet again the geothermal heat flux.The upper thermal boundary is now defined by data of the year 2005 from all GST loggers, repeatably applied for a spin-up time of 10 years.Numerical tests showed that quasi-equilibrium is reached after 10 years of spin-up and that repeating the first modeled year (2005) provides stable initial conditions and a smooth transition to transient model simulation, even though 2005 is a comparatively cold year.The spin-up procedure showed a quite high influence on the results and had to be handled with care, especially since the initial ice content has the potential to influence the model results substantially due to the latent heat reservoir.

| Time-Stepping and Mesh Size
The time-stepping follows the adaptive time-stepping scheme proposed by the COMSOL Multiphysics [35] solver.The maximum time step size is 1 day in order to make sure that each data point from the boundary data is modeled.The mesh is composed of triangular elements and consists of a total of 11,137 nodes.The element size is increasing with depth, and the range of mesh size is indicated in Table 1.The comparatively fine mesh in Domain 2 allows for an accurate representation of the airflow.The applied configuration is a result of different mesh tests, revealing the applied mesh size as a good compromise between accuracy and computational effort, whereby a coarser mesh tends to be less accurate with respect to the validation data, and a finer meshing leads to a higher computational cost and at some point to deteriorating results.

| Validation Data
Data for validation are used from the two boreholes on site (cf.Section 2.2).The validation data are thus completely independent from the forcing data.The ground temperature data from the two boreholes DRE_0104 and DRE_0204 are aggregated to daily means to compare measured and modeled data.For borehole DRE_0104, data are missing from 2013-02-27 to 2014-07-16, and for DRE_0204, data are only available until 2016.

| Modeled Subsurface Temperatures
Figures 7 and 8 show modeled and measured temperature at the two boreholes along the profile for the three permeability values of 1 * 10 −7 , 1 * 10 −6 , and 3 * 10 −6 m 2 .At the upper borehole DRE_0104 (Figure 7), it is clearly seen that the smallest permeability value yields far too high temperatures.From this, it becomes clear that convective heat transfer is needed to cool the ground to the observed temperatures, which is not the case for a permeability of 1 * 10 −7 m 2 with conduction dominated heat transfer.Larger permeability values of 1 * 10 −6 and 3 * 10 −6 m 2 lead to results in the range of the measured temperatures.It is important to note that ice aggregation is not included in the porous talus material in the model and, therefore, latent heat effects cannot be represented.This leads to a discrepancy between modeled and observed temperatures at depths of 2.5 and 6 m, where no zero curtain period is simulated in the model.At depths > 10 m, the virtual borehole in the model reaches morainic, water-saturated substrate, but the modeled temperatures are roughly 1 K too high to aggregate ice in the virtual borehole.Further downslope in the model domain, where the temperature is lower, ice aggregates (Figure 9) and allows to represent the observed periods of permafrost.
In the lower borehole DRE_0204 (Figure 8), the missing integration of latent heat becomes apparent.Aggradation of ice in the borehole is also clearly seen by the zero-curtain period, which is absent in the model as porous talus domain is simulated without water in the model.Nevertheless, the subzero temperatures and the observed temperature amplitude, if not masked by latent heat effects, are well captured by the model representing the colder microclimate at the base of the talus slope.

| Temperature and Air Flow Field
Figure 9 shows the simulated circulation and temperature distribution at the end of September just before start of the cooling period (2012-09-30).The 0°C isotherm is shown in blue, indicating that the frozen part is restricted to a small zone in the lower part of the talus, where the model is able to aggregate ground ice in the morainic material.This coincides with the findings from Morard, Delaloye, and Lambiel [8], schematically illustrated in Figure 2. The circulation is a uniform gravitational downflow due to the descending cold air inside the talus.In contrast, Figure 10 shows a situation at the end of January (2012-01-30) where the convective ascending air circulation has been well established.an upwarddirected convective flux can be observed due to the relatively warmer air in the talus.Different convection cells develop, and the location of the cells is hereby partly influenced through the location of the GST Loggers and the associated temperature gradient along the slope.The frozen part at the bottom of the slope is a remnant from the previous year.The (re-) freezing of the morainic material starts about mid-February for a permeability value of 1 * 10 −6 m 2 , whereas higher permeability values lead to stronger and earlier (re-)freezing of the morainic layer as well as more pronounced convection cells.

| Permafrost Occurrence
Figure 11A,B shows the temporal evolution extent of the frozen area in Domain 2 (the morainic material below the porous talus slope) as a function of permeability and porosity.According to the definition, a model region that stays frozen for two or more years can be interpreted as (short-term) permafrost.The size of this frozen area depends on the permeability value of the porous talus above (Figure 11A) and on the porosity of Domain 2 (Figure 11B), which corresponds to the water (respectively, ice) content in the morainic material.The permeability value proved to be a decisive parameter for the presence or absence of permafrost conditions in the simulation.A low permeability value of 1 * 10 −7 m 2 yields almost no freezing in the morainic material underlying the talus slope.The medium value of 1 * 10 −6 m 2 leads to an occurrence of short-term permafrost in certain years.The timing of the occurrence of this short-term permafrost corresponds to the periods of observed freezing in borehole DRE_0104 at a depth of 11.5 m (Figure 7).This is also the depth, where freezing occurs in the model, whereby it is observed slightly more downslope (Figures 9 and 10).A high ground permeability value of 3 * 10 −6 m 2 yields a continuous existence of permafrost in the morainic domain below the talus (Figure 11A).

| Sensitivity
Besides the sensitivity on the permeability values, which is an integral part of the model approach and presented along with the results, two more process-related sensitivities are presented in this section: (i) the influence of preferential flow paths and the resulting seasonal difference in permeability and (ii) the water content in the sub-talus morainic layer that influences the spatio-temporal evolution of the frozen part of the talus.First, Figure 12 shows a reduced data set of a model run activated seasonal model switch.In this configuration, permeability stays constant throughout the year.Consequently, convective airflow acts in the same way as during winter, and modeled temperatures are far too warm from the start onward (the spin-up is performed with the same constant permeability) and never reach temperatures that are close to the observations.As mentioned above, the exact process leading to the observed asymmetry of summer and winter convection is not fully understood, but the high sensitivity in the model to this seasonal switch gives further evidence of its existence.The causes and effects are further discussed in Section 5.2.
Second, the model results are sensitive to the water content in the morainic material below the porous talus slope.The borehole data with the occurrence of permafrost with a long zero curtain period in the temperatures at this depth suggest the existence of ground ice and therefore the presence of water.In Figure 11B, we showed the sensitivity of the results to the porosity of the morainic layer corresponding to the water content.This sensitivity is summarized in Figure 13, showing the minimal, maximal, and mean extent of the frozen surface over the whole modeling period according to the relative permeability and porosity.It is important to note that a higher water content leads to more stable permafrost conditions.The strong cooling in the lower part   of the talus can be more efficiently stored in form latent heat, whereas a smaller water content may lead a bigger extent of permafrost but is temporally less stable as the reservoir of latent heat is smaller.

| Short-Term Permafrost
Our simulations have shown that convection has the potential to cool the ground to an extent that permafrost can develop, but that it can also vanish again in the following years.The occurrence of permafrost is therefore controlled by the complex interaction of the thermal gradient (between air temperature above the ground and talus temperature), the permeability, and the groundwater (respectively, ice) content.Our results for the low-elevation talus slope Dreveneuse clearly show that convection can indeed build up ground ice and thus permafrost conditions, which may persist for some years in the ground, whereas without convection (i.e., for a low-permeability model run), no ground ice is built up.Validation data confirm these findings, but they do only exist pointwise at the borehole locations, which would in general not rule out the possibility that larger and more persistent ground ice occurrences exist in other parts of the slope.In addition, our model approach, even if solving for the full physical process, is still quite conceptual in many aspects.Therefore, we cannot conclude which permeability values can best represent the whole modeling domain and to what extent in time and space we can expect permafrost conditions.
In the context of climate-induced changes of the permafrost dynamics, such temporal and spatial variability of permafrost occurrence is of relevance.Although aggregation of ice due to convective cooling can be studied best close to or even below the lower limit of permafrost (e.g., [4,6]), it is very likely that the process is of importance at any site with a highly permeable ground.At low-altitude sites with a relatively low ground ice content, the model results, supported by observational data, indicate that permafrost degradation is reversible and represents a quite direct response to current cold or warm climatic conditions.The long-term trend is rather described by the frequency of occurrence of cold years than with the evolution of mean climate conditions, the former having decreased in recent years, as no permafrost was observed in the last decade.Therefore, the occurrence of short-term permafrost periods can also be interpreted as a phase of transition towards non-permafrost conditions.However, the high temporal dynamics of these low-altitude convective systems makes it harder to observe a clear long-term trend as can be seen in, for example, in temperature measurements in coarse blocky, but ice-rich, permafrost landforms such as the Murtèl rock glacier [19].
The permafrost distribution and the techniques used for predicting it are affected by this convective process.At the lower margin of permafrost, as observed at Dreveneuse, convective cooling is the tipping point for the existence of permafrost.At higher elevation, where permafrost is more widespread, convection may also play a crucial role in the distribution and the persistence of permafrost [42].Kenner et al. [43] differentiate between ice-poor permafrost, which is mainly the result of the atmospheric forcing, and ice-rich permafrost as a result of mass wasting processes incorporating snow and ice into the landform.This explains why massive ground ice is often found at the base of a slope.Our results confirm that it is likely to find ground ice at the base of slopes, but they also show that this can be the result of air convection in the ground.Unconsolidated debris originating from mass wasting processes are highly permeable porous media favoring the occurrence of air convection in the ground, thus having the potential to control the occurrence and extent of permafrost, at least at its lower margin.

| Seasonal Asymmetry
The exact physical explanation of the lower amplitude of the summer circulation is still lacking, but it has been mentioned in the literature [21,23,44,45].Higher airflow velocities were measured during the autumn aspiration phase than during the summer expulsion phase [22].The asymmetry of winter and summer circulation has also been observed and described in caves [46], although their interpretation is strongly linked to the cave shape and cannot directly be transferred to a talus slope.Further, specifically at Dreveneuse, the location of areas with the observed minimum winter and summer temperatures does not match, pointing to a spatial asymmetry in the convective airflow [21,23].It can thus be assumed that winter and summer air circulation do not use the same flow paths.These temporal differences in flow paths are not represented in the model but approximated by a change in permeability.They can be explained by the nature of the flow, which is directed to the atmosphere when ascending, thereby being less constrained by the talus topography, and directed to the foot of the slope due to gravity when descending.In the latter case, the often slightly concave talus geometry could lead to convergence and preferential flow paths within the talus, thereby channeling the air flow and increasing its strength.This would be in good agreement with observations of predominantly punctual cold air expulsion at the base of the talus, whereas cold air aspiration seems to be more diffuse [45].Other influencing factors for either permeability or seasonal surface boundary conditions include ground heterogeneities that are not captured by the ERT surveys, seasonal variations in snow cover, and ice content within the talus as well as vegetation cover, which all could seasonally seal or alter pore space structures locally.However, this seasonality and its parametrization through a permeability switch are not of relevance for landforms with a slope angle < 15°.In summer, convection takes place only as unicellular gravitative downward flow in the talus and is thus affected by this parametrization.Bories and Combarnous [47] showed in an experimental study, which was later reworked by Guodong et al. [31], that the temperature field from a flat and a sloping layer differs only from an angle of ~15° onwards.This is corresponding to Wicky and Hauck [15], who showed that natural multicellular air convection during winter is fundamental to the understanding of the ground thermal regime of rock glaciers.

| Latent Heat
Latent in the model is represented in the morainic layer, freezing of water takes place and latent heat can be stored and released.This influences the temperature gradients, especially during summer.A prolonged zero-curtain period with the release of latent heat influences the temperature gradient and leads thus to stronger internal air circulation.In this regard, the model has major shortcomings as water is only present in the morainic layer, and all other domains are assumed to be dry.No ice can be built up within the talus, whereas the borehole data show such a seasonal ice formation in this layer.However, it is of a small extent and does not persist throughout the year at borehole DRE_0104, whereas temperatures in the porous domain at the (lower) borehole DRE_0204 remained at or below 0°C for more than 1 year (Figure 8).
In a talus that is ice-free during summer, the energy balance of latent heat within the talus is net-zero over a year.In this case, neglecting latent heat within the talus has only an influence on the convection strength, as the thermal gradient between interior (talus) and exterior (air) temperature depends on the potential existence of ice.Sawada, Ishikawa, and Ono [10] observed the refreezing of meltwater and identified this, besides winter air convection, as the main driving factor of the local cold microclimate of low-elevation talus permafrost.They also stated that the meltwater does not origin locally from the snow on-site but is supplied from groundwater, which is the way the water is represented in our model approach.
A further shortcoming in the model concerns the neglection of incoming convective heat flux through meltwater and liquid precipitation.A previous study from an Alpine site concludes that air convection is of much higher importance for the ground thermal regime than the heat input of meltwater [48], but attention has to be drawn to that issue in future studies.To our knowledge, air convection in the porous ground has so far never been modeled with three phases-rock matrix, water, and air on a two-dimensional spatial scale.Such modeling is challenging, as the boundary conditions are very hard to define because they are usually not fully known (e.g., melt rates, interaction with groundwater, and flow paths of water).

| Conclusions
We applied a physically based numerical model approach in order to better understand and assess the influence of convective cooling in a low-altitude cold talus slope.The numerical model solves for conduction and also explicitly models convention with a Darcy approach in porous media.The following main conclusions can be drawn: • The model approach reproduces the observed seasonally reversing air circulation as a result of the density difference of cold and warm air, leading to a buoyancy-driven convective flow within the talus.The seasonal reversal of the convective flow leads to a pronounced cooling of the lower part of the talus slope.
• Permafrost, respectively, ground ice, can be formed and maintained in sedimentary material beneath a porous talus slope thanks to convective cooling.This is in line with the conclusions drawn from the interpretation of observational data and is now confirmed by the possibility to represent the process in a physics-based model.
• The model results show that the spatio-temporal dynamics of the ground ice depend on the forcing at the upper boundary, the permeability of the porous talus, and thus the strength of the convective heat transfer as well as on the water content in the ground.
• The stronger air circulation during the cooling phase needs to be parameterized, and the full physical process is not known yet.
• The ground heterogeneity was successfully approximated by incorporating electrical resistivity data from ERT surveys.

FIGURE 1 |
FIGURE 1 | (A) Map of the Western part of Switzerland with the location of the study site Dreveneuse (copyright swisstopo).(B) Picture of the dwarf tree area.The black rectangle shows a person for scale (Photo: Ch.Hilbich).(C) Limestone blocks in the lower part of the talus slope (Photo: J. Wicky).

FIGURE 2 |
FIGURE 2 | Overview sketch of the Dreveneuse site with the location of temperature loggers (mean ground temperature in brackets), boreholes, and the start (E1) and end (E48) of the geoelectrical ERT monitoring profile.The arrows mark the winter air circulation situation.([19, 20], modified).

FIGURE 3 |
FIGURE 3 | (A) ERT inversion results for two measurements along the talus slope in October 2007 and 2015.The black vertical lines mark the upper (DRE_0104) and lower (DRE_0204) boreholes.(B) RST inversion results (extent marked in lower panel of A) from September 2020 along the same measurement line.(C) Borehole temperature measurements for the two measurement dates shown in (A).

FIGURE 5 |
FIGURE 5 | Model geometry with three domains: porous talus material (1), a layer of morainic material (2) and bedrock (3), and the meshing.The boreholes are indicated with the purple and pink lines, and the line lengths indicate the borehole depth.The GST loggers defining the forcing at the upper boundary are indicated in green.

FIGURE 7 |
FIGURE 7 | Modeled and measured temperatures at different depths at the location corresponding to the DRE_0104 borehole (Figure 5).Colors represent different ground permeability values in the porous talus (Domain 1).The black line represents the measured data at borehole DRE_0104.

FIGURE 8 |FIGURE 9 |
FIGURE 8 | Modeled and measured different depths in a virtual borehole at the location of borehole DRE_0204.Colors represent different permeability values in the porous Domain 1, and the black line represents the measured data.

FIGURE 10 |
FIGURE 10 | Temperature field (color) and airflow vector field (black arrows) at 2012-01-30 for a permeability value of 1 * 10 −6 m 2 .The °C isotherm is marked as blue contour line showing the extent of the frozen domain for a winter cooling situation with predominantly upwards convective flow.

FIGURE 11 |FIGURE 12 |
FIGURE 11 | Frozen area (m 2 ) in the morainic material beneath the talus (Domain 2).The two panels illustrate the sensitivity of the extent of the frozen area to permeability (A) and porosity (B) and thus water/ice content in Domain 2. The color lines indicate the respective parameter combination.

FIGURE 13
FIGURE 13Mean (dot), minimal, and maximal (range indicated through bars) extents of the frozen area, summarizing the data shown in Figure11A,B.The parameter combination is indicated on the x-axis, respectively.

TABLE 1 |
Summary of the main properties of each domain, depicted in Figure5.