The surface energy balance in a cold-arid permafrost environment, Ladakh Himalaya, India

Abstract. Cryosphere of the cold-arid trans-Himalayan region is unique with its significant permafrost cover. While the information on the permafrost characteristics and its extent started emerging, the governing energy regimes of this cryosphere region is of particular interest. This paper present the results of Surface Energy Balance (SEB) studies carried out in the upper Ganglass catchment in the Ladakh region of India, which feed directly to the River Indus. The point SEB is estimated using the one-dimensional mode of GEOtop model from 1 September 2015 to 31 August 2017 at 4727 m a.s.l elevation. The model is evaluated using field monitored radiation components, snow depth variations and one-year near-surface ground temperatures and showed good agreement with the respective simulated values. The study site has an air temperature range of −23.7 to 18.1 °C with a mean annual average temperature (MAAT) of −2.5 and ground surface temperature range of −9.8 to 19.1 °C. For the study period, the surface energy balance characteristics of the cold-arid site show that the net radiation was the major component with mean value of 28.9 W m−2 followed by sensible heat flux (13.5 W m−2) and latent heat flux (12.8 W m−2), and the ground heat flux was equal to 0.4 W m−2. The partitioning of energy balance during the study period shows that 47 % of Rn was converted into H, 44 % into LE, 1 % into G and 7 % for melting of seasonal snow. Both the study years experienced distinctly different, low and high snow regime. Key differences due to this snow regime change in surface energy balance characteristics were observed during peak summer (July–August). The latent heat flux was higher (lower) during this period with 39 W m−2 (11 W m−2) during high (low) snow years. The study also shows that the sensible heat flux during the early summer season (May, June) of the high (low) snow was much smaller (higher) −3.4 W m−2 (36.1 W m−2). During the study period, snow cover builds up in the catchment initiated by the last week of December facilitating the ground cooling by almost three months (October to December) of sub-zero temperatures up to −20 °C providing a favourable environment for permafrost. It is observed that the Ladakh region have a very low relative humidity in the range of 43 % as compared to, e.g., ~ 70 % in the Alps facilitating lower incoming longwave radiation and strongly negative net longwave radiation averaging ~ −90 W m−2 compared to −40 W m−2 in the Alps. Hence, the high elevation cold-arid region land surfaces could be overall colder than the locations with more RH such as the Alps. Further, it is apprehended that high incoming shortwave radiation in the region may be facilitating enhanced cooling of wet valley bottom surfaces as a result of stronger evaporation.



