Ground surface temperature variability and permafrost distribution over mountainous terrain in northern Mongolia

ABSTRACT Permafrost plays an important role in numerous environmental processes at high latitudes and in high mountain areas. The identification of mountain permafrost, particularly in the discontinuous permafrost regions, is challenging due to limited data availability and the high spatial variability of controlling factors. This study focuses on mountain permafrost in a data-scarce environment of northern Mongolia, at the interface between the boreal forest and the dry steppe mid-latitudes. In this region, the ground temperature has been increasing continuously since 2011 and has a high spatial variability due to the distribution of incoming solar radiation, as well as seasonal snow and vegetation cover. We analyzed the effect of these controlling factors to understand the climate–permafrost relationship based on in situ observations. Furthermore, mean ground surface temperature (MGST) was calculated at 30-m resolution to predict permafrost distribution. The calculated MGST, with a root mean square error of ±1.4°C, shows permafrost occurrence on north-facing slopes and at higher elevations and absence on south-facing slopes. Borehole temperature data indicate a serious wildfire-induced permafrost degradation in the region; therefore, special attention should be paid to further investigations on ecosystem resilience and climate change mitigation in this region.


Introduction
Both simulations and field investigations at different locations in the Northern Hemisphere report the intensification of enhanced warming and its effects on highand mid-latitude permafrost and frozen ground dynamics (Pepin et al. 2015;Biskaborn et al. 2019). Changes in frozen ground have important impacts on the surface energy balance, hydrological cycle, and carbon exchange between the atmosphere and ground surface (Zhang et al. 2005). Therefore, the response of permafrost to climate change is of major concern, but present knowledge is still limited, especially in remote areas where long-term environmental monitoring is scarce (Azócar, Brenning, and Bodin 2017). In Mongolia, the mean annual air temperature has increased by 2.07°C since 1940, with higher intensification over the mountainous regions (Ministry of Environment and Green Development 2014). Such an amplified warming has been found to increase the active layer thickness, eventually leading to permafrost degradation in northern Mongolia (Ishikawa et al. 2018).
Seasonal freezing and thawing total degree days of air and ground surface (FDD a/s and TDD a/s ) are important factors that reflect changes in the cryosphere (Nelson 2003). These factors have been used to predict seasonally frozen ground and permafrost distribution at regional and global scales (Pepin et al. 2015;Obu et al. 2019). However, for smaller geographic scales, including northern Mongolia, the model outputs have been less accurate due to a lack of precise data sets, such as vegetation and snow cover and microclimatic conditions (Riseborough et al. 2008;Hasler et al. 2015). In mountainous terrain, air temperature, surface characteristics, and radiation balance are major factors in the variability of surface net energy balance (Williams and Smith 1989), as well as the extent of mountain permafrost, especially in semi-arid regions. Several studies have been conducted in the Khuvsgul Mountains of northern Mongolia regarding the distribution of mean ground surface temperature (MGST) to predict the permafrost distribution , active layer thickness , and thawing n-factor (the ratio between air temperature [T a ] and surface temperature [T s ] during the thawing season) in different land covers (Anarmaa et al. 2007).
In regions with rugged topography, such as Sugnugur catchment in the Khentii Mountains, the distribution and dynamics of permafrost have not been well documented. Permafrost in this catchment plays an important role in the hydrological cycle and acts as an impermeable layer allowing lateral subsurface flow that sustains river discharge during the dry seasons (Kopp, Lange, and Menzel 2016). Therefore, it is important to assess the distribution of permafrost with a high spatial resolution. Negative MGST is a good indicator of permafrost presence over mountain terrains where the thermal properties of soil can differ significantly Luo, Jin, and Bense 2019).
This study focuses on understanding the effects of snow and vegetation cover, potential solar radiation (PSR), and topography on the spatial variability of T s and calculating MGST over the Sugnugur catchment to predict permafrost distribution with a 30-m resolution. The calculated MGST is validated against observed in situ T s measurements and borehole data. The results contribute to a better understanding of the temporal and spatial variability of T s and permafrost occurrence in this rugged mountainous terrain of northern Mongolia.

