Meteorological impacts of a novel debris-covered glacier category in a regional climate model across a Himalayan catchment

Many of the glaciers in the Nepalese Himalaya are partially covered in a layer of loose rock known as debris cover. In the Dudh Koshi River Basin, Nepal, approximately 25% of glaciers are debris-covered. Debris-covered glaciers have been shown to have a substantial impact on near-surface meteorological variables and the surface energy balance, in comparison to clean-ice glaciers. The Weather Research and Forecasting (WRF) model is often used for high-resolution weather and climate modelling, however representation of debris-covered glaciers is not included in the standard land cover and soil categories. Here we include a simple representation of thick debris-covered glaciers in the WRF model, and investigate the impact on the near-surface atmosphere over the Dudh Koshi River Basin for July 2013. Inclusion of this new category is found to improve the model representation of near-surface temperature and relative humidity, in comparison with a simulation using the default category of clean-ice glaciers, when compared to observations. The addition of the new debris-cover category in the model warms the near-surface air over the debris-covered portion of the glacier, and the wind continues further up the valley, compared to the simulation using clean-ice. This has consequent included in high-resolution regional climate models over debris-covered glaciers.

included in high-resolution regional climate models over debris-covered glaciers.

K E Y W O R D S
debris-covered glaciers, Himalaya, mountain meteorology, WRF

| INTRODUCTION
Debris-covered glaciers consist of ice covered by a layer of loose rock, over the full width and at least part of the length of the glacier (Kirkbride, 2011;Anderson and Anderson, 2016). In the Himalaya, 16.9% of the glaciers are debris-covered, and this is set to increase with future glacier melt (RGI Consortium, 2015;Scherler et al., 2018). Much of this debris cover is thick debris, for example on Lirung Glacier and Ngozumpa Glacier, both in Nepal, debris thickness is measured up to 2.3 m and 8 m, respectively (McCarthy et al., 2017;Nicholson et al., 2018). Correctly simulating the meteorology above debris-covered glaciers is necessary for accurate inputs into hydrological and glaciological models, and for local-scale forecasting, which are essential for predicting future energy and water security in the Himalaya (Eriksson et al., 2009).
The substantial and varying effects of glacier debris cover on near-surface meteorological variables have been recognized in recent years through a series of observational studies. Higher surface temperatures, and a larger amplitude in the diurnal cycle of near-surface temperature have been found over debris-covered glaciers, compared to clean-ice glaciers (Steiner and Pellicciotti, 2016;Yang et al., 2017;Nicholson and Stiperski, 2020). It has been shown that after heating during the day, debris cover greater than 0.05 m returns the energy to the atmosphere at night, rather than heating the ice below (Reznichenko et al., 2010;Nicholson and Benn, 2013). Differences have also been found in near-surface wind speed and incoming short wave radiation between debris-covered and clean-ice glaciers (Yang et al., 2017). High-resolution turbulence-resolving models have also been used to find that increasing roughness lengths over debris-covered glaciers increases spatial heterogeneity of the surface heat flux (Bonekamp et al., 2019), but such models are not suitable for regional climate modelling. Observational and turbulence-resolving modelling studies such as these are vital to give an indication of how debris cover affects near-surface meteorology on different glaciers. However, regional atmospheric models are also needed, as they allow for a more direct comparison between the effects of debris-covered and clean-ice glaciers on the atmosphere than observations, through sensitivity experiments. Regional models are also useful to compare a wider range of meteorological and surface parameters than are available from observations, and over a much wider area.
There has only been one attempt so far to include a representation of debris-covered glaciers in a regional atmospheric model . This was introduced as part of a coupled atmosphere-glacier model, using the Weather Research and Forecasting (WRF) model for the atmosphere modelling. As such, the surface fluxes were computed by the glacier model (Collier et al., 2014), and debris characteristics such as thickness were needed as a model input. The model was tested over the Karakoram and compared to satellite data, and  found that near-surface air temperatures were substantially improved, and up to 20 C higher at some elevations in the model run with the addition of debris cover, compared to a corresponding clean-ice glacier model run. In addition, wind speed was reduced and the daytime wind direction changed from katabatic in the clean-ice model run to anabatic in the model run with debris cover. Coupled models are considerably more computationally expensive to run than atmosphere-only models, and therefore in addition to the model proposed by , there is a need for a simple representation of debris-cover within the existing land and soil categories in the WRF model. This study presents a method for adding a simple representation of thick debris-cover to the WRF model with the Noah Land Surface Model used as the land surface scheme. The model results are compared to those from a model with clean-ice glaciers, against observations at a weather station on Changri Nup Glacier (in the Dudh Koshi River Basin, Nepal). Furthermore, the two model runs are compared over the entire Dudh Koshi River Basin, to establish more generally the effects of adding debris-cover to the WRF model on the near-surface meteorology.