Introduction
The Himalayan cryosphere is important for sustaining the flows in the major rivers originating 35 tested (Cao et al., 2019;Fiddes et al., 2015) elsewhere for further investigations in the Ladakh region, where only little data on ground temperatures and permafrost are currently available. It 95 will also help to interpret differences in the relationships of air and shallow ground temperatures (surface offset) observed in Ladakh (Wani et al., 2020) and other permafrost areas (Boeckli et al., 2012;Hasler et al., 2015;PERMOS, 2019).
The specific objectives of this study are to (a) quantify the point Surface Energy Balance (SEB) and its components in a cold-arid Himalayan permafrost environment, (b) evaluate the quality Ladakh is a Union territory of India and has a unique climate, hydrology and landforms. Leh 110 is the district headquarter, where long-term climate data is available (Bhutiyani et al., 2007).
Long-term mean precipitation of Leh (1908Leh ( -2017, 3526 m a.s.l.) is 115 mm (Lone et al., 2019;Thayyen et al., 2013) and the daily minimum and maximum temperatures during the period (2010 to 2012) range between -23.4 to 33.8 °C (Thayyen and Dimri, 2014). The spatial area of  The catchment lies in the Ladakh mountain range and is part of the main Indus river basin.

Meteorological data used 135
The automatic weather station (AWS) in the catchment is located at an elevation of 4727 m a.s.l. at South-Pullu ( Figure 1). It is located in the wide deglaciated valley trending southeast.
The site has a local slope angle of 15°, and the soil is sparsely vegetated. Weather data has been collected by a Sutron automatic weather station from 1 September 2015 to 31 August 2017.
The study years 1 September 2015 to 31 August 2016 and 1 September 2016 to 31 August 140 2017 hereafter in the text will be designated as 2015-16 and 2016-17 respectively. The variables measured include air temperature, relative humidity, wind speed and direction, incoming and outgoing shortwave and longwave radiation and snow depth ( Table 1). The snow depth is measured using a Campbell SR50 sonic ranging sensor with a nominal accuracy of ±1 cm ( Table 1). The analysis of data was performed using R (R Core Team, 2016; Wickham, 145 2016Wickham, 145 , 2017Wickham and Francois, 2016;Wilke, 2019). To reduce the noise of the measured snow depth, a six-hour moving average is applied. Near-surface ground temperature (GST) is measured at a depth of 0.1 m near the AWS using miniature temperature data logger (MTD) manufactured by GeoPrecision GmbH, Germany. GST data was available only from 1 September 2016 to 31 August 2017 and is used for model evaluation, only. All the four solar 150 radiation components, i.e., incoming shortwave (SWin), outgoing shortwave (SWout), incoming longwave (LWin) and outgoing longwave (LWout) radiation were measured. Before using these data in the SEB calculations, necessary corrections were applied Oerlemans and Klok, 2002): (a) all the values of SWin < 5 Wm −2 are set to zero, (b) when SWout > SWin (3 % of data understudy), it indicates that the upward-looking sensor was covered with 155 snow (Oerlemans and Klok, 2002). The SWout can be higher than SWin at high elevation sites such as this one due to high solar zenith angle during the morning and evening hours (Nicholson https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. et al., 2013). In such cases, SWin was corrected by SWout divided by the accumulated albedo, calculated by the ratio of measured SWout and measured SWin for a 24h period (van den Broeke et al., 2004). 160 In high elevation and remote sites, the snowfall measurement is a difficult task with an under 165 catch of 20-50% (Rasmussen et al., 2012). At the South Pullu station, daily precipitation including snow was measured using a non-recording rain gauge. It is a known fact that the snow water equivalent measurements in the mountainous region using collectors have significant errors due to under catch (Yang et al., 1999). In this high elevation area, an under catch of 23% of snowfall was reported earlier  [Unpublished work]. To 170 improve the data quality and match the temporal resolution of precipitation data with other meteorological forcing's, we adopted the method proposed by Mair et al. (2016), called Estimating SOlid and Liquid Precipitation (ESOLIP). This method makes use of snow depth https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. and meteorological observations to estimate the winter precipitation. In this method, the subdaily solid precipitation is estimated in terms of snow water equivalent (SWE). To estimate the 175 SWE of single snowfall events from snow depth measurements, an identification of the snow height increments of the single snowfall events and an accurate estimate of the snow density are necessary. The fresh snow density was estimated based on air temperature and wind speed as below (Jordan et al., 1999):  (Endrizzi et al., 2014) 185 was used for the modelling of point surface energy balance. GEOtop is a physically-based fully distributed model for modelling of water and energy balances at and below the soil surface. It represents the combined ground heat and water balance, the exchange of energy with the atmosphere by taking into consideration the radiative and turbulent heat fluxes. The model also has a multi-layer snowpack and solves the energy and water balance of the snow cover. 190 Furthermore, the temporal evolution of snow depth and its effect on soil temperature are simulated. The GEOtop also simulates the highly non-linear interactions between the water and energy balance during soil freezing and thawing (Dall'Amico et al., 2011;Endrizzi et al., 2014). The model solves Richard's equation in three or one dimensions, and the heat equation in one dimension (1D) (Endrizzi et al., 2014). It can be applied in high mountain regions with 195 complex terrain and makes it possible to account for topographical and other environmental https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License.
variability (Fiddes et al., 2015;Gubler et al., 2013). The model takes into account the effects of complex topography in the estimation of radiation components (Endrizzi et al., 2014), such as (i) the incoming solar radiation is partitioned into direct and diffuse components according to Erbs et al. (1982), (ii) taking into account the solar incidence angle and shadowing of direct 200 incoming solar radiation by topography, (iii) the effects of topography on diffuse radiation coming on the surrounding terrain (Iqbal, 1983).
The model can be operated in two configurations, either in pointwise (1D) or distributed mode (2D) and the processes of interest can be controlled through parameters (Endrizzi et al., 2014).
Generally, the surface energy balance (SEB) (Eq. 3) is written as a combination of net radiation (Rn), sensible (H) and latent heat (LE) flux and heat conduction into the ground or to the snow (G) and must balance at all times (Oke, 2002;Williams and Smith, 1989): where Fsurf is the resulting energy flux at the surface. During the summertime when the melting conditions prevail, the Fsurf is positive and is the energy available for melting snow; otherwise, Fsurf is equal to zero.
But in the GEOtop (Endrizzi et al., 2014) the equations of SEB are described separately. The 220 equations and the key elements of GEOtop are explained in Endrizzi et al. (2014), and here, where , the temperature of the surface, is an unknown in the equation. The Qs is a function of the temperature at the surface ( ), which is an unknown in the equation. Other terms in Eq.
4 which are a function of include , H and LE. In addition, the LE also depends on the 230 soil moisture at the surface ( ), linking the SEB and water balance equations. In Eq. 4, the sum of and is equal to the net radiation ( ) (Oke, 2002). The sign convention adopted is as energy is considered as a gain for the surface or system if Rn is positive and negative for H and LE. Conversely, energy is considered as loss for the surface or system, if Rn is negative and positive for H and LE (Endrizzi, 2007;Oke, 2002). 235 The term in Eq. 4 is equal to the difference between the incoming solar radiation ( ) coming from the atmosphere and the reflected shortwave radiation ( ) (Oke, 2002). The is given by multiplied by the broadband albedo. The (Eq. 5) on a flat ground surface is the product of shortwave radiation at the top of the atmosphere ( ), atmospheric transmissivity ( ) and cloud transmissivity ( ) as follows: 240 The albedo in GEOtop (Endrizzi et al., 2014) is considered as per the ground surface conditions such as, for the snow-free ground, the albedo varies linearly with the water content of the https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License.
topsoil layer, and for snow-covered surfaces, the albedo is estimated according to the Biosphere 245 Atmosphere Transfer Scheme (BATS) (Dickinson et al., 1993).
Also, in Eq. 4 is equal to the difference between the incoming longwave radiation ( ) coming from the atmosphere and the outgoing longwave radiation ( ) radiated by the surface (Oke, 2002). received at the surface is the combined result of the radiations radiated at different heights in the atmosphere with different temperatures and gas 250 concentrations. The clear sky incoming longwave radiation ( , ) is calculated with one of the nine parameterizations present in the model (see Endrizzi et al., 2014). Generally, these parameterisations apply the Stefan-Boltzmann law (Eq. 6), which depends on near-surface air temperature (Ta) in Kelvin, effective atmosphere emissivity ∈ (-) a function of Ta (K) and water vapour pressure (bar), as below: 255 , = ∈ ( , ). .
where is the Stefan-Boltzmann constant equal to 5.67x10 8 Wm −2 K −4 . The parameterisations in GEOtop mainly differ in the value of effective atmosphere emissivity (∈ ). In GEOtop (Endrizzi et al., 2014), the cloudy conditions are treated according to the formulation of 260 (Crawford and Duchon, 1999) as below (Eq. 7): where ∈ is the cloudy sky atmosphere emissivity and is the shortwave radiation cloud transmissivity. The radiated by the surface is also estimated using the Stefan-Boltzmann 265 law (Eq. 8), as below: where is the surface temperature (K) and ∈ is the surface emissivity.
The turbulent fluxes (H and LE) are driven by the gradients of temperature and specific 270 humidity between the air and the surface, and due to turbulence caused by winds as main transfer mechanism in the boundary layer (Endrizzi, 2007). GEOtop estimates the turbulent heat fluxes H (Eq. 9) and LE (Eq. 10) using the flux-gradient relationship (Brutsaert, 1975;Garratt, 1994) as below: where is the air density (kg m -3 ), is the wind speed (m s −1 ), the specific heat at constant pressure (J kg -1 K -1 ), the specific heat of vaporisation (J kg −1 ), and * are the specific humidity of the air (kg kg −1 ) and saturated specific humidity at the surface (kg kg −1 ) respectively, and is the aerodynamic resistance (-). The and are the coefficients 280 that take into account the soil resistance to evaporation, and only depend on the liquid water pressure close to the soil surface. They are calculated according to the parameterisation of Ye and Pielke (1993), which considers evaporation as the sum of the proper evaporation from the surface and diffusion of water vapour in soil pores at greater depths. The aerodynamical resistance is obtained applying the Monin-Obukhov similarity theory (Monin and Obukhov, 285 1954), which requires that values of wind speed, air temperature and specific humidity are available at least at two different heights above the surface.
The input meteorological data required for running the 1D GEOtop model include time series of precipitation, air temperature, relative humidity, wind speed, wind direction and solar radiation components and the description of the topography (slope angle, elevation, aspect 290 angle, and sky view factor) for the simulation point. Also, the latitude and longitude of the https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. study area have to be defined to allow the model to calculate the solar zenith angle, which is important for shadowing estimations.

Heat equation
The equation (Endrizzi et al., 2014) representing the energy balance in a soil volume subject to 295 phase change is given below (Eq. 11): where is the volumetric internal energy of soil (J m −3 ) subject to phase change, t(s) time, ∇· the divergence operator, G the heat conduction flux (W m −2 ), is the energy sink term (W m −3 ), is the mass sink term (s −1 ), (J kg −1 ) the latent heat of fusion, the density of 300 liquid water in soil (kg m −3 ), is the specific thermal capacity of water (J kg −1 K -1 ), T (°C) the soil temperature and (°C) the reference temperature at which the internal energy is calculated. The detailed description of the heat conduction equation used in GEOtop can be found in Endrizzi et al. (2014).

Modelling of snow depth 305
The snow cover buffers the energy exchange between the soil and atmosphere and critically influences the soil thermal regime. In GEOtop, the equations for snow modelling are similar to the ones used for the soil matrix (Endrizzi et al., 2014). The snow processes are solved in a particular order such as (i) solving the heat equation for snow, (ii) metamorphism of the snowpack, (iii) water percolation in the snow and (iv) accumulation due to snow precipitation 310 (Endrizzi et al., 2014).

Model setup and forcing
The 1D GEOtop simulation was carried out at South-Pullu ( Figure 1). The soil column is discretised into 19 layers, with thickness increasing from the surface to the deeper layers. The top 8 layers close to the ground surface were resolved with thicknesses ranging from 0.1 to 1 m, because of the higher temperature and water pressure gradients near the surface (Endrizzi et al., 2014), while the lowest layer is 4.0 m thick.
The snowpack is discretised in 10 layers, which are finer at the top with the atmosphere and at the bottom with the soil, than in the middle. At the top and bottom regions of the snowpack, the vertical gradients are high because of the interactions with the atmosphere at top and the 320 soil at the bottom (Endrizzi et al., 2014).
The model was initialised at a uniform temperature of -0.5 ºC. The soil column in the model is 10 m deep. The model is initialised by repeatedly modelling the soil temperature down to 1 m (2 years*25 times = 50 years), then using the modelled soil temperatures as an initial condition to repeatedly simulate soil temperature down to 10 m (50 years) and finally simulating soil 325 temperatures down to 10 m depth. This initialisation technique has been successfully applied in earlier work (Fiddes et al., 2015;Gubler et al., 2013;Pogliotti, 2011). Preliminary tests show that the minimum number of repetitions required to bring the soil column to equilibrium was 25 ( Figure S1).
In this study, the data recorded by the AWS was used as model forcing, and the forcing data 330 consist of hourly air temperature, wind speed and direction, and global incoming shortwave radiation. The ESOLIP estimated hourly precipitation was also used as forcing to the model.
The model runs at an hourly time step corresponding to the measurement time step of the meteorological data.