Study area
The study area is the Sugnugur catchment, which is situated in the discontinuous permafrost region of northern Mongolia at the southern boundary of Eurasian permafrost (Figure 1). Thus, it includes the sensitive transition zone between the steppe and the boreal forest (Kopp, Lange, and Menzel 2016). This catchment generates an essential headwater stream of Kharaa River, which eventually drains into Lake Baikal through Selenga River (Menzel, Hofmann, and Ibisch 2011). Therefore, the seasonal frozen ground and permafrost in the catchment are of ecological and socioeconomic importance for the approximately 147,000 inhabitants living in the Kharaa River basin (Karthe et al. 2014).
The elevation of the catchment ranges from 1,062 m. a.s.l. to 2,703 m.a.s.l. at the peak of the Khentii Mountains. Sugnugur catchment has rugged topography, with steep south-facing and relatively gentle north-facing slopes. Mean air temperature (MAT) for the last five years is slightly lower than −1°C and mean annual precipitation is about 350 mm (Kopp, Minderlein, and Menzel 2014). In the study area, north-facing slopes are covered by boreal forest with thick moss layers. South-facing slopes are sparsely vegetated, and the valley bottom is dominated by shrubs and short grasses. Since 2009, about one third of the Figure 1. Location of the study area, its topography, and the positions of temperature measurement transects across the Sugnugur Valley (indicated by T1-T4; see details in Table 1). The inset shows the PPI (at 1-km resolution), which indicates the likelihood of permafrost existence at a global scale (Obu et al. 2019). The hydroclimatic station (green dot) is situated at 1,193 m.a.s.l. BH = the position of the two boreholes located on the north-facing slope near transect 3. forest cover has been affected by wildfires (Kopp, Lange, and Menzel 2016). Permafrost in northern Mongolia mostly exists on north-facing slopes and is absent on south-facing slopes Ishikawa et al. 2008). Depending on the climatic conditions and surface characteristics, it can also sometimes be found in the valley bottoms. Distinct hydrological behaviors of runoff generation processes on burned (quick flow predominant) and unburned (preferential flow predominant) north-facing slopes are evident in this catchment (Lange et al. 2015). The soil in the shrub-dominated bottom of the valley has relatively high moisture contents throughout the summer months due to interflow from the north-facing slopes .

Data
In situ measurements Ongoing hydroclimatic observations have been monitored since 2010 at the bottom of the Sugnugur Valley on a permanent monitoring site (referred to as the hydroclimatic station site in this article; Figure 1). More than forty parameters are continuously measured with high temporal resolutions, including the individual radiation components, snow depth, T s , and T a (for a detailed list of sensors and parameters, see Minderlein and Menzel 2014). Before our monitoring program started, the study region was lacking any environmental data. We calculated the summer temperature lapse rate (LR) for summers 2011 and 2012 using daily mean T a data from three temporary climatic stations distributed along the Sugnugur Valley and the permanent hydroclimatic station. T a was also measured at four locations for 2016-2017 (including complete thawing T a > 0°C and freezing T a < 0°C seasons) using temperature measured from iButtons (Maxim integrated) with ±0.25°C accuracy (Way and Lewkowicz 2018). The iButtons were installed in radiation shields and placed 1.5 m above the surface and recorded temperature at 3-hour intervals. The measurement positions included the bottom of the valley and a mountain summit at 2,020 m.a. s.l. Temperature inversions, which frequently occur during winter, were also studied from the iButton measurements.
We also used iButtons to measure T s at twelve locations at valley transects, in addition to the hydroclimatic station. Measurements were taken at a depth of 5 cm in diverse surface and topographic settings. Two iButtons were operated in parallel at each location in order to prevent any data loss. iButton measurements were evaluated using the T s record with that of a Campbell 107 temperature probe installed at the same depth at the hydroclimatic station site. The iButtons provided slightly higher temperatures with regard to the standard probe, possibly as a result of a greenhouse effect created by the plastic shelters in which the iButtons were housed. The overall Pearson's correlation coefficient between the two separate records was r = 0.997, allowing us to calibrate the iButton data with the Campbell 107 temperature probe using a linear regression.
In addition, two boreholes were dug manually in a severely burned forest (SBF, in December 2015) and in an unburned forest (UBF, in September 2016). The boreholes were each equipped with temperature and moisture sensors at depths ranging from 0 to 290 cm for the SBF site and from 0 to 120 cm for the UBF site. Deeper placement of sensors at the UBF site was not possible because we reached the permafrost table at a depth of 120 cm. Temperature data at 5 cm depth were used as an additional validation of MGST. The relative distance between the boreholes is roughly 600 m and both are located on the same north-facing slope. All temperature sensors were connected with HOBO-U12 data loggers; recordings were made at 4-hour intervals. Both iButton and borehole T s measurement sites are described in Table 1.
Soil thermal conductivity and biomass measurements were also conducted during the thawing season in September 2017 at all iButton and borehole locations. Soil thermal conductivity was measured at depths of 10, 20, and 40 cm using a Decagon KD2 Pro thermal conductivity meter. Biomass samples were collected from 50 × 50 cm plots and weighed and then dried in an oven at 100°C for approximately 2 hours until the weight stabilized. The late winter accumulated snow depth (ASD) was manually measured during a winter field campaign (in the beginning of March 2017).

Remote sensing data
The spatial distribution of vegetation cover was obtained from the soil-adjusted vegetation index (SAVI; Abramov, Gruber, and Gilichinsky 2008;Boeckli et al. 2012), which was derived from six Landsat images (2012)(2013)(2014)(2015)(2016)(2017). These images were taken in either July or August, the months assumed to represent peak vegetation cover (Pan et al. 2019). The snow cover area (SCA) over the Sugnugur catchment was derived from forty remote sensing images from Landsat 7 and 8 and Sentinel 2 satellites . A digital elevation model (DEM) with 30-m spatial resolution from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM 2) was used to derive the spatial variability in T a .

Methodology
This section describes the different processes in the framework of the applied methodology ( Figure 2). First, temperature regime and seasonal n-factors were calculated to understand the permafrost dynamics and the effect of snow and vegetation cover on T s variability (Lunardini 1978;Smith and Riseborough 1996). We then calculated the spatial variability of MGST for 2016-2017, assuming that negative values can indicate the presence of permafrost. MAT for 2016-2017 was derived from the distribution of PSR and annual LR. Surface offset (SO), which is the difference between MAT and MGST, was calculated from in situ temperature measurements in various topographic settings with seasonal snow and vegetation cover and then added to the MAT to obtain MGST based on the derived SCA and SAVI. To reduce unreasonably high or low values, focal statistics were used for smoothing by averaging the values of the surrounding 3 × 3 pixels. The final MGST was validated using iButton surface temperature measurements at twelve locations and temperature data from the two boreholes near transect 3 ( Figure 1) and were also compared to the recently published permafrost probability index (PPI) map by Obu et al. (2019). The equations, parameters, and abbreviations are shown in Table 2.
Temperature variability Temperature regime. To analyze the effects of topography on T s , we placed iButtons along transects from the south-facing to the north-facing slopes of the valley, including the valley bottom. In total, four transects (see Figure 1) were developed at different elevations. The surface characteristics and a description of each temperature measurement site are shown in Table 1. The Figure 2. Flowchart of the methodology.   iButtons were positioned at the same elevation on both slopes in each transect to reduce the effect of the adiabatic T a gradient on T s . T a was observed at the bottom of the valley and on the mountain top, above the T s measurements (Figure 3b). The summer (June-July-August) and winter (December-January-February) LRs were used to derive an annual LR, which was later used as an adjustment factor in the spatial distribution of MAT to 2016-2017. The frost number (FN), which indicates the likelihood of permafrost occurrence, was calculated for all T a measurement locations ( Table 2). The FN ranges from 0 to 1, with a threshold of >0.5 indicating the presence of permafrost at a regional or global scale ). However, it does not confirm the existence of permafrost, because this cannot be determined without knowing the effects of snow and vegetation covers and the thermal properties of the soil (e.g., Smith and Riseborough 2002). To understand the general temperature regime in this region, the TDD s for six seasons and FDD s for seven seasons were also calculated at the hydroclimatic station site.
Effect of surface cover. Seasonal snow and vegetation cover insulates the heat flow between the atmosphere and ground surface, resulting in higher during winter and cooler T s during summer compared to T a . Snow cover decouples the T s from T a due to its low thermal conductivity (Ishikawa 2003;Harris et al. 2009). Using the T s and T a measurements at the hydroclimatic station site, we evaluated the interannual variability of the thawing (n t ) and freezing (n f ) n-factors. Furthermore, the relation between biomass and TDD s was analyzed across the valley transects. The SO was also calculated at four locations where both T s and T a exist. The equations can be found in Table 2. The seasonally persistent SCA-that is, the snow that uninterruptedly covers the ground-was analyzed by determining the Normalized Difference Snow Index with a value of >0.4, indicating the presence of snow cover, as recommended by Hall et al. (2002). Based on the SAVI, the study area was classified into vegetated and nonvegetated surfaces. To simplify the variability in SO over the study area, it was assumed to be different for the vegetated and persistent snow cover areas as well as the nonvegetated, mostly rock debris surfaces on steep slopes based on the calculated n-factors.

Regionalization of air temperature
In mountainous regions with strong topographic gradients, incoming solar radiation plays a significant role in the spatial distribution of both T a and T s . The relationship between PSR and either T s or T a has been studied in permafrost regions (Lacelle et al. 2016). T a in heterogeneous terrain mainly depends on topography and other microclimate factors such as cloud cover, precipitation, and surface characteristics (Strachan and Daly 2017). Temperature and solar radiation follow similar patterns as they change with time, but variances in temperature lag behind variances in solar radiation. In this study, empirically derived MAT was considered as a good indicator of MGST because T a usually follows the diurnal variation in T s due to the reflected radiation from the ground surface (Chou, Sheng, and Zhu 2012). Therefore, we analyzed a five-year time series of incoming solar radiation and T a data recorded at the hydroclimatic station and found that 30 days of lag is the best fit to remove the time lag between the two variables. Similar time lags were also found and used for modeling temperature in previous studies at both local and global scales (Douglass, Blackman, and Knox 2004;Huang et al. 2008). This suggests that if the study area is assumed to have homogenous atmospheric conditions, T a can be predicted as a fraction of PSR after a series of adjustments.
Adjustment of T a includes (1) aggregation of the observed daily values (from transects 2 and 3) to ten-day averages, (2) removal of the variation using discrete Fourier transform (FT), and (3) removal of the thirty-day time lag to improved the correlation with solar radiation. The FT uses two phases, which indicate the absolute values within the annual time series. This results in a smoothing of the T a by sinusoidal oscillations around the mean value. To meet the time step of the adjusted T a , PSR was also calculated for a whole year as a ten-day running average using the PSR model using ArcGIS spatial analysis tools (Fu and Rich 2002) and the DEM. The calculation is the sum of direct and diffuse radiation repeated for each feature location or every location on the topographic surface. Then, MAT for 2016-2017 was obtained from the calculated PSR using a linear regression analysis and adding the annual LR with the assumption that a clear sky was homogenous over the study area. The equation is provided in Table 2.

Temperature regime
Data from the surface temperature iButtons indicate that the freezing season extends from the beginning of October to the end of March at all observed locations, with prolonged periods below 0°C at higher elevations . The surface freezing period on north-facing slopes starts a few days earlier and ends later compared to that on south-facing slopes due to extended snow cover, duration of phase change, surface cover, and less incoming radiation. FN values calculated from the observed T a measurements ranged from 0.5 to 0.6, indicating the potential occurrence of permafrost in the study area (Etzelmüller, Berthling, and Sollid 1998). However, MGST is more sensitive to changes in n-factors, solar radiation, and soil moisture Juliussen and Humlum 2007).
The observed MGST on south-facing slopes was 3-4°C higher than that on north-facing slopes (Figure 3a). The temperature difference between the slopes is more pronounced during the thawing season, especially on hot days with clear skies.
Observations of T a and T s at the four sites indicate that T s is initially warmer than T a at all sites, with notably less fluctuation in winter due to the insulating effect of snow. The difference in T s between the opposite slopes is the lowest during winter. It increases again in late winter because T s on the south-facing slope increases rapidly with T a , whereas T s on the north-facing slope gradually increases due to the extended snow cover and required latent heat energy to warm the soil during phase change. Therefore, the variability in MGST across the valley is high.
The observed MAT at the hydroclimatic station site was −1.06°C. The study period 2016-2017 was unusually warm with low FDD s and high TDD s due to a thicker snowpack in winter and a dry summer, respectively. Moreover, T s observation revealed a continuous increase in both thawing and freezing seasons over the last seven years (Figure 4) as representative of the general temperature regime in the study area. Unusual wildfire events occurred in the region in summer 2017 as a result of the extreme dryness and increased temperatures. Thus, the continuously increasing T s may have contributed to the deepening of the active layer, particularly in areas that were seriously affected by wildfires.
We analyzed the relationship between T s and T a at different elevations using linear regression analysis. The correlation coefficient (r 2 ) between the two variables decreased as elevation increased ( Figure S1, Supplementary material) and can be attributed to the extended snowpack at higher elevations and to slower temperature changes at the ground surface compared to in the air (Peng et al. 2015). Strong seasonal variability in T a exists in the study area as well, as indicated by the mean summer LR of −0.011°C m −1 (found from the four climate stations) and winter inversion of 0.006°C m −1 . Nevertheless, the annual LR was found to be −0.005°C m −1 , or similar to the commonly

Surface covers as insulating materials
To analyze the influence of the snow cover on T s , the snow depth observed at the hydroclimatic station was examined in conjunction with other corresponding parameters ( Table 3). The seasonally persistent snow cover lasts about 130-140 days from early November to late March, and its duration is substantially longer at higher elevations. The mean insulating effect of seasonal snow cover on T s , or n f , was relatively consistent or roughly 0.5 for the last five winters at the hydroclimatic station site (Figure 3).
The development of a seasonal snow cover plays a significant role in the variability in T s ; a gradually weakening insulating effect with increasing snow depth and air temperature toward the end of the freezing season was evident ( Figure 5). During this time, T s became less sensitive to variations in T a due to the greater snow depth. The low T s below the snowpack tended to remain for an extended period as a consequence of the low thermal conductivity of snow restricting the loss of heat from the ground.
In general, snow depth appears to be independent of elevation based on the field measurements (Table S1, Supplementary material). The percentage of persistent snow cover area throughout winter was 94 percent in the Sugnugur catchment . We conclude that n f is homogenously 0.5 in the areas with persistent snow cover. Steep south-facing slopes were snow-free most of the time due to prevailing winds and relatively high incoming solar radiation, which resulted in snow redistribution and greater sublimation rates, respectively.
Forested north-facing slopes were the richest in biomass, followed by the valley bottom shrublands and sparsely vegetated south-facing slopes. Near transect 3, Table 3. Mean and maximum snow depth (cm) at the hydroclimatic station site, together with ground surface temperature (T s ,°C), air temperature (T a ,°C), and incoming solar radiation (ISR, Wm −2 ) in winters since 2012. 2012-2013 2013-2014 2014-2015 2015-2016 2016-2017 2017-2018   for instance, a biomass value of 0.07 kg/m 2 on the south-facing slope, 0.825 kg/m 2 on the bottom of the valley, and a value of 4.935 kg/m 2 under dense forest were found. Notably, the amount of biomass in the SBF was 2.430 kg/m −2 , or roughly half that of the UBF; this distinction demonstrates the effects of wildfires on biomass. Yoshikawa et al. (2002) explained the implications of wildfire and reduction in vegetation cover on the ground thermal regime in Interior Alaska. The soil in the SBF showed the highest thermal conductivity due to the reduction in the moss layer and increased water content, which is similar to that found in Interior Alaska (Yoshikawa et al. 2002). The reduced vegetation cover resulted in a higher TDD s in the SBF compared to the adjacent UBF.
To evaluate the effect of vegetation cover on T s , we correlated the sampled dry biomass measurements against the surface TDD s at all T s measurement sites, except for high-elevation rock debris surfaces ( Figure  S2, Supplementary material). There was a significant negative relationship between the two variables, indicating that greater biomass protects the underlying surface from summer overheating. The lowest TDD s values were observed at the valley bottom site near transect 4, which is characterized by dense shrubs located near the Sugnugur riverbank. We speculate that this is probably due to increased soil moisture content, which reduces surface heating as a result of an interrelation between surface and subsurface flows from the north-facing slope and the shading effect of the dense shrubs.

N-factors and surface offset
The relationship between surface and air temperatures has not been previously described for the Khentii Mountains. The seasonal snow-induced freezing n f was relatively consistent, ranging between 0.5 and 0.6 at all sites, and thawing n t ranged from 0.8 at the top of the mountain to 1.2 in short grassland at lower elevations (Table 4). The n t values exceeding 1.0 indicate that T s was warmer than T a . Unlike at other sites, T s at the rock debris surface was lower than T a , probably due to high albedo resulting in less net radiation. However, other site characteristics such as canopy closure, dryness, and organic layer thickness can also influence this relationship (Karunaratne and Burn 2004). In winter, T s at all locations was higher than T a due to the snow insulation effect.
The period with temperature measurements was not long enough to perform a comprehensive evaluation of the interannual variations in the n-factors. Nevertheless, we believe that seasonal n-factors at different elevations and slopes were identified for this region because the annual variations in the n-factors were relatively stable at the hydroclimatic station.
Considering the similar n f and n t factors except for the rock debris surface, we conclude that SO is +4°C for vegetated and persistent snow cover area and +1.5°C for nonvegetated rock debris surface. These results are in line with Boeckli et al. (2012), who assumed that the SO of the vegetated area is 2°C greater than that of nonvegetated rock debris surfaces in the European Alps, further suggesting that our results are reasonable.
The heat exchange between air and surface varied significantly in early spring and late autumn (the transition periods) when n f and n t converge. Therefore, the freezing season n-factors during the transition periods are important components of the overall seasonal n f . The calculated mid-season n f values were closest to the seasonal n f . Like n f , mid-season n t values were closest to the seasonal n t and during that time T s and T a were similar ( Figure 5). However, the seasonal n t may depend on the length of the thawing season. Figure 5 illustrates the ratios between daily T s and T a in the four transects. Note that calculations were only made for the bottom of the valley and on the mountain top, where both T a and T s measurements exist. There were some days with colder T s than T a during winter at transects 3 and 4 in the upper valley, which was probably due to the inversion effect. Anarmaa et al. (2007) studied n t under various vegetation covers in a valley of the Khuvsgul Mountains in northern Mongolia using a T a measurement from a single location. This would mean that T a is constant across the valley and that the only influence on T s is the overlying vegetation. Such calculation may lead to relatively low n t on north-facing slopes and high n t on south-facing slopes. We suggest that T a can vary at smaller scales depending on both exposure and elevation. Therefore, the variability in T s across the Note: The annual surface offsets at vegetated sites were relatively similar but differed significantly for rock debris surfaces. The values at transect 1 are the mean of the last five years.
Sugnugur Valley was considered as a result of the spatial distribution of snow and vegetation cover, incoming solar radiation, topography, and the length of the thawing season.
In addition, we checked the associations between the calculated parameters and MGST at all T s measurement sites using Pearson's correlation analysis (Table 5). The parameters include the measured snow and biomass, as well as 30-m gridded data including PSR, elevation, and mean land surface temperature (LST) derived from the Landsat images that were used to retrieve SAVI in this study. To avoid high collinearity, the parameters that require T a measurement were not considered due to (1) the limited number of T a measurements in the valley and (2) its strong relationship with T s . The results show that the spatial distribution of MGST might be better explained by the variability in TDD s (r = 0.94, p < .01) rather than FDD s (r = 0.68, p < .05). This indicates that though n f is relatively consistent over the study area, n t plays a major role in the variability in MGST over the Sugnugur catchment. Both the PSR and biomass showed statistically significant correlations with MGST. Due to the high temperature differences across the valley transects during hot summer days, we also tested the relationship between the Landsat-derived LST and the corresponding parameters. These results showed significant negative correlations with biomass and elevation and a positive relationship with TDD s , indicating that the greater the biomass and the higher the elevation, the lower the T s during the thawing season. The ASD did not show any significant relation with its associated parameters.

Derivation of mean air temperature
The MAT for 2016-2017 was derived empirically from the PSR and annual LR. Figure 6a illustrates the tenday running average curves of PSR and T a . Figure 6b shows the linear regression analysis between the two variables and the obtained regression parameters that were used to calculate MAT over the study area. Following the thirty-day time lag adjustment, the r 2 Note: *p < .05. **p < .01. Figure 6. (a) Distribution curves of ten-day average potential radiation and observed and adjusted air temperatures at transect 3. After time lag and FT adjustments, the variance between potential radiation and temperature was significantly reduced. (b) The linear regression on the right shows the correlation between the adjusted T a and potential radiation at transects 2 and 3. The linear equation was later applied for the estimation of MAT for 2016-2017.
increased from 0.76 to 0.96. The result indicates that the distribution of MAT can be improved by considering the influence of incoming solar radiation in addition to the local adiabatic LR. Consequently, the derived MAT showed a significant difference between the north-and south-facing slopes at the same elevations. In a similar landscape, Dashtseren et al. (2014) observed a much lower T a in a forested north-facing slope in both summer (by 2-3°C) and winter (by 1-2°C) seasons, relative to the opposite southfacing slope. However, the linear relationship with PSR may lead to uncertainties on steep slopes with high inclinations where the incoming radiation is very high or low. Therefore, additional T a measurements on different slopes, placed at short horizontal distances to minimize the effect of adiabatic air, may reduce the uncertainties.

Determination of mean ground surface temperature
MGST was estimated for a one-year period from September 2016 to September 2017, including the complete freezing and thawing seasons. Although the study period was short, we believe that the spatial pattern of MGST and permafrost distribution has been successfully derived (Figure 7). The comparison between the produced permafrost map and the recent MODIS-derived PPI using the temperature on top of permafrost model (Obu et al. 2019) shows overall agreement but with significant disagreement particularly on south-facing slopes.
These discrepancies were expected due to the scale difference between the two maps (30 m and 1 km).
Most of the T s measurement sites at each transect, for instance, are within a single 1-km pixel. It is also worth mentioning that the modeled PPI is underestimated in this region because permafrost is covered by forest and moist organic layers, which favor the existence of permafrost (Obu et al. 2019). The findings show good agreement with the T s validation sites (inset in Figure 7), with the exception of a few outliers. In general, south-facing slopes indicate the absence of permafrost (MGST > 0°C), whereas north-facing slopes show the presence of permafrost (MGST < 0°C). The MGST for the study period ranges from +7.5 to −12°C with a root mean square error of ±1.4°C (r 2 = 0.68). The two outliers with underestimations as much as 3°C represent the sites on strongly inclined northfacing slopes. In contrast, the outlier with an overestimation of 2°C comes from the valley bottom at transect 4, where a mixture of shrub and forest persists with shallow snow cover (only 12 cm). The areas in the valley bottom are mostly shaded by dense shrubs and show a relatively high soil moisture content (Kopp, Minderlein, and Menzel 2014). This may result in high evapotranspiration rates, reduced summer heating, and enhanced winter heat loss, thus favoring the existence of permafrost. We believe that such outliers appear due to the empirical linear relationship between PSR and MAT, as well as the shallow ASD and high moisture content. Both the estimated MGST and PPI showed the possible existence of permafrost at the two borehole sites on the north-facing slope near transect 3. The calculated MGST in the SBF was 0.34°C (observed value 1.46°C) with a PPI of 0.41, and that in the UBF was −0.21°C (observed value 0.15°C) with a PPI of 0.44. Borehole temperature data between 2016 and 2018 in these two land covers indicate a significant increase in ground temperature in the SBF ( Figure  S3, Supplementary material). It is important to note that the combination of wildfire and an increased surface temperature in both thawing and freezing seasons over the last seven years may have triggered permafrost degradation in the study area. In addition to the borehole data, electrical resistivity surveys across the valley transects may help to further validate our results and draw a strong conclusion on permafrost distribution.

Conclusions
The distribution of permafrost in mountainous regions is a result of the heat balance between the ground and the air and can be predicted more precisely by determining the effects of vegetation and snow cover, microclimatic conditions, topography, and thermal conductivity of the soil. Based on in situ temperature measurements, we observed an MGST difference of 3-4°C on the north-and south-facing slopes within a short horizontal distance, seasonal T a variability with elevational gradient including winter inversion, and generally increasing T s in both freezing and thawing seasons in the study area. Moreover, the seasonal persistent snow cover insulates the T s from T a by roughly 50 percent as indicated by the average n t of 0.5. The sampled net biomass shows clear spatial variability across the valley transects and has a strong effect on the ground temperature regime and permafrost occurrence. However, short grass and sparse vegetation cover did not show any cooling effect (n t < 1) on T s , probably due to the characteristics of the soil thermal conductivity and the summer overheating at the ground surface.
The overall SOs were roughly 4°C in the vegetated and persistent snow cover area and 1.5°C in the nonvegetated rock debris surfaces. We conclude that MAT can be approximated from PSR and LR after applying several adjustments and contribute to deriving the spatial variability of MGST. Finally, MGST was estimated for 2016-2017 with a 30-m resolution to predict permafrost distribution. The data limitation of the subsurface components did not allow us to derive the thermal offset (difference between MGST and temperature on top of the permafrost) for different soil types, a requirement for determining the temperature on top of the permafrost. The produced MGST indicated permafrost distribution on north-facing slopes and in higher elevations, whereas south-facing slopes were mostly permafrost-free. The result shows good agreement with the most recent PPI map by Obu et al. (2019), except for south-facing slopes. In accordance with both of the maps, the borehole temperature data revealed a presence of permafrost in the UBF and degraded permafrost in the adjacent SBF.
Further investigations will focus on the stability of permafrost and the continuation of both T a and T s in situ measurements. Given that about one third of the forest in this headwater catchment is affected by wildfire, special attention should be paid to further investigations on ecosystem resilience and climate change mitigation in this region. Overall, we conclude that the main controlling factors of MGST and permafrost existence, such as topography, distribution of incoming solar radiation, and seasonal snow and vegetation cover, have been well investigated in this data-scarce mountainous permafrost region.

Author contributions
Munkhdavaa Munkhjargal analyzed the data, carried out the fieldwork, and wrote the article. Gansukh Yadamsuren contributed to the fieldwork and data analysis. Lucas Menzel provided climate data, supported the fieldwork, and supervised the study. Jambaljav Yamkhin carried out the snow measurements and provided crucial suggestions to improve the article.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
Munkhdavaa Munkhjargal received funding from the Federal Ministry of Education and Research (BMBF) through the German Academic Exchange Service (DAAD; award number 57156376) during the preparation of this article. We acknowledge financial support by Deutsche Forschungsgemeinschaft within the funding program Open Access Publishing, by the Baden-Württemberg Ministry of Science, Research and the Arts and by Ruprecht-Karls-Universität Heidelberg. This work was also supported by the Deutscher Akademischer Austauschdienst (57156376).