| METHODS
In this study, the WRF model is run at 1 km resolution with hourly output, for July 2013 over the Dudh Koshi River Basin, Nepal (Figure 1). The model is run twice, with and without debris-cover over part of the glacierized region. The first run uses the default clean-ice glacier mask ("clean-ice run"), and the second includes the new debris-cover land cover and soil category ("debris-cover run") as described below. The model output is compared to near-surface observations of air temperature, wind speed, and relative humidity from an automatic weather station (AWS) on the debris-covered area of Changri Nup glacier, located at 27.98 N, 86.87 E ( Figure 1). The data at Changri Nup were collected by the glacioclim group (Université Grenoble Alpes), and more information can be found at https://glacioclim.osug.fr/spip.php?article75. Observations are missing from the AWS for each morning, from about 02:00 to 07:00 LT.
To compare the model output against the data at the AWS at Changri Nup, the nearest model grid point to the AWS is taken (less than 500 m away from the AWS both latitudinally and longitudinally). The model temperature at 2 m is adjusted to account for the difference between the model elevation (5,432 m asl [meters above sea level]) and actual AWS elevation (5,363 m asl), using the average 2 m air temperature lapse rate calculated over the entire model domain at each hour (the slope of a linear regression between air temperature at 2 m and elevation). The observed wind is adjusted from the observational height of around 2 m to the model height of 10 m using a logarithmic profile, using a roughness length of 0.02 m and assuming neutral stability (Whiteman, 2000). This increases the maximum average daily wind speed by approximately 0.5 ms −1 .

| Atmospheric model description
The WRF model version 3.8.1 (Skamarock et al., 2008) is run with four nested domains, with resolutions of 27, 9, 3 and 1 km ( Figure 1). There are 50 vertical levels in all model domains, from the surface to 50 hPa. The topographical input comes from the Shuttle Radar Topography Mission (Jarvis et al., 2008), at 90 m resolution. The permanent snow and ice (i.e., the glaciers) in the WRF model have been updated using the Randolph Glacier Inventory (RGI) (RGI Consortium, 2015), as suggested by , in both the clean-ice run and the debris-cover run. At present, the RGI does not record debris-covered glacier outlines separately from clean-ice glaciers. The outer domain of the model is forced using ERA-Interim at the lateral boundary (Dee et al., 2011), every 6 hr. Both model runs are run from 1st to 31st July 2013, with a 2-week spin-up period from 16th June. The model setup is the same as that described in Potter et al. (2018), with the exception that the land surface scheme is changed to the Noah LSM (land surface model) scheme (Chen and Dudhia, 2001), instead of the Noah MP (Multi-Physics) scheme. While Noah MP provides a better representation of snow and ice processes (Niu et al., 2011), Noah LSM consists of a simpler framework in which to implement the debriscover addition. The full model setup can be found in Table S1. The default landcover categories used in this study are provided by the U.S. Geological Survey (USGS) dataset, which contains 24 land cover categories (grassland, mixed forest, permanent snow or ice ["glaciers", in this instance], etc.). There are also 16 soil/rock texture categories (silt, clay, sand, etc.). In the Noah LSM scheme in the WRF model, values for surface parameters such as albedo, surface roughness, and vegetation parameters are assigned to each of the land cover categories in the VEGPARM.TBL table and fixed for summer and winter. These are then modified by the Noah LSM scheme based on factors such as snowfall during the model run. In addition, subsurface soil texture properties such as porosity and hydraulic conductivity are assigned through the SOILPARM.TBL table (note that here "soil" includes any subsurface properties, including bedrock and debris). In the Noah LSM scheme, there is a separate regime for calculating surface fluxes over glaciers (for example, it stipulates a maximum surface temperature of 0 C over the glacier), which is switched off for the debris-cover addition. Different land surface schemes in the WRF model calculate surface fluxes in different ways, and may not use the VEGPARM.TBL and SOIPARM.TBL values. For example, the method used here, which is applicable to the Noah LSM scheme, is not appropriate for use with the Noah-MP scheme.
For this simple representation of debris cover in the WRF model, the debris is assumed to occupy the full soil depth in the WRF model, which extends to 255 cm below the surface (i.e., there is no representation of the ice underneath). As such, this method best represents thick debris cover, as in general thick debris cover is thought to provide insulation between the air and the glacier ice below the debris cover (Reznichenko et al., 2010;Nicholson and Benn, 2013).
As such, there are two stages necessary to add a debris-covered glacier category to the WRF model. First, a debris mask is created to input as a land cover and soil category, which is added to the WRF Preprocessing System (WPS). Second, values are chosen for the land cover and soil parameters which represent debris-covered glaciers, and these values are added to the VEGPARM.TBL and SOILPARM.TBL tables.
2.2 | Creating the debris-covered glacier category 2.2.1 | Defining the debris cover mask Areas exceeding a certain brightness threshold associated with clean-ice are automatically identified using a cloud free satellite image from Landsat 8, from April 2017. This creates an approximation for the clean-ice area of the glaciers, which is superimposed onto (and cropped by) the glacier area depicted in the RGI. Thus the glacier area is split into a clean-ice area, with the remaining part being debris-covered area. These two new glacier outlines replace the original glacier mask in the 24 USGS land cover categories. A few pixels were manually adjusted to match the surface type, where it was known. Due to the relatively coarse resolution of the model, compared to the width of the glaciers, we did not adjust fine details in the debris-cover mask as these would not be visible at the 1 km resolution. In addition, differences in the glacier cover between recent years would not be apparent at this resolution. See the supporting information section 2 ( Figure S1) for more details.