Model performance evaluation 335
The evaluation of point SEB was done based on three variables such as radiation components, snow depth and the GST. These variables were chosen because they represent different physical processes influencing the surface energy balance at the ground. For example, (a) the radiation components are the main input driving the surface energy balance, (b) the melt-out date of the snow depth is a good indicator showing how good the surface energy balance is simulated and (c) the GST is the result of all the processes occurring at the ground surface such as radiation, turbulence, latent and sensible heat fluxes (Gubler, 2013). Model performance is evaluated based on the measured and the simulated time series (Gubler et al., 2012). Typically, a variety of statistical measures are used to evaluate the model performance because no single measure enclose all aspects of interest. In this study also, different model evaluation statistics 345 were used for (a) radiation components, and (b) GST and the snow depth as described below.

Performance statistics for evaluation of radiation components
For the evaluation of radiation components, we prefer the statistics mean bias difference (MBD) and the root mean square difference (RMSD) (Badescu et al., 2012;Gubler et al., 2012;Gueymard, 2012). These statistics indicate model prediction accuracy (Stow et al., 2003). The 350 MBD (Eq. 12) is a simple and familiar measure that neglects the magnitude of the errors (i.e. positive errors can compensate for negative ones) (Gubler et al., 2012): Here, is the modelled output variable, and * is the corresponding measured variable. The 355 MBD ranges from -∞ to ∞. The perfect model is the one with an MBD value equal to 0. The RMSD (Eq. 13) is calculated as (Gubler et al., 2012): The RMSD takes into account the average magnitude of the errors and puts weight on larger errors, but neglects the direction of the errors (Gubler et al., 2012). The RMSD ranges from 0 360 to ∞. The perfect model is the one with an RMSD value equal to 0. The formulae (Eq. 12 and 13) used for estimation of MBD and RMSD respectively provide dimensionless quantities https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. since their right-hand side has been divided by the mean of the measured variable (Badescu et al., 2012). Hence, they are expressed in per cent throughout the manuscript for clarity. Also, we use the coefficient of determination (R 2 , Eq. 14) defined as the squared value of the 365 coefficient of correlation which indicates the amount of variation in the modelled variable predictable from the measured variable: 2.5.2 Performance statistics for evaluation of snow depth and near-surface ground temperature 370 For the evaluation of near-surface ground temperature and snow depth apart from the coefficient of determination (R 2 ), different statistical measures have been used such as mean bias (MB), and root means square error (RMSE).
The MB (Eq. 15) provides a good indication of the mean over or underestimate of predictions (Carslaw and Ropkins, 2012). MB is in the same units as the variables being considered. 375 The optimal value of MB is equal to zero. The positive and negative MB values indicate model over-estimation and under-estimation bias, respectively.
The RMSE (Moriasi et al., 2007) is a commonly used statistic that provides a good overall 380 measure of how close modelled values are to predicted values and is given below (Eq. 16):

Meteorological characteristics 385
A summary of the meteorological conditions measured at South-Pullu (4727 m a.s.l.) study site is given in Table 2 to provide an overview of the prevailing weather in the study region. The daily mean air temperature (Ta) throughout the study period varies between -19.5 to 13.1 °C with a mean annual average temperature (MAAT) of -2.5 °C (Figure 2A). The Ta shows significant seasonal variations and instantaneous hourly temperature at the study site range 390 between -23.7 °C in January and 18.1 °C in July. During the two year study period, sub-zero mean monthly temperature prevailed for seven months from October to April in both the years 13.1 and -13.7 °C, for post-winter months (March and April), mean monthly Ta was -5.8 and -8 °C, respectively and for summer months (May to August), the respective monthly mean Ta was 6.6 and 5.5 °C. A sudden change in the mean monthly Ta characterises the onset of a new season, and the most evident inter-season change was found between the winter and summer with a difference of about 16 °C during both the years. 400 The mean daily GST recorded by the logger near the AWS available for one year (1 September 2016 to 31 August 2017) is also plotted along with air temperature (Figure 2A). The mean daily GST ranges from -9.1 to 15.1 °C with mean annual GST of 2.1 °C. The instantaneous hourly GST at the study site range between -9.8 °C in December and 19.1 °C in July. The GST followed the pattern of air temperature, but during winter, the snow cover dampened the 405 pattern. The GST was higher than the Ta except for a short period during snowmelt.
Mean relative humidity (RH) was equal to 43% during the study period ( Figure 2B  The precipitation estimated by the ESOLIP approach at the study site equals 92.2 and 292.5 mm w.e. during the years 2015-16 and 2016-17 respectively. The comparison between observed precipitation (mm w.e.) and the one estimated by the ESOLIP approach is given in (Table S1). 445

Observed radiation components and snow depth
The observed daily mean variability of different components of radiation, albedo and snow depth from 1 September 2015 to 31 August 2017 at South-Pullu (4727 m a.s.l.) is shown in Figure 3. Daily mean SWin varies between 24 and 378 W m -2 (Table 2). Highest hourly instantaneous short wave radiation recorded during the study period was 1358 W m -2 . Such 450 high values of SWin are typical of a high elevation arid-catchment (e.g., MacDonell et al., 2013). Persistent snow cover during the peak winter period for both the years extending from January to March resulted in a strong reflection of SWin radiation ( Figure 3A). During most of the non-free period, mean daily SWout radiation ( Figure 3A) remain more or less stable below 100 W m -2 . Daily mean SWout varies between 2.4 and 262.6 W m -2 with a mean value of 83.3 455 W m -2 ( Table 2). The daily mean LWin shows high variations and ranges between 109 and 345 W m -2 with an average of 220 W m -2 ( Figure 3B, Table 2). Whereas LWout was relatively stable and varied between 211 and 400 W m -2 with an average of 308 W m -2 ( Figure 3B, Table 2).
The LWout shows higher daily fluctuations during the summer months as compared to the core winter months. The daily mean SWn during the study period ranges between 2.5 and 319 W m -460 remain more or less constant with a mean value of -88 W m -2 ( Figure 3C). The mean daily observed Rn ranges from -80.5 to 227.1 W m -2 with a mean of 39.4 W m -2 (Table 2). During 465 both the years 2015-16 and 2016-17, the Rn was high in summer and autumn but low in winter and spring. When the surface was covered with the thick snow during January to early April in 2015-16 and during January to early May in 2016-17, the Rn rapidly declined to low, or even became negative values ( Figure 3D). Albedo (α) is calculated as the ratio of SWout to SWin and is of particular importance in the SEB and in the Earth's radiation balance that dictates the rate 470 of heating of the land surface under different environmental conditions (Strugnell and Lucht, 2001). The daily mean observed α at the study site ranges from 0.1 to 0.48, with a mean value of 0.2 ( Table 2). The daily mean α was low in summer and high in winter and increased significantly when the ground surface was covered with snow ( Figure 3E).

Modelled surface energy balance
The mean daily variability of modelled surface energy balance (SEB) components is shown in Figure 4. The sign convention adopted is as energy is considered as a gain for the surface or 500 system if Rn is positive and negative for H and LE. Conversely, energy is considered as loss for the surface or system, if Rn is negative and positive for H and LE (Endrizzi, 2007;Oke, 2002). The average daily simulated Rn ranges between -35.4 to 136.9 W m -2 with a mean value of 28.9 W m -2 . The Rn shows the seasonal variability and decreases as the ground surface gets covered by seasonal snow cover during wintertime, and increases as the ground surface become 505 snow-free ( Figure 4A The seasonal variation in H points to a larger temperature gradient in summer than in winter. The daily mean LE ranges between -7.2 to 71 W m -2 with a mean value of 12.8 W m -2 . During the snow-free freezing period (October to December) of both the years, the LE decreases (from positive to zero) due to the freezing of moisture content in the soil and also fluctuates close to 515 zero. Furthermore, when the seasonal snow is on the ground, the LE is positive, indicating sublimation and keeps increasing (more positive) after snowmelt indicating evaporation is taking place.
The heat conduction into the ground G remains relatively a smaller component in the SEB ( Figure 4C). The mean daily G ranges between -34.5 to 67.2 W m -2 with a mean value of 0.4 520 W m -2 . The sign of the G, which shifted from positive during summer to negative during winter, is a function of the annual energy cycle. The heat flux available at the surface for melting (Fsurf) https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. ranges between -21.8 to 77.1 W m -2 with a mean value of 2.1 W m -2 (Table 3). During the summer, when snow melting conditions were prevailing, the Fsurf turns positive as a result of energy available for melt ( Figure 4C). Otherwise, the Fsurf become zero. 525    Figure S4). During the representative days of both the years, the significant difference was the changing amplitude of the energy fluxes.
In the 2015-16 year ( Figure S3), the amplitude of Rn and the G during pre-winter, post-winter and summer season were the largest and smallest in winter. The G peaks earlier than those of the LE and H during the pre-winter, post-winter and summer season. The LE and H show strong 545 seasonal characteristics such as (a) during pre-winter season, the magnitude of diurnal variation of H was greater than LE depicting lesser soil moisture content because of freezing conditions at that time, (b) during winter season, the amplitude of LE was slightly greater (sublimation process) than H, (c) during post-winter of low snow year, the amplitude of LE and H was more or less equal and (d) during summer season, the amplitude of LE was greater (higher soil 550 moisture content and evaporation rate) and comparable to H, which is opposite to that of the pattern seen during pre-winter season. However, during winter season the G was close to zero.

Comparison of seasonal distinction of SEB during low and high snow years
The seasonal distinction of modelled SEB components (SWn, LWn, Rn, LE, H, G and Fsurf) for the low and high snow years of the study period is analysed. The seasons were defined as winter (Sep-April) and summer (May-Aug) ( Table 4). These seasons were further divided into two 565 sub-seasons each such as winter but without snow (Sep, Oct, Nov and Dec) and winter with snow (Jan, Feb, Mar and Apr). Similarly, the summer season was divided into two sub-seasons called early summer (May and June; some years with extended snow) and peak summer (July and August).  From the inter-year comparison, it was found that the extended snow sub-season of the 2016-610 17 (high snow year) showed the major differences in energy fluxes between the years.

Model evaluation
In this section, the capability of the 1D GEOtop model to reproduce the point-scale SEB is evaluated. The model was evaluated based on observed radiation components, snow depth and one-year GST. In this study, the simulation results are based on the standard model parameters 615 obtained from the literature (Table 2 and 3, Gubler et al., 2013) and were not improved by trial and error and the same simulation results are used for model evaluation.

Evaluation of radiation components
The first step in our model evaluation was to test the radiation components estimated by the  Furthermore, the simulated melt-out date during the 2016-17 year was slightly delayed by a few days. Overall during the study period, the h was slightly overestimated with MB value of 6 mm. The RMSE value between observed and simulated h was 55 mm. The scatter plot showing the comparison between the simulated and measured h is given in supplementary material ( Figure S5). 650 Furthermore, the performance of the ESOLIP estimated precipitation was evaluated by running the GEOtop model twice, (a) first model run was made with precipitation data measured in the field, and (b) second model run was made with the ESOLIP estimated precipitation as input.
The difference in the simulation of snow melt-out date shows the overall quality of the precipitation data. In comparison to the observed precipitation data as model input, the timing, evolution of snowpack and its melt-out are very well reproduced by the GEOtop model when ESOLIP estimated precipitation was used as input ( Figure 6). Hence, we can say that the ESOLIP is a good approach for precipitation estimation, where snow depth and basic meteorological measurements are available. GST. However, the thermal influence of the winter snowpack was simulated reasonably well by the 1D GEOtop. For example, the start of dampening of temperature fluctuations by the snowpack was reasonably well estimated by the model; however, towards the melt-out of the 675 snowpack, the simulated zero-curtain was little longer as compared to the observation. The possible reason for this behaviour might be due to the estimated precipitation used as input to the model.
In conclusion, given that the GEOtop can simulate the point snow depth and the GST properly, hence, we believe that the model is reliable enough to make robust calculations of SEB in 680 complex topographies. To understand the critical periods of meteorological forcing and its effect on modelled SEB fluxes, we will discuss the diurnal variation of modelled SEB only for one season, i.e., postwinter, which shows major differences in the amplitude of energy fluxes (Figure 8). During 695 the pre-winter, winter and summer seasons ( Figure S3, S4 whereas, all other components (LE, H and G) were of almost zero amplitude ( Figure 8B). The smaller amplitude of LE, H and G is due to the smaller input (solar radiation) and the extended seasonal snow on the ground. Therefore, we can say that the different SEB characteristics during these two years' is a reaction to the forcing of precipitation via snowfall. 705 During the study period, the proportional contribution shows that the net radiation component dominates (80%) the SEB followed by H (9%) and LE fluxes (5%). The G was limited to 5% of the total flux, and 1% was consumed for melting the seasonal snow. The percentages of the energy fluxes were calculated by following the approach of Zhang et al. (2013).  16 and 2016-17, respectively. However, the LE was lower (5%) during a low snow year in comparison to 6% modelled during the high snow year. The H was 9 and 8% during the low and high snow years, respectively. The G was found to be 4.7 and 4.5% for the 2015-16 and 2016-17 years, respectively. During the low and high snow years, the energy utilised for melting seasonal snow was 0.4 and 1.7%, respectively. The mean monthly modelled SEB 720 components for both the years are given in Table S2.
Furthermore, during the study period, the partitioning of energy balance show that 47% (13.5 W m -2 ) of Rn (28.9 W m -2 ) was converted into H, 44% (12.8 W m -2 (Figure 9). In the early winter (September to December), the mean monthly Rn and LE were, respectively, the same for both the years. The main difference can be seen from March to May ( Figure 9A) for Rn, and LE from February to May ( Figure 9B) when the Rn and LE during low snow year were larger than the high snow year. During both the years, the mean monthly H was more or less equal from September to 735 April ( Figure 9C), but during the low snow year, the H was higher for rest of the months.
However, the mean monthly G does not show a clear pattern ( Figure 9D) but the modelled G during both the years follows the pattern of Rn. We found that April, May and June are the most critical months by experiencing higher Rn, LE, and H during low snow year as compared to high snow years ( Figure 9). Also, we found that the mean annual Rn, LE, H and G were 740 larger during low snow years as compared to high snow years.
In both low and high snow years ( Figure 9E, F), the mean monthly H and LE heat fluxes show a seasonal characteristic, such as the H was higher than the LE from September to December for both the years. But for the low snow year, the H was lower than the LE from January to April and after that from May to August again the H was higher than LE. In the high snow 745 year, the LE was higher than the H from January to July and was positive pointing towards the sublimation (January to May) and evaporation processes (June to August).  In early October, the LE began to weaken up to the December for both the years as the seasonally frozen ground began to freeze. Therefore, the seasonal freezing/thawing of the 755 ground affect the LE causing its rapid decrease/increase. Early in October/November, when the seasonally frozen ground began to freeze, the H did not show any significant variability. plateau (Gu et al., 2015;Yao et al., 2011).

Influence of snow cover on near-surface ground temperature
Snow cover affects the ground thermal regime by altering the surface energy balance due to its unique characteristics such as, (a) high albedo, (b) high absorptivity, (c) low thermal conductivity, and (d) high latent heat due to snowmelt that is a heat sink (Goodrich, 1982;765 Gruber, 2005;Zhang, 2005). In comparison to other natural ground surface materials, the low thermal conductivity allows the snow cover to act as an insulator between the atmosphere and the ground. To analyse the effects of snow cover on GST, we plotted the relationship between observed snow depth and GST during the seasonal snow period from 1 January to 23 April 2017 at South-Pullu (4727 m a.s.l.) ( Figure 10). For the shallow snow depth, the GST was 770 smaller, and as the depth of snowpack starts increasing, the GST also starts increasing towards 0 °C. The GST varied from -10 °C to about -2 °C under 40 and 900 mm of snow, respectively.
During the low snow year, the modelled snow depth show no or small insulating effect on GST ( Figure S7). However, during the high snow year, the variations of GST are dampened with increasing snow depth. Furthermore, during the high snow year, only the modelled snow depth 775 greater than 350 mm shows an insulating effect on GST.
The timing of snow cover start and its duration has a non-linear influence on the ground surface temperatures (Bartlett et al., 2004). In the early winter, a thin snow cover can cool the ground, whereas a thick snow cover insulates the ground from cold air temperature variations (Keller and Gubler, 1993). During both the years, the snowfall in the catchment occurred by the last 780 week of December facilitating the ground cooling by almost three months (October to December) of sub-zero temperatures up to -20 °C. This could be a key factor in controlling the https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License. thermal regime of permafrost in the area. Extended snow cover during the high snow year insulates the ground from warmer temperature during May. Therefore, the observed snow cover processes in the study area during these two years have a positive influence on permafrost 785 at the study site. However, long-term studies and larger number of direct snow depth measurements at different locations are required for a better understanding of these processes. Figure 10 Relationship between daily mean observed snow depth (mm) and near-surface ground temperature (° C) from 1 January to 23 April 2017 at South-Pullu (4727 m a.s.l.). The 790 vertical dashed line shows the peak of the snow accumulation, and after that, snowpack starts melting.

Comparison with other environments
In this section, the SEB components from our cold-arid catchment in Ladakh, India are compared with other cryospheric systems, globally (Table 5). Although aiming to represent 795 differing permafrost environments, this comparison also includes SEB studies on glaciers for lack of other data. In most of the studies referred here, the radiation components are measured, and the turbulent (H and LE) and ground (G) heat fluxes are modelled.
The mean SWin measured (210 W m -2 ) in this study site was comparable with values reported from the Tibetan Plateau (200 W m -2 ) but lower than the Andes (239 W m -2 ) and significantly 800 higher than the values reported from other studies such as the Alps (136 W m -2 ). The different surface albedo (α) values help to distinguish the surface characteristics. The mean α for all the sites where radiation balance is measured either on bedrock or tundra vegetation was smaller than those measured over firn or ice during summer with few exceptions. Albedo ranges for glacier ice from 0.49 to 0.69 and for tundra/bedrock from 0.25 to 0.54. The mean LWin (220.4 805 W m -2 ) observed at Leh station is comparable with values observed at Tibetan Plateau (221 W m -2 ) and smaller than the other studies except for Antarctica (184.1 W m -2 ). In the present study, the mean SWn is the largest source of energy gain (127 W m -2 ) and LWn the largest energy loss and strongly negative (-87.6 W m -2 ) and both were higher than in other studies.
However, the SWn (127 W m -2 ) observed at Leh station was comparable with the values 810 observed in the tropical Andes (123 W m -2 ). The mean measured RH (43 %), and the mean annual precipitation are both smaller in this study than in the other areas compared.
Based on the comparison of measured radiation and meteorological variables with other, betterinvestigated regions of the world (Table 5), it was observed that our study area have a very low RH (40% compared to ~70% in the Alps) and cloudiness, leading to (a) Reduced LWin and 815 strongly negative LWn (~90 W m -2 on average, much more than in the Alps). This will lead to surfaces being overall colder than at a similar location with more RH, (b) Increased SWin. This will mean that sun-exposed slopes will receive more radiation and shaded ones less (less diffuse radiation) than in comparable areas, and (c) Increased cooling by stronger evaporation in wet places such as meadows. Where there is enough water, you can cool the ground significantly.   (Bintanja et al., 1997) https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License.

Conclusion
In the high-elevation, cold-arid regions of Ladakh significant areas of permafrost occurrence are highly likely (Wani et al., 2020) and large areas experience deep seasonal freeze-thaw. The present study is aimed at providing first insight on the surface energy balance characteristics 830 of this permafrost environment. The one-dimensional mode of GEOtop model was used to estimate the surface energy balance at South-Pullu (4727 m a.s.l.) in the upper Ganglass catchment from 1 September 2015 to 31 August 2017 using in-situ meteorological data. The model performance was evaluated using measured radiation components, snow depth variations and one-year near-surface ground temperatures which shows good agreement. For 835 the period under study, the surface energy balance characteristics of the cold-arid site in the Indian Himalayan region show that the net radiation was the major component with a daily mean value of 28.9 W m -2 , followed by sensible heat flux (13.5 W m -2 ) and latent heat flux (12.8 W m -2 ), and the daily mean ground heat flux was equal to 0.4 W m -2 . During the study period, the partitioning of surface energy balance show that 47% (13.5 W m -2 ) of Rn (28.9 W 840 m -2 ) was converted into H, 44% (12.8 W m -2 ) into LE, 1% (0.4 W m -2 ) into G and 7% (2.1 W m -2 ) for melting of seasonal snow. Among the two observation years, one was low snow year (maximum snow depth of 258 mm and duration of 120 days) and the another was high (maximum snow depth of 991 mm and duration of 142 days). During these low and high snow years, the energy utilised for melting seasonal snow was 3% and 15%, respectively. Key 845 differences in surface energy balance characteristics were observed during early (May-June) and peak (July-August) summer season of the high snow year. For example, the latent heat flux was higher (39 W m -2 ) during the peak summer of high snow year compared to the low snow year (11 W m -2 ). However, the sensible heat flux during the early summer season of the high snow year was much smaller (-3.4 W m -2 ) compared to the low snow year (36.1 W m -2 ). The 850 diurnal variation of surface energy balance components shows that the extended snowfall https://doi.org/10.5194/tc-2019-286 Preprint. Discussion started: 9 March 2020 c Author(s) 2020. CC BY 4.0 License.
during the high snow year affects surface energy balance characteristics at the study site. The air temperature throughout the study period varies between -23.7 to 18.1 °C with a mean annual average temperature of -2.5 °C, whereas, the near-surface ground temperature ranges from -9.8 to 19.1 °C with a mean annual value of 2.1 °C. During both the low and high years, the 855 snowfall in the catchment occurred by the last week of December facilitating the ground cooling by almost three months (October to December) of sub-zero temperatures up to -20 °C.
The extended snow cover during the high snow year also insulates the ground from warmer temperature until May. Therefore, the late occurrence of snow and extended snow cover could be the key factors in controlling the thermal regime of permafrost in the area. A comparison of 860 observed radiation and meteorological variables with other regions of the world show that the study site/region at Ladakh have a very low relative humidity (RH) in the range of 43% compared to, e.g. ~70% in the Alps. Therefore, rarefied and dry atmosphere of the cold-arid Himalaya could be impacting the energy regime in multiple ways. (a) This results in the reduced amount of incoming longwave radiation and strongly negative net longwave radiation, 865 in the range of -90 W m -2 compared to -40 W m -2 in the Alps and therefore, leading to colder land surfaces as compared to the other mountain environment with higher RH. (b) Higher global shortwave radiation leads to more radiation received by sun-exposed slopes than shaded ones in comparable areas and wet places such as meadows, etc. experience increased cooling as a result of stronger evaporation. However, sun-exposed dry areas could be warmer, leading 870 to significant spatial inhomogeneity in permafrost distribution. The current study gives a firstorder overview of the surface energy balance from the cold-arid Himalaya in the context of permafrost, and we hope this will encourage similar studies at other locations in the region, which would greatly improve the understanding of the climate from the region.