Impact of lateral groundwater flow on hydrothermal conditions of the active layer in a high arctic hillslope setting

Modeling the physical state of permafrost landscapes is a crucial addition to field observations in order to understand its feedback mechanisms within a warming climate. A common hypothesis in permafrost modeling is that vertical heat conduction is most relevant to derive subsurface temperatures. While this approach is mostly applicable to flat landscapes with little topography, landscapes with more topography are subject to lateral flow process as well. With our study, we want to contribute to the growing body of evidence that lateral surfaceand subsurface processes can have a significant impact on permafrost tem5 peratures and active layer properties. We use a numerical model to simulated two idealized hillslopes with inclinations that can be found in Adventdalen, Svalbard, and compare them to a flat control case. We find that ground temperatures within the active layer uphill are generally warmer than downhill in both slopes (up to 1.2◦C in the steep, and 0.7◦C in the medium slope). Further, the slopes are found to be warmer in the uphill section and colder in the very bottom of the slopes compared to the flat control case. As a result, maximum thaw depth increases by about 5 cm from the flat (75 cm) to the steep slope (80 cm), while 10 the medium case does not exhibit a deepening in thaw depth (75 cm). Uphill warming on the slopes is explained by additional energy gain through infiltration and lower evaporation rates due to a overall drier environment. The major governing process causing the cooling on the downslope side is heat loss to the atmosphere through evaporation in summer and enhanced heat loss in winter due to wetter conditions and resulting higher thermal conductivity. On a catchment scale, these results suggest that temperature distributions in hilly terrain can vary considerably compared to flat terrain, which might change the response 15 of subsurface hydrothermal conditions to ongoing climate change.


Introduction
Permafrost is defined as ground that remains below 0 • C for at least two consecutive years. It covers approximately 24% of the exposed land area in the northern hemisphere (Zhang et al., 1999) and stores about 1030 Pg of organic carbon in the upper 3 meters of soil (Hugelius et al., 2014). With increasing air temperatures in the Arctic, this carbon stock gets thawed out of the 20 permafrost, exposing it to microbial decomposition and displacement. How much carbon gets released from the permafrost is strongly influenced by the depth of the active layer, the part of the soil that seasonally thaws out (e.g., Biskaborn et al., 2019).
The correlation between increasing air temperature and depth of the active layer is well established (e.g., Zhang et al., 1997;Isaksen et al., 2007;Frauenfeld et al., 2004). Especially high summer temperatures in dry environments have a direct impact on to the subsurface. Snow and ice are also subject to melting and ponding and can provide a source of water infiltration and/or surface runoff. Unfrozen water flow on the surface follows the Manning equation Painter et al. (2016).
In the subsurface, conductive heat transport follows Fourier's law, with an effective thermal conductivity based on the material properties and accounting for the phase state of the pore-filling fluid (as ice, liquid or air) (Painter, 2011). Advective 95 heat transport occurs as heat carried by water movement in the porous media. Subsurface flow of water is governed by the extended Darcy law for partially saturated flow, where phase transitions follow the Clausius-Clapeyron relationship accounting for latent heat transfer. Soil moisture retention curves, adopting a van Genuchten formulation, are used to describe effective permeability in the variably saturated pore space, accounting for the presence of air and ice, where ice is considered an immobile phase, causing a reduction in available porosity . Furthermore, ATS adopts a flux-conserving 100 finite volume solution scheme and supports unstructured meshes, thus can conveniently be used for applications in 1D, 2D and 3D, accounting for vertical and lateral processes in all dimensions considered.

Field data
The observational weather data to drive the model (hereinafter referred to as the forcing dataset) is derived from a automatic weather station located in Adventdalen (78.2 • N 15.87 • E). The station is operated by the University Center in Svalbard and 105 captures all data needed for the daily surface energy balance, except for precipitation, in the time from 2013 to 2019. Precipitation data was retrieved from the long-term weather station at Longyearbyen airport (78.24 • N 15.51 • E) operated by the Norwegian Meteorological Institute. Mean values for every day of the year (day-of-year average) between 2013 and 2019 were calculated to represent current average weather conditions. Further data processing involved the classification of precipitation as rain if mean daily air temperatures were above 0 • C and as snow if air temperatures were below 0 • C. An adjustment for 110 precipitation undercatch in Svalbard has been suggested to be 1.85 for snow and 1.15 for rain (Førland and Hanssen-Bauer, 2000). Precipitation is multiplied by these factors. The resulting sum of rain (160 mm) and snow (170 mm w.e., total precipitation = 330 mm) was then redistributed to equal daily amounts during the rain-and snow period, respectively. The mean annual air temperature for the calculated averages is -2.8 • C. A representation of all variables within the forcing dataset can be found in the supplementary material (Fig. S1). To inform the model, the same forcing dataset is used for the entire model domain 115 without accounting for temperature lapse rates between the lower and upper part of the transect.

Simulation configurations
Three model cases are considered; a steep case a steep case with 22 • inclination, a moderate case with 11 • inclination, and a flat case with 0 • inclination. The flat case is used primarily as reference to evaluate effects of inclination and to normalise quantities for analysis. The model cases are identical in all respects other than inclination. Note that the elevation difference 120 between the uppermost and lowermost part of the slopes is 10 and 20 m for the medium and steep slope, respectively, but temperature does not change depending on altitude in this setup.
The inclinations are based on slopes as they can be found in Adventdalen and its southern tributaries mostly below 200 m elevation. Geologically, the slopes are located within the Carolinefjellet formation, which mainly consists of shale, siltstone  Areal image of Adventdalen overlooking the same area as in the maps in panels (a) and (b). The picture was taken from a helicopter by A.
Skosgslund (Norwegian polar institute). and sandstone. All hillslope areas greater than 10 • inclination in the area of question are shown in Fig. 1. The flat control case 125 corresponds to areas with no considerable inclination as they can be found in the valley bottom in Adventdalen. These areas are characterized by holocene glaci-fluvial deposits (Norwegian polar institute).
Each model case has a corresponding surface and subsurface mesh. The surface mesh is a 2D layer which extends 50 m in x-direction and 1 m in y-direction, and the subsurface mesh extends 50 m in x-direction, 1 m in y-direction and 20 m in the z-direction (Fig. 2). Both have a lateral resolution of 2 m yielding 25 mesh elements along the x-direction. Only one element 130 with unit width is assigned in the transverse y-direction. Thus, the subsurface elements are 3D volumes and yield volumetric flow quantities, but the model setup effectively represents a 2D transect of the surface-subsurface system with unit width. In   (Schuh et al., 2017). has found to be the temperature at 19 m depth in a borehole in Endalen, one of Adventdalen's tributaries (Hanssen-Bauer et al., 2018). The borehole is located on a slope and therefore assumed to be representative for other slopes in Adventdalen.
As the borehole temperature experiences a linear increasing trend throughout the 2013-2019 period, the mean value of the same period is used. The surface is subject to hydro-meteorological conditions measured on-site (the forcing dataset), which effectively drives the dynamics of heat and water flow through the model system. Precipitation is added as snow and rain 145 on the surface, which allows for infiltration, and heat is supplied by the surface energy balance. Water can leave the system via evapotranspiration and the surface allows for snow and ice accumulation as well as water ponding. Snow distribution is intentionally disabled, in order to yield the same snow accumulation on the surface of the model domain. Model initialisation and spin-up is conducted with a 3-step procedure following previously established routines for permafrosthydrological modelling (Frampton et al., 2013;Karra et al., 2014;Painter et al., 2016). First, a single 1D column is used to establish hydrostatic conditions with the water table at a target depth, using pressure boundary conditions for the top and bottom faces of the model. Second, the soil and water in the column is cooled from below with an assigned sub-zero bottom 165 temperature, until the column is fully frozen and reaches a cryotic steady-state. In the third step, the forcing dataset is used to bring the thermal hydraulic conditions of the column model into an annual steady state. The day-of-year average yearly cycle created from the weather data is repeated for 50 years to create the forcing dataset. After 50 years, this yields an inter-annual temperature differences throughout the column of less than 0.01 • C.
This above procedure is necessary to obtain a physically consistent system which can be used as initial condition for main between the different inclination cases, and there are also temperature differences between the uphill and downhill observation locations. Additionally, timing of thaw and freeze-up varies between cases. To enable a systematic study of the impact of the different hillslope inclinations, we consider daily temperature differences ∆T I between the steep slope and flat case (steep-flat), as well as between the medium slope and flat case (medium-flat). We also consider daily temperature differences ∆T E between uphill and downhill observation points (uphill-downhill), corresponding to different elevations along a hillslope (Fig. 3).  There is significant variability in these temperature differences over the year, with most pronounced differences occurring during the warm season, typically including a peak at onset of thaw and another peak before freeze-up, indicating greatest differences occurring during these times. Between the uphill and downhill side in the steep and medium slope (Figures 3a,b), it can be seen that the uphill side is warmer than the downhill side throughout the year (positive ∆T E ), with two short exceptions during thaw and freeze-up (negative ∆T E ). The warming is strongest in summer with two peaks developing just after thaw and  Temperature differences ∆T I between the steep slope and the flat case (Fig. 3c,d) and between the medium slope and the flat 195 case (Fig. 3e,f) show that on the uphill side, the slopes are warmer in summer, colder during freeze-up and very similar to the flat case in winter. On the downhill side, the slopes are slightly colder than the flat case in winter and significantly colder during summer and just after freeze-up. During freeze-up however, the downhill sides are slightly warmer than the flat reference case.
An overview of the yearly maximum differences is given in Tables 2 and 3 (Table 2). Also, the steep slope case exhibits greater differences than the medium slope case for all depths in the active layer.
Note however that the differences for the steep slope case are not quite double that of the medium slope case, despite the inclination being doubled.
Considering the temperature difference between slope inclination cases (Table 3) This inversion of temperature difference indicates that different processes are responsible for these differences in subsurface temperatures at various depths.

Spatial analysis of ground temperatures
The greatest difference in temperature along the subsurface transect is between the two outermost columns (at x=0 m and 210 x=50 m), corresponding to the two locations farthest apart along the hillslope. To better visualise the ground temperature differences between cases throughout the subsurface domain, the temperature difference between the steep and the flat case (Fig 4a), and the medium and flat case (Fig 4b) in the upper 1.2 m of the subsurface are considered. Note that Fig 4 shows cell-based temperature differences between cases; thus slope inclination is not depicted.
The upper three panels in each figure show snapshots of temperature differences during thaw (June) and summer (July, Au-215 gust), and the lower three panels show temperature differences during freeze-up (October, November) and winter (December).
In both cases, the dates are separated by 20 days. For each day, the 0, • C isotherm(s) from the steep and medium case respectively is (are) represented as black dotted line(s). The average temperature in this volume of the subsurface (upper 1.2 m) is given in Table 4.
Ground temperatures in the sloped cases are generally warmer than in the flat case during thaw and summer (red shades).

220
The temperature differences are greatest near the progressing thaw front, i.e. near the 0 • C isotherm, as well as on the uphill side (x=0 m), but a gradual change towards same temperatures as the flat case (red to white) can be observed in the lateral direction  In the medium slope case, only the first 10 m (x=0-10 m) of the slope show warmer temperatures in the active layer on that day, compared to the rest of the slope. In both cases, it can be seen that the last column (x=50 m) remains significantly cooler 235 than the rest of the slope, causing the active layer to be shallower there.  To highlight differences in the developing thaw depth between cases, spatially averaged thaw depths over the entire transect for each case are calculated (Fig. 5a). As noted above, the steep slope has a warming effect on the uphill section of the transect, which causes its maximum thaw depth (Fig. 5a, dark blue line) to be greater; spatial mean active layer depth on the date of maximum active layer depth is -0.77 m. On the other hand, the medium slope (cyan line) does not develop a substantially deeper 240 maximum active layer thickness than the flat case, as the temperature difference is not enough to increase it ( -0.73 m maximum thaw depth in both cases). However, it can be seen that the medium slope case reaches maximum active layer thickness slightly earlier than the flat case (yellow) and freezes up slower, reflecting the slightly warmer temperatures.
Three snapshots of the thaw progression throughout the subsurface during the summer and freeze-up seasons are shown in -0.55 m) for lateral distances between x=10-48 m (5b, solid lines). Between x =0 to 10 m the thaw depth in the medium case in July is deeper and the last two meters (x =48 to 50 m) shows a shallower thaw depth as compared to the flat case.
This observation on thaw depths for July 20 is consistent with the previously discussed temperature differences in Figure 4; temperatures are generally warmer uphill and cooler downhill for the sloped cases, and the cooling effect is primarily observed A further observation is that two-sided freezing occurs. This is evident by the remaining thaw bulb with encroaching freezing front both from above as well as below (dotted lines). Subsequently, the active layer is fully frozen again on November 20 for the steep slope case, on November 19 for the medium slope case, and finally on November 18 for the flat case.

265
In summary, we observe that the steep slope case has a notable influence on thaw propagation and active layer thickness, which we attribute to an increase in ground temperatures compared to the flat case, observed primarily in the center-uphill active layer is not necessarily deeper in the medium slope. Hence, the warming effect due to slope inclination does not only play a role in the vertical soil profile, but also in the timing of freeze and thaw.

Saturation and thermal conductivity
Due to lateral gravitational water flow during the warm period, moisture accumulates on the downhill side, and yields greater liquid saturation there when compared against the flat reference case which is not subject to lateral flow (Fig. 6, first column).

280
This leads subsequently to greater ice saturation on the downhill side during the frozen period ( Fig. 6 second column). The increased ice and liquid saturation in the slopes, and consequently reduced air saturation (Fig. 6 third column), results in a considerably greater effective thermal conductivity during winter and slightly greater effective thermal conductivity during summer ( Fig. 6 fourth column). Considering the little snow cover in winter (max. 0.01 m, see Fig. S2), the effect should be an enhanced heat loss (cooling of the ground) during winter, and slightly enhanced heat gain (warming of ground) during summer, 285 when compared against the flat reference case.
Recall the previous discussion on temperature differences between the sloped and flat cases (Fig. 3d,f). The downhill side of the sloped cases experience cooler winter temperatures, especially shortly after freeze-up in late November. This is consistent with the differences in effective thermal conductivity, as an increased thermal conductivity during cold periods enables an enhanced ground heat loss, yielding cooler winter ground temperatures. However, summer temperatures also exhibit significant  cooling when compared to the flat case (Fig. 3d,f). This is not consistent with the increased effective thermal conductivity summertime, as it should enhance heat uptake to the ground, leading to warmer ground temperatures. Thus we conclude changes in effective thermal conductivity alone does not suffice to explain the temperature changes incurred on the downhill side of the domain for the two hillslope cases.
Next, consider the uphill side of the domain for the hillslope cases. They are slightly drier at depths 0.2 m, 0.4 m and 0.75 m, 295 both for summer with less liquid saturation, and winter with less ice saturation (Fig. 6, first and second columns, respectively).
This slightly reduces effective thermal conductivity with respect to the flat case at those depths, mainly in winter and slightly discernible also in summer (Fig. 6, fourth column). Thus, when compared against the flat reference case, the uphill side of the inclined cases should exhibit warmer ground temperatures during winter due to reduced thermal conductivity (greater

Combined mass and energy fluxes
To further understand how much energy is carried by seeping water, we compare the lateral advected energy flux on the left or 360 right faces of the CVs alongside the water mass flux on the same faces, and compare the timing of peaks ( Fig. 9; a complete presentation of the mass fluxes across all faces is provided in Fig. S4 in the supplementary material). Note that units between advective heat flux and mass flux are different and that the following interpretation focuses on the shape of the curves, rather than absolute values. As can be seen (Fig. 9), advective heat flux (blue) peaks before September in both uphill and downhill CVs and declines shortly after. Mass flux (yellow) also has its first peak before September, but with prolonged duration of flow 365 and declines gradually.
For the uphill CV, it can be seen that advective heat flux is close to zero already by the end of September, while mass flux reaches zero only by mid November. The downhill CVs exhibit a second, less distinct mass flux peak just before and during freeze-up in the end of October, which is however not associated with a peak in advective heat flux. This indicates that heat is being carried with water flow during the warm season, corresponding to mid-thaw period, but little advective heat is being 370 transported by the end of the thaw season. This is caused by the permafrost acting as a significant heat sink and reservoir for cooling of the soil column above. Infiltrating water from the surface gets cooled down rapidly causing it to attain equilibrium with its surroundings. Then, although water seepage and flow occurs, it does not contribute much to advective heat transport, as the flowing water is at the same temperature as its surroundings. lateral cryosuction.

Impact of changes in precipitation
In order to investigate the effects of reduced or increased precipitation, two additional wetness scenarios are considered for each hillsope; a dry scenario (S0R0) and scenario with increased wetness (S2R2). Snow (S) and rain (R) precipitation rates are set to 0 for S0R0, resulting in a completely dry climate, and the rates are multiplied by 2 in the S2R2 scenario, resulting in 380 a climate that is twice as wet as the current climate. We compare the scenarios with regard to temperature differences, active layer thickness, timing of freeze-up, and advective energy flux.
Firstly, we find that both slopes and the flat case are notably warmer in the no-precipitation scenario and colder in the doubled precipitation scenario. Relative temperature differences between the slopes are generally in a similar range as in the original precipitation scenario with the exception of the downhill side (last column of the domain) in the S2R2 scenario being 385 even colder. On October 28, the downhill side is as much as -1.3 • C colder in the steep case compared to the flat case, while the original precipitation scenario showed a maximum difference of -0.25 • C on the same day. In the dry S0R0 scenario the difference is only -0.09 • C on the same day between the steep and the flat case. A similar pattern can be seen consistently throughout the year.
Active layer thickness further support these findings. Maximum active layer thickness is deepest in the scenario with no  Fig. S5 and Fig. S6). Note that the difference in absolute maximum active layer thickness between the medium and flat slope is very small averaged throughout the transect. Due to the temperature difference, however, the medium case experiences an earlier thaw and delayed freeze-up in the sensitivity scenarios as well as in the original scenario.

395
The timing of freeze-up is different throughout the inclinations in each scenario. In the original scenario, all cases are fully frozen again on November 20. While the scenario with no precipitation (S0R0) also shows the latest freeze-up on November 20, in S2R2 the last day with unfrozen subsurface-cells is November 3, about two weeks earlier than in the other two scenarios.
Overall, the scenarios show that the higher the amount of recharge added through precipitation on the surface, the lower the ground temperatures will be in the both the sloped cases as well as the flat case. Note that multiplying snow by a factor of 2 400 still did not result in a snow cover significant enough to have an insulating effect on the subsurface.
These results are consistent with previously observed cooling effect of precipitation on the active layer. Wen et al. (2014) and Wu and Zhang (2008) both documented a cooling of the active layer in response to rainfall on the Tibetan Plateau. In contrast, e.g., Douglas et al. (2020) and Mekonnen et al. (2021) found a warming effect of summer precipitation on active layer temperatures. However, those studies do not account for the influences of topography.

405
As for heat fluxes along the faces of the CVs, we find that most fluxes follow the overall observed patterns in the original simulation depending on moisture content. For instance, diffusive heat flux in summer near the surface is lowest in the dry case S0R0 due to lower saturation, highest in the wetter case S2R2, and with the original precipitation scenario located in between.