| Adding the debris cover properties
For the purposes of this study, the parameters for debris cover in the VEGPARM.TBL table are chosen assuming that the debris cover is entirely barren of vegetation, and thus the parameters relating to vegetation are chosen to match those for snow and ice (as there are no other categories that assume the land is entirely barren). The other values used in VEGPARM.TBL have been altered based on measurements reported in previous literature on debris-covered glaciers. A summary of these values, and the reasoning for them, is given in Table 1. The values for SOILPARM.TBL have also been chosen based on previous literature where possible, however few soil properties have been measured over debris-covered glaciers. Those that have are also given in Table 1. Full justification for the parameter choices are given in supporting information section 3. The full VEGPARM.TBL and SOILPARM. TBL can be found at https://github.com/Empott/debris_ cover_in_WRF. Note that it is not possible to compare these values directly to those over the clean-ice glacier in the WRF model, due to the separate glacier scheme used by the Noah LSM scheme. However, the albedo is likely to be substantially lower, and the roughness length higher, in the new debris-cover category, compared to the clean-ice glacier.

| Comparison to observations
Both of the model runs show a cold bias in the daytime average air temperature at 2 m compared to the observations at Changri Nup weather station, but the debriscover run shows an improvement over the clean-ice run, and the two runs are significantly different at 95 %, both in the day and at night. There is a well-defined maximum daily air temperature of just over 7 C in the observations, and a less clear maximum daily temperature of approximately 4 C in the debris-cover run, and under 2 C in the clean-ice run (Figure 2a). The debris-cover run shows very little bias at night, unlike the clean-ice run which is approximately 2 C too cold. While the amplitude of the diurnal cycle appears to have been reduced in the debris-cover run compared to the cleanice run, this is likely to be due to this unrealistic cold bias at night in the clean-ice run at Changri Nup.
The WRF model fails to capture the diurnal cycle in wind speed at 10 m (Figure 2b), although the debriscover run has the lowest RMSE when compared with observations (0.89 ms −1 , compared to 1.01 ms −1 in the clean-ice run). However, the low wind speed and the strong dependence of the wind on very local topography makes comparisons between the model and observations difficult. In addition, the wind comes from a range of directions which the model does not perfectly capture, which may be a partial cause of the inaccurate model wind speed (not shown).
At night the relative humidity at 2 m is relatively well represented in both model runs (Figure 2c). During the day, both model runs are too humid, and the peak minimum daytime relative humidity is approximately 3 hr early. However, the debris-cover run shows a statistically significant improvement over the clean-ice run, with an RMSE of 13.07 %, compared to 20.13 % in the clean-ice run. The model humid bias may also be partially related to the differences between the model and observed wind speed and direction. The debris-cover run also substantially improves the representation of surface incoming and outgoing longwave and shortwave radiation during the day, compared to the clean-ice run, at Changri Nup AWS (see supporting information section 4 [ Figure S2]).
No attempt was made to tune the parameters in the debris-cover category to the observations at Changri Nup, to avoid overcompensating for underlying biases in the WRF model, unrelated to the surface type. However, a number of sensitivity tests were performed to examine the sensitivity of the WRF model to variations in albedo, emissivity and surface roughness over the debris-covered area. Varying each of these parameters in the debriscover run made only a small difference to the model output, compared to the difference between the debris-cover run and clean-ice run (see supporting information section 5 [ Table S2]). As such, it is likely that many of the improvements described above over the debris-covered glacier result from removing the flag which indicates to the Noah LSM that the land cover is a glacier, and hence removing the constraint that the glacier surface cannot be above 0 C, for example. In addition, there are many small snowfall events in the debris-cover (and clean-ice) run, so the debris at Changri Nup is covered with a thin layer of snow for most of the model run. This is likely to diminish the effects of altering albedo, emissivity and surface roughness.

| Effects of adding debris cover to the WRF model
The output from the debris-cover run and the clean-ice run are further compared across all the glacierized regions of the Dudh Koshi River Basin. Note that Figures  3-7 only show areas where the debris-cover run is significantly different (at 95 %) from the clean-ice run, as measured by a t test and a bootstrap test (see Potter et al. (2018)). The debris-cover run substantially increases the surface temperature compared to the clean-ice run over the debris-covered areas of the glacier (Figure 3). The largest increases occur during the day, where the surface temperature is increased by over 18 C in some areas ( Figure  3a). On average, the surface temperature increases from −0.81 C in the clean-ice run to 10.63 C in the debriscover run, over the debris-covered areas of the glaciers during the day (values averaged over all the daytime hours in July, and all the pixels where there is debris F I G U R E 2 The average diurnal air temperature at 2 m (a), wind speed at 10 m (b) and relative humidity at 2 m (c) for July 2013 at the Changri Nup AWS, in local time (LT). The observations are shown (red solid line), along with the two model runs: the original clean-ice category (clean-ice; green dotted line), and the new debris-cover category (debris; blue dashed line). Note that data are missing from the observations for each morning, from about 02:00 to 07:00 LT. Shading indicates one standard deviation from the mean. The root mean square error (RMSE) values between each model run and the observations are also shown cover in the debris-cover run, for the two runs). There is also an increase at night (from an average of −7.53 C in the clean-ice run to 0.23 C in the debris-cover run, over the areas that are debris covered in the debris-cover run) suggesting that the debris cover retains heat past sunset (Figure 3b), and that overall the amplitude of the diurnal F I G U R E 3 The significant differences in surface temperature between the debris-cover run and the clean-ice run (i.e., debris-cover run minus clean-ice run). The grey outline shows the edge of the glacier in the model, and the stippling shows the clean-ice areas of the glaciers. The average daytime (06:00-18:00 LT) values over July 2013 are shown in (a), average nighttime (19:00-05:00 LT) values are shown in (b). Note that the domain has been cropped to show only the glaciers and their immediate surrounding areas, where almost all the significant differences occur F I G U R E 4 As for Figure 3, but for the significant changes to wind velocity (green vectors). For clarity, every second vector in the south-north direction is shown. The underlying contours show the model topography cycle is increased, compared to clean-ice glaciers. Cleanice glaciers can only reach a maximum temperature of around 0 C in the Noah LSM scheme, so this increase is to be expected.
In the clean-ice run, the average wind at 10 m is upvalley during the day, with the strongest winds at the bottom of the valley, becoming weaker over the higher, glacierized regions (not shown). At night, there is a similar pattern with much weaker up-valley winds. The increase in surface temperature over the debris-covered areas of the glacier in the debris-cover run is consistent with a temperature gradient between the warmer surface and the cooler air at the same elevation away from the surface, and therefore a horizontal pressure difference.
F I G U R E 5 As for Figure 3, but for the air temperature at 2 m F I G U R E 6 As for Figure 3, but for the column integrated water vapour Over clean ice glaciers, this temperature gradient is likely to be much weaker, or the other way around. The temperature gradient in the debris-cover run is likely to lead to the simulated increase in upslope winds over the debris-covered ice seen in Figure 4. The increase in upslope winds is largest during the day (Figure 4a), with an average daytime increase in wind speed of 1.8 ms −1 in some locations.
The increase in surface temperature is also seen in the air temperature at 2 m ( Figure 5), however on average the temperature at 2 m increases more at night than during the day (2.89 C and 1.86 C over the debris-covered areas of the glaciers, respectively) with the addition of the debris cover, which is in contrast to the changes in surface temperature. This is also seen in the diurnal cycle at Changri Nup (Figure 2a). This may be due to the increase in near-surface wind during the day causing more mixing than at night, transferring heat away from the near surface.
The increase in up-valley winds is consistent with a transport of water vapour further up the valley. There is an increase in column-integrated water vapour (calculated by multiplying the dry air density by the vapour mixing ratio and integrating over the column [Orr et al., 2017;Cossu and Hocke, 2014]) over parts of both the debris-covered and clean ice glaciers, in the debris-cover run compared to the clean-ice run during the day ( Figure  6a), particularly evident at the boundary between the debris-cover and the clean-ice glaciers, approximately where there is a large increase in daytime winds.
There is also an increase in column-integrated total (liquid + ice) water path (calculated as for vapour, then by summing cloud droplets, rain particles, ice crystals, snow particles and graupel particles) over some of the glacierized areas during the day, in the debris-cover run compared to the clean-ice run (Figure 7a). The average increase is small, but in some places over the clean ice areas the increase is over 0.16 kg m −2 , which is over 50% of the clean-ice daytime average at that location ( Figure  7a). This is likely due to both the increase in water vapour, and the convergence of the horizontal wind occurring between the debris-covered and clean-ice areas of the glaciers during the day. This leads to an increase in vertical wind (not shown), and therefore cloud particle formation. It is not clear what causes the very small increase in water vapour and total water path over the unglacierized areas of the valley. While these may represent a real phenomenon, they may also be due to the sensitivity of the WRF model to small changes.

| DISCUSSION AND CONCLUSIONS
This study presents a new method for adding a representation of thick debris-covered glaciers into the WRF model with the Noah LSM land surface scheme. A model run with the new debris-covered glacier land cover and soil category (debris-cover run) is compared to a model run with the default clean-ice glacier land cover category (clean-ice run), for the topographically-complex Dudh Koshi River Basin, for July 2013. The debris-cover run is an improvement over the clean-ice run, when outputs are compared to observations at Changri Nup of air temperature at 2 m. Although there are biases in the WRF model in relative humidity at 2 m and wind speed at 10m, the debris-cover run does show an improvement, which we would expect to remain if the underlying WRF model biases were not present. To fully test the new debris-cover category, the model should be tested over a range of different areas, and over a longer time period (in particularly for all seasons), which is beyond the scope of this study. It would be useful to compare the simple method described here directly with the more computationally-expensive coupled model proposed by , to determine whether the coupled model reduces the biases seen here compared to observations, or whether these come from the WRF model itself. However, the results presented here suggest that the addition of the debris-cover category improves the ability of the WRF model to accurately represent the weather and climate over glacierized valleys in the Himalaya. As such, a simple uncoupled model may be adequate to represent the effects of debris-cover on the near-surface meteorology when only the atmospheric data are needed, as opposed to information about sub-surface glacier melting.
Output from the two model runs are compared more widely over the Dudh Koshi River Basin. The addition of the debris cover has an effect not just over the debris-cover glacier areas, but also over some of the clean-ice areas (where the land cover category is unchanged). Surface temperature is increased over the debris-cover areas in the debris-cover run compared to the clean-ice run. These findings agree with those of , who compared debris-covered and clean-ice versions of the WRF model as part of a coupled atmosphere-glacier model over a river basin in the Karakoram region. Both the present study and that of  find an increase in surface temperature of around 20 C in some locations, between clean-ice and debris-cover runs. An increase in surface temperature was also found by Reid et al. (2012), who examined the effects of adding debris-cover to a glacier model over a glacier in Switzerland.
In this study, the increase in surface temperature caused by adding debris cover to the model leads to an increase in the upslope wind speed at 10 m over the debris-covered areas. This is in contrast to the findings of , who found that the increase in roughness length over the debris-cover slowed the wind speed at 10 m. This difference in findings could be due to the different representations of debris-cover, or to the differing magnitudes of the opposing effects of increased temperatures and increased surface roughness. It could also be due to differing initial climatic conditions, as  found sustained katabatic winds in the clean-ice run, which are not present in this study. As such, debris-cover may affect the wind speed over glaciers differently in different regions, particularly over glaciers with different topographies or climatic regimes. The increase in upslope wind has consequent effects on other meteorological variables, by driving water vapour further up the valley and increasing the column-integrated total water path over the clean-ice areas to the north of the valley. While the effects discussed here are representative of the monsoon season in the region, the effects of adding debris cover are likely to differ in the dry seasons.
This study demonstrates the importance of including debris-covered glaciers in high-resolution atmospheric models over debris-covered glacier areas, and presents a simple method for doing so. The effects of adding debris cover to regional atmospheric models should be investigated further, over a range of different glaciers and valleys in the Hindu Kush Karakoram Himalaya region.