Introduction

Among the most impacted regions around the world by recent climate change is the central Himalaya. In recent decades, this region has experienced rapidly increasing temperature1,2,3 and decreasing monsoon precipitation.4 These meteorological trends are likely major factors in the rapid retreat of glaciers in the region.5,6,7 These glaciers are depended upon as reservoirs of moisture, which provide water for agriculture, hydropower, and household use as they melt annually. Melting of snow and ice provides a third of annual discharge in the central Himalaya.8 More broadly, the central Himalaya is the source of some of the world’s great rivers, including the Ganges and Brahmaputra, so that precipitation in the region is critical to water resources in the Indian subcontinent.

Reanalyses show an upper-tropospheric circulation trend over Asia in recent decades with an increasing influence of westerlies at low latitudes in recent decades (Fig. 1a), as found by Mölg et al.9 This westerly trend is linked to a cyclonic trend over West Asia and and an anti-cyclonic trend over East Asia. The center of the anti-cyclonic trend approximately corresponds to one of the centers of action of the major teleconnection pattern over Eurasia in the boreal summer known as the “Silk Road Pattern,” with the late 1990s being identified as a shifting point.10 This regional trend represents a disturbance to the upper-tropospheric circulation over Asia during the Indian summer monsoon since 1979.

Fig. 1: Locations of the WRF domains and significant daily averaged trends in JJA.
figure 1

a Statistically significant CFSR trends in JJA from 1979 to 2014 of the indicated variables (wind vectors are plotted where either the u or v trend is statistically significant) and showing the locations of the two nested WRF domains; b, c the area covered by the inner WRF domain showing significant trends from the WRF downscaling of the indicated variables in JJA during the same period. All panels also show the 3-km elevation contour from the relevant dataset indicating the location of the Tibetan Plateau and Himalaya. The rectangle in b, c is the study region, which is plotted in more detail in d, showing elevation (grays, lighter shades indicating higher elevation), glacier coverage (blue), political boundaries, the location of Everest, and the grid spacing of the inner domain (bottom right of panel).

To understand the impact of the regional trend on mesoscale circulations in the Himalaya, a 36-year dynamical downscaling was performed11 from 1979 to 2015 at 6.7-km grid spacing over a region covering the Tibetan Plateau, the Himalaya, and the Karakoram (Fig. 1—domain boundaries shown in panel a and inner domain plotted in panels b, c). The downscaling employed the Weather Research and Forecasting (WRF) model,12 with Climate Forecast System Reanalyses (CFSR13) for initial and boundary conditions. The downscaling was performed to investigate the contrasting behavior between the central Himalaya and Karakoram glaciers,11 where the latter are not uniformly retreating.14,15,16,17 The downscaling shows the central Himalaya has experienced surface warming and drying since 1979 (shown for summer in Fig. 1b, c).

This dynamical downscaling was previously used to calculate trends from seasonally averaged variables.11 However, the Indian summer monsoon, during which almost all of the annual precipitation over the central Himalaya falls,18 is associated with a pronounced diurnal cycle over the Himalaya.19,20,21,22,23 Daytime heating of the slopes generates a horizontal potential-temperature gradient, resulting in an upslope acceleration of the winds; at night, radiative cooling of the slopes creates an opposite-signed temperature gradient, resulting in a downslope acceleration of the winds.24 This diurnal cycle maximizes precipitation over the peaks in the daytime and in valleys at night.19,23 To highlight this local circulation, the analysis in the current study is restricted to a small subregion spanning 86–92°E and 26–29°N (Fig. 1d, location shown in Fig. 1b, c), corresponding to the glaciated areas of the central/eastern Himalaya, just during summer months (JJA). This subregion was chosen for analysis because of the co-location of significant drying (Fig. 1b) and warming (Fig. 1c) trends over the Himalaya during JJA.

This study examines trends in the diurnal cycle in this subregion in JJA. The diurnal cycles of precipitation, temperature, and winds using this WRF configuration over the Himalaya in JJA have been evaluated,23 as have the daily averaged trends (1979–2014) identified by the downscaling.11 In both studies, the level of agreement with observations was sufficient to justify an analysis of the trends in the diurnal cycle using the same downscaling. Here, we hypothesize that the warming and drying over the central Himalaya demonstrated by the downscaling (Fig. 1b, c) are influenced by a complex response of the diurnal cycle of the mountain–valley circulation to the regional summer trend.

Results

The diurnal distribution of precipitation and temperature trends

The climatological diurnal cycle over the Himalaya during the summer monsoon is characterized by upslope (southerly) winds in the afternoon, maximizing precipitation over peaks (Fig. 2, left panels). At night, downslope (northerly) winds generate weaker maxima of precipitation over the foothills (Fig. 2, right panels). However, there is still nocturnal precipitation at high elevations, albeit of lower magnitude than during the afternoon and more evenly distributed over the study region (cf. Fig. 2c, d). This diurnal cycle over the Himalayan peaks and slopes is consistent with that described in the literature based on observations19 and has been described in more detail for a single year of this WRF downscaling.23 The full diurnal cycle of precipitation transitioning between the daytime and nocturnal distributions is shown in Supplementary Fig. 1.

Fig. 2: Climatological values from WRF in JJA over the period 1979–2014 at 1200–1500 (left) and 2100–0000 (right).
figure 2

Top panels show 10-m meridional wind speed (in general, positive values indicate upslope wind and negative values indicate downslope wind) and bottom panels show 3-h precipitation accumulation. All panels also show elevation (contours every 1 km). The given times are UTC + 6, roughly local time.

Regarding trends, the warming at high elevations along the central Himalaya (Fig. 1c) is predominantly during daytime hours, and applies minimally at night (Fig. 3, top row). By contrast, the negative precipitation trend along the central Himalayan arc (Fig. 1b) is predominantly nocturnal (Fig. 3, second row). There are also significant trends in the 10-m winds: southerly (upslope) by day and northerly (downslope) by night (Fig. 3, top row). The downslope trend is by far the more remarkable, being of greater magnitude and covering a greater area. A nocturnal cooling trend is observed at low elevations (Fig. 3b), likely associated with the downslope winds, bringing air from high elevations. The current study is concerned with trends over the mountains and so this cooling trend is not further investigated.

Fig. 3: WRF trends of temperature, winds, and precipitation at specific times of day in JJA.
figure 3

Top two rows: Statistically significant WRF trends in JJA from 1979 to 2014 at 1200–1500 (left) and 2100–0000 (right). Top row: 2-m temperature (colors) and 10-m winds (vectors, plotted where either the u or v trend is statistically significant). Second row: 3-h precipitation accumulation. T and v trends are calculated by averaging between the two given time slices to enable direct comparison with precipitation over the same 3-h period. All panels also show elevation (grays, lighter shades indicating higher elevation). The given times are UTC + 6, roughly local time. Bottom row: The mean diurnal cycle of e 2-m temperature, f 10-m v, and g 3-h precipitation up to the given time, each averaged over the first and last 10 years of the downscaling. Temperature and precipitation are averaged over all grid points >3000-m elevation in the area shown in top panels to highlight the trends at high elevation; v is averaged over all grid points 3000–4000 m, to highlight upslope and downslope winds. Blue dots indicate a statistically significant difference between 1979–1988 and 2005–2014 at the given time of day.

Averaging over the area and just at high elevations, we show that 2-m temperature in fact exhibits warming at all times of day (Fig. 3e), but that the nocturnal trends are not statistically significant. By contrast, changes to 10-m winds and precipitation are non-uniform in the diurnal cycle: winds have become more southerly (upslope) by day and more northerly (downslope) by night (Fig. 3f), while the majority of precipitation decrease is from 1800–0300 local time and statistically significant mostly just in the evening (Fig. 3g). These trends in the diurnal cycle are somewhat simplified by averaging over a large number of grid points, but are representative of the trends in the diurnal cycles of temperature and winds (Supplementary Fig. 2) and precipitation (Supplementary Fig. 3) at individual grid points. The nocturnal trends of enhanced downslope winds and decreased precipitation are particularly remarkable in terms of the fractional change (quantified in Supplementary Fig. 4), given that climatological daytime winds and precipitation are both of greater magnitude than at night.

The remainder of this study will investigate the changing diurnal cycle of the region, specifically: (1) the daytime warming trend at high elevations and upslope trend in near-surface winds; and (2) the nocturnal negative precipitation trend at high elevations and downslope trend in near-surface winds.

The amplification of diurnal upslope and downslope winds

The change in the daytime mountain circulation is shown relative to the climatology in a cross-section at 90°E in Fig. 4 (left panels). Note that this figure shows acceleration of the winds, [\(\frac{\partial v}{\partial t}\), \(\frac{\partial w}{\partial t}\)], as opposed to the winds themselves, [v, w]. This is because the horizontal temperature gradient between the slope and the adjacent atmosphere generates an upslope acceleration of winds,24 which does not necessarily imply an upslope wind (i.e., winds may be downslope prior to the onset of the temperature differential and just become weaker downslope winds). Therefore, the potential-temperature contours and acceleration of winds best illustrate the dynamic response to the thermal forcing.

Fig. 4: A cross-section at 90°E at 0900–1200 (left) and 1800–2100 (right), showing climatology and trends from WRF of the mountain–valley circulation in JJA from 1979 to 2014.
figure 4

Top row: Climatological potential temperature, θ (contours every 1 K), and wind vectors with x and y components of \(\frac{\partial v}{\partial t}\) and \(\frac{\partial w}{\partial t}\), respectively, representing acceleration of cross-section winds. Second row, \(\frac{\partial \theta }{\partial y}\), illustrating the temperature differential between the slopes and free atmosphere. Third row: Statistically significant trends of θ (colors) and of [\(\frac{\partial v}{\partial t}\), \(\frac{\partial w}{\partial t}\)] (vectors, plotted where either the x or y component is significant), illustrating trends in the acceleration of cross-section winds. Bottom row: The horizontal gradient in the y-direction of the θ trend, illustrating the differential warming between the slopes and free atmosphere. In the panels with wind vectors, the vertical component is re-scaled according to the relative scales of the y- and z-axis. In the second and bottom rows, the large values of alternating sign over the peaks arise from the undulations in the θ contours over jagged terrain and are not the focus.

In the climatology, the warming of the slopes in the daytime, relative to the adjacent free troposphere at the same altitude accelerates the winds toward the mountain (Fig. 4a), as in conventional depictions of anabatic winds24 with the warming of the slopes more clearly illustrated by the horizontal potential-temperature gradient, \(\frac{\partial \theta }{\partial y}\) (Fig. 4c). The trends clearly show an enhancement of this diurnal cycle: there is greater warming over the slope than the adjacent free troposphere, hence there is a trend in acceleration of winds up the mountain between about 800 and 500 hPa (Fig. 4e), with the greater warming of the slopes more clearly shown by the horizontal gradient of the θ trend (Fig. 4g).

By contrast, in the nocturnal climatology, the slopes cool faster than the adjacent free troposphere at the same altitude, accelerating winds down the mountain (Fig. 4b), as in conventional depictions of katabatic winds,24 with the cooling of the slopes more clearly shown by \(\frac{\partial \theta }{\partial y}\) (Fig. 4d). As in the daytime, the trends show an enhancement of the climatological circulation: there is less of a warming trend over the slope than the adjacent free troposphere, causing a trend in acceleration of winds down the mountain between about 800 and 600 hPa (Fig. 4f), with the lesser warming of the slopes more clearly shown by the horizontal gradient of the θ trend (Fig. 4h).

The period after sunset is unique in having a downslope trend in acceleration of the winds, that is, around sunset the transition is made from an upslope trend in the daytime to a downslope trend at night (Supplementary Fig. 2). Similarly, the period after sunrise is unique in having an upslope trend in acceleration. Because, climatologically, the acceleration of the winds arises from a horizontal temperature gradient, it is evident that the mountain–valley circulation has been enhanced by an increasing horizontal temperature differential near the slopes, after both sunrise and sunset.

Factors associated with daytime warming of the Himalaya

Why have mountain slopes warmed at a greater rate than the surrounding free troposphere in the daytime and how does this relate to enhanced anabatic winds (Fig. 4e)? The warming at high elevations is explained most simply by global warming, which feeds into the WRF simulations from the reanalyses via spectral nudging (see ref. 11 for more details of the downscaling methodology). Temperature, among other variables, was nudged in the WRF simulations from the reanalyses at all vertical levels. In the WRF simulations, the upper troposphere in the study region has warmed during the downscaling period (Fig. 4e, f). Viewed as horizontal maps in the upper and mid-troposphere in both CFSR and WRF over the full WRF inner domain, it is evident that the distribution of significant warming trends in WRF is derived from that in the reanalyses, albeit with some mesoscale structure from internal variability in the model (Supplementary Fig. 5). Because the warming in our study region is just from about 650 hPa upward (Fig. 4e, f), only high elevations experience surface warming, explaining the distribution of afternoon 2-m warming trends in Fig. 3a.

According to the WRF simulations, there is no radiative contribution to the warming, given that there are near-neutral trends in both shortwave and longwave net (downward minus upward) radiation over the mountains in the early afternoon (Fig. 5a, b). There are just scattered negative trends of net shortwave and positive trends of net longwave, which are not comparable to the positive trends in 2-m temperature that are widespread at high elevations at the same time of day (cf. Fig. 3a). There is in fact a negative trend in downward shortwave radiation (i.e., not accounting for the reflected radiation) over the central Himalayan peaks (Supplementary Fig. 9b). This negative trend in shortwave radiation is related to the enhanced anabatic winds, described in the previous section.

Climatologically, upslope winds induce convergence over peaks and divergence over valleys; thus, the upslope trend is convergent over peaks and divergent in valleys (Fig. 5c; and shown in a zoomed-in plot that illuminates individual peaks and valleys in Fig. 5d), and accordingly, there are some increasing trends in afternoon cloudiness over the peaks and decreasing trends in valleys (Fig. 5c, d; full diurnal cycle of cloudiness shown in Supplementary Figs. 6–8). Hence, incoming shortwave is reduced over the peaks.

Fig. 5: Statistically significant afternoon (1200–1500) WRF trends in JJA from 1979 to 2014.
figure 5

a Net (downward minus upward) shortwave radiation at the surface; b net (downward minus upward) longwave radiation at the surface; c cloud fraction at 460 hPa (chosen because it is just above the high peaks in the WRF domain, so shows trends in low-level cloud over the central Himalaya) and 10-m winds; d as in c but in a zoomed-in plot, with the location indicated by the box in c. The actual radiation received at the surface is also influenced by slope, but since this is constant over time it does not affect the sign of trends. In c, d, vectors are plotted where either the u or v trend is significant. All panels show elevation (grays, lighter shades indicating higher elevation).

The near-neutral trends in net shortwave (Fig. 5a), despite these decreases in downward shortwave illustrate the effect of the Himalaya becoming less reflective, that is, decreased albedo (Supplementary Fig. 9a). Decreasing reflectivity of the peaks is largely a function of decreasing snow cover, which has been confirmed from observations.25,26 The decreasing snow cover is a result of both warming (Fig. 3) and a negative trend in pre-monsoon snowfall (not shown, but captured by the downscaling over Himalayan peaks). These negative influences on snow cover are well documented by observations.26 Therefore, despite enhanced afternoon cloud cover resulting from enhanced upslope winds, the reduced albedo of the Himalaya’s high elevations imply that a greater fraction of the incoming radiation is absorbed by the surface, so that the downscaling exhibits near-neutral net shortwave trends.

Factors associated with nocturnal drying of the Himalaya

Why have the slopes warmed at a slower rate than the surrounding free troposphere at nighttime, resulting in enhanced downslope winds (Fig. 4f)? The downscaling exhibits significant negative trends in both upper-tropospheric (Fig. 6a) and mid-tropospheric cloud fraction (Fig. 6b) in the evening. These trends in cloud cover are widespread through the day (Supplementary Figs. 6 and 7; shown as cross-sections at 90°E in Supplementary Fig. 8). However, at sunset the reduced cloudiness has the unique effect of increasing top-of-atmosphere outgoing longwave radiation (OLR, Fig. 6c). Climatologically, as the sun sets, longwave radiative cooling of mountain slopes occurs faster than cooling of the adjacent free troposphere,24 resulting in the downslope acceleration depicted in Fig. 4b. The negative OLR trend illustrates that this effect has been amplified in recent decades, explaining why the slopes exhibit less of a warming trend than the adjacent free troposphere at sunset (Fig. 4f, h). In particular, clearer skies have increased the radiative cooling of the slopes. The co-location of the OLR and downslope trends in the early evening (Fig. 6c, d) illustrates the response of the katabatic winds to this radiative cooling. Enhanced downslope winds and reduced cloudiness also explain the negative nocturnal precipitation trend (Fig. 3d). This argument does not apply in the daytime because the enhanced anabatic winds have enhanced convergence over peaks.

Fig. 6: Statistically significant nocturnal (2100–0000) WRF trends in JJA from 1979 to 2014.
figure 6

a Cloud fraction at 300 hPa, b cloud fraction at 460 hPa (chosen because it is just above the high peaks in the WRF domain, so shows trends in low-level cloud over the central Himalaya); c outgoing longwave radiation at the top of the atmosphere and 10-m winds; d as in c but in a zoomed-in plot, with the location indicated by the box in c. In c, d, vectors are plotted where either the u or v trend is significant. All panels show elevation (grays, lighter shades indicating higher elevation).

It is interesting that this negative precipitation trend is restricted to above about 3-km elevation (Fig. 3d), which is likely explainable by the downslope trend. Zooming further into the study region, the trend of downslope winds is convergent in valleys and divergent over peaks (Fig. 6d). At low elevations, the enhancement of surface convergence due to the downslope winds prevents drying. This contrast in precipitation trend between low and high elevations is also evident in observations that show no summer precipitation trend averaged across stations.26 Almost all of these stations are at low elevations (<3000 m). The same observational dataset also shows a contrast between north (the high Himalaya) and south (the foothills) Nepal in summer, with no significant trends to the south, but with some significant drying trends to the north. The analysis afforded by our high-resolution multi-decadal dataset is in a unique position to explain trends with both mesoscale and synoptic-scale influences: summarily, an enhancement of nocturnal downslope winds has reduced precipitation at high, but not low, elevations.

Discussion

Because the WRF simulations are forced by CFSR, we subsequently interpret these reanalysis trends and discuss their reliability. The location of the positive trend in 200-hPa geopotential height (Fig. 1a) is one of the centers of action of the “Silk Road Pattern,” whose shift in the late 1990s explains most of the shift in precipitation and surface temperature since that time over large parts of Eurasia, including the Tibetan Plateau, according to reanalyses.10 This shift may also be identified as a positive shift in the “Karakoram Vortex,”27,28 approximately corresponding to the cyclonic part of the trend, over southwestern Asia (Fig. 1a). A cyclonic anomaly in that location indicates the positive phase of the Karakoram Zonal Index, which in summer months has become stronger and more frequent in recent years,27 consistent with our results. It is beyond the scope of this study to explain the source of this large-scale trend, but, as discussed in the introduction, these trends may be monsoon related, or from an enhancement of westerlies, and could result from warming of the Tibetan Plateau, all of which are likely inter-connected.

Along with the Interim European Center for Medium-Range Weather Forecasts (ECMWF) ReAnalysis (ERA-Interim), CFSR is one of the best performing reanalyses over the Tibetan Plateau, evaluated by 3-D observations.29 ERA-Interim30 shows a similar regional trend in JJA, as do the Modern-Era Retrospective analysis for Research and Applications 2 (MERRA 231) and the Japanese 55-year Reanalysis (JRA-5532) (Supplementary Fig. 10). There are variations between reanalyses in the magnitude and area over which the anti-cyclonic trend is significant, but there is agreement that the anti-cyclonic trend is significant over most of High Mountain Asia, extending southward to the Himalaya. These comparisons indicate that the reanalysis trends downscaled by WRF are not unique to CFSR.

Regarding the reliability of the WRF simulations themselves, the amplification of the diurnal pattern by the regional trend has occurred by contrasting mechanisms between day and night, both of which involve the interaction of dynamics, clouds, and radiation. The dynamics associated with the trends should be reliable, provided that sufficient resolution is used. Although the grid spacing employed in this study (6.7 km) is more coarse than the <4-km or even ~1-km usually considered necessary for resolving orographic flow patterns,33 greater resolution was not computationally feasible for 36 years of simulation. Importantly, no fundamentally different patterns in the diurnal cycle over the Himalaya are simulated with WRF with a 2.2-km over a 6.7-km domain, despite some differences in the magnitude of winds and precipitation.23 It is likely that some events during the downscaling would have been simulated differently had convective parameterization been employed in the simulations (see Methods), which could have impacted some of the trends. However, the central conclusions of the study appeal to dynamical, rather than convective, arguments, and should not be sensitive to whether convective parameterization was employed.

To simulate clouds and radiation, this downscaling employed the Thompson microphysics34 and Rapid Radiative Transfer Model for GCMs (RRTMG) scheme,35 where the microphysics interacts with the dynamics, and the radiation depends on the microphysics. The RRTMG scheme uses cloud water, cloud ice, and snow concentrations from the microphysics when calculating radiative fluxes, which performs well when coupled with Thompson microphysics.36 The interannual variability of cloudiness in JJA over the central Himalaya was compared between the WRF downscaling and the Moderate Resolution Imaging Spectroradiometer Cloud Product, with the level of agreement sufficient to justify an analysis of the WRF trends.11 It is likely that the daytime cloud trends are more difficult for the model to represent accurately than the nighttime trends, given the greater convective influence in the daytime and the relatively coarse grid spacing. Therefore, the daytime cloud cover may not have increased as shown in Supplementary Figs. 7 and 8, although this is consistent with the enhanced anabatic upslope trend in which we have more confidence from the model, given that the anabatic winds are dynamically, rather than convectively, forced. However, the daytime cloud trends are not central to our conclusions, whereas we argued that the nocturnal decreasing trends in cloud cover shown by the model have led to the enhanced downslope winds and reduced precipitation. Given the lesser dependence of the nocturnal cloud on convective processes, these simulated trends should be more reliable.

Errors may also have arisen from the initialization of the surface properties by the reanalyses. This is relevant to the negative trend in albedo shown in Supplementary Fig. 9a, which we argued has caused more incoming shortwave radiation to be absorbed. Because our simulations run from March of one year to April of the following year (see Methods), the simulated winter snowpack is not effective during the monsoon. Hence, all surface properties were initialized from the reanalyses, but with 3 months of spin up prior to the simulation of the monsoon each year. The surface properties are therefore more realistic than those derived from the reanalyses (which, particularly for end-of-winter snow and glacier coverage, are very poorly represented, due to the low resolution), but less than if, for example, a full year of spin up had been employed. We argue that the negative trend in snow cover during the monsoon is in fact probably of greater magnitude than that captured by the downscaling because of a negative trend in winter precipitation captured by the downscaling,11 but not effective in the summer simulation. This effect would lead to an even greater amount of downward shortwave radiation being absorbed, hence a greater daytime warming trend than was simulated by our downscaling (Fig. 3a), hence an even greater trend in the upslope anabatic winds (Figs. 3, 4, and 5).

Note that the above argument is dependent on the balance between incoming shortwave radiation and surface albedo, and the approximate balance between the decreasing incoming shortwave and decreasing surface albedo (i.e., a greater fraction is absorbed) exhibited in this study (Supplementary Fig. 9) may be specific to the model configuration and parameterizations employed. Therefore, we refrain from making firm conclusions on trends in net shortwave radiation, but present the trends captured by the model to illustrate the interplay between these two effects.

The results of this study are highly relevant to glaciers in the central Himalaya, which exhibit some of the most rapid retreat rates on earth in recent decades.5,6 Although the causes of this retreat are probably not purely meteorological, the warming and drying identified by this WRF downscaling over the central Himalaya in both winter and summer11 are both consistent with and likely factors in the retreat of the glaciers. Given that glacier melt occurs primarily in summer months, warming in the summer, particularly in the daytime as captured by our analysis, indicates that a greater fraction of glacier mass has melted over time. Although winter and spring precipitation trends would also be impactful to glaciers due to the larger fraction of solid precipitation, the most robust precipitation decreases over the central Himalaya in the downscaling are shown in the summer,11 which is when that region receives almost all of its precipitation.18 Because, even in summer, some of that precipitation is solid, summer trends are also relevant for glacier mass balance, and more broadly for water resources downstream. Therefore, we conclude that the summer trends in both precipitation and temperature have likely played the largest role in affecting water resources in regions influenced by glacier melt. This study has identified that this warming and drying in summer are in the daytime and nocturnal, respectively. Therefore, our results suggest that the amplified diurnal cycle of daytime upslope and nocturnal downslope winds has been a large factor in glacier retreat and general disruption to water resources in recent decades.

Finally, we note that the mesoscale and orographic trends identified by this study highlight profound shortcomings of global and even regional models over the Himalaya, whose twenty-first century projections have been used to drive hydrological models and hence project future glacier mass balance and runoff over this region.8,37,38,39,40 Using this approach, in the representative concentration pathway (RCP) 2.6 scenario, 10.6% of glaciers over the Himalaya and Karakoram are projected to have disappeared by the end of the twenty-first century, whereas in the RCP 8.5 scenario, 27% will have disappeared in the same period.39 Some river catchments are particularly vulnerable, for example, the Langtang catchment in the Nepalese Himalaya, which is projected to lose one third of its ice under the RCP 4.5 scenario by 2100 and two-thirds under the RCP 8.5 scenario, and the Baltoro catchment in the western Himalaya, which is projected to lose a half of its ice under the RCP4.5 scenario and two-thirds under the RCP8.5 scenario.38 As the glaciers retreat, increases in runoff are projected throughout the Himalaya until about 2050, after which runoff begins to decrease.8,38

These hydrological projections are dependent on the temperature and precipitation from GCMs, which are both projected to increase through the twenty-first century throughout the Himalaya.37,38,39,40 The temperature projections exhibit relatively little spread between models, generally within about a quarter of the magnitude of warming.40,41 Temperature projections should therefore be relatively robust, but it is well documented that higher elevations are projected to warm at greater rates,42 and GCMs do not capture the high elevation of the glaciers, which may lead to underestimates of warming across models. Moreover, the ability of GCMs to represent the relationship between snowpack and radiative warming, illustrated to be important for daytime temperature trends by the current study, is questionable. The greatest shortcoming of GCM projections over complex terrain, however, is in precipitation. Based on GCM output used to drive a high-resolution glacio-hydrological model over the Himalaya, the largest contribution of uncertainty in runoff arises from the uncertainty in precipitation projections, with model spread of >100% of the magnitude of the trend for a given forcing scenario.38 Because of the lack of representation of nocturnal precipitation at high elevations in the Himalaya in GCMs, and even with WRF at 10-km grid spacing,23 it is quite possible that twenty-first century precipitation trends over the Himalaya will not even be within this large spread of GCM projections. In particular, GCMs near uniformly project wetting, but there may be drying in parts of the Himalaya if nocturnal precipitation continues to decrease via a strengthening of downslope winds as documented in this study.

Evidently, the uncertainty in projected precipitation and temperature trends implies similar uncertainty in glacier mass balance and runoff over the Himalaya. Although GCM projections are currently the best available data for twenty-first century temperature and precipitation, and their use to drive hydrological models and project water resources is insightful, GCMs fail to represent crucial processes that determine temperature and precipitation trends in the Himalaya, highlighted in the current study. As computational capabilities improve, it will become possible to use high-resolution meteorological projections over the Himalaya and other complex terrain that take into account the orographic flow patterns that respond uniquely to regional climate change. Consequently, a greater understanding will be gained of how water resources will be affected through the twenty-first century.

Methods

WRF downscaling

This study utilizes the 36-year dynamical downscaling with the Advanced Research WRF 3.7.112 from April 1979 to March 2015 over High Mountain Asia.11 Initial and boundary conditions were provided by CFSR,13 and downscaling was performed with nested model grids of 20- and 6.7-km grid spacing (locations shown in Fig. 1a; inner domain plotted in Fig. 1b, c). The inner domain was chosen to cover almost the full Tibetan Plateau, and Himalayan and Karakoram ranges, and the boundaries of the outer domain were chosen following tests of model stability, both in the inner and outer domains. Although the southern boundaries of the two domains are close, tests revealed that no major differences in precipitation features in the inner domain were captured when the outer domain extended farther south. The downscaling consists of 36 separate simulations of 13 months initialized at 00 UTC (Coordinated Universal Time) 1 March and run through 00 UTC 1 April of the following year, with the initial March discarded as spin up. Thus, the JJA period from each year of simulation, analyzed in this study, follows 3 months of spin up, which is relevant for the land surface model (LSM). The Noah multiparameterization LSM was used for the simulations, which more realistically simulates snow and ice processes than other LSMs,43 with all surface properties initialized by the reanalyses. The simulations were performed with no convective parameterization in the inner domain and with parameterizations and other details of the model configuration previously detailed.11 We acknowledge that our chosen grid spacing is in the “gray zone” in which the explicit model dynamics are almost but not fully capable of resolving convective features that are parameterized at coarser scales.44 In some cases, it may be preferable at 6.7-km grid spacing to employ convective parameterization and there may be some convective–orographic systems during the downscaling period that could have been better resolved with convective paramterization. Although 6.7 km is coarse in terms of explicitly resolving convection, no fundamentally different patterns in the diurnal cycle over the Himalaya were simulated with a 2.2-km domain for a few representative days during the monsoon.23 This issue should be considered for further downscaling studies and considered when interpreting this downscaling’s results, particularly in the daytime.

Calculation of diurnal trends

All relevant two- and three-dimensional (3-D) variables were archived at 3-hourly resolution, permitting an analysis of the diurnal cycle. The analysis in this study is performed by averaging a given variable at a given time of day (i.e., one of the eight output times) across JJA of one year at each grid point (i.e., taking the mean of 92 days). During JJA, there is very little variability in solar declination over this subregion (over Mount Everest, sunrise varies between 0715 and 0750 and sunset varies between 2034 and 2112—https://www.timeanddate.com), so that each 3-hourly diurnal time slice may be taken to be effectively the same part of the day throughout summer. At each grid point and for each time of day, this approach generated a 36-element time series, where each element represents the mean over one summer of the given variable at the given time of day. For 3-D variables, output was first interpolated onto pressure levels and then a time series generated for a given variable at a given pressure level and grid point.

These time series were then used to perform a two-tailed Monte Carlo significance test. This test consists of taking the given 36 values and randomly shuffling the values 1000 times, calculating the trend for each replication. If the magnitude of the true trend is greater than the magnitudes of the trends of 95% of replications, then the trend is taken to be significant. Grid points with missing data in the figures are where the given trend was not statistically significant, based on this test. Figures showing climatological values are where the mean was calculated of the 36 elements in the given time series at the given time of day. Fractional trends (Supplementary Fig. 4) were obtained by dividing the significant trends by the climatological values.

For precipitation, trends apply to the accumulation during the given 3-h period. Other variables were averaged over the given 3-h period for direct comparison with the precipitation accumulation at the same time. For example, for the 0000–0300 time slice the precipitation trend applies to the precipitation accumulated during those 3 h, whereas for other variables the same time slice was obtained by interpolating between the instantaneous values at 0000 and 0300 and then calculating the trend based on the interpolated values. Trends were also calculated using various percentile regressions, for example, calculating the diurnal trends based on the wettest 10, 20, 30, or 45 days of each JJA, but using all 92 days of JJA returned the most robust trends.

Certain trends were also calculated for WRF data based on daily averaged fields (Figs 1b, c and Supplementary Fig. 7a). For these figures, the mean value for the given year was calculated from all daily output times in JJA, yielding a 36-element time series upon which trends were calculated the same as for the diurnal trends.

For the temperature and wind line plots in Fig. 3e, f, there was no interpolation performed, that is, the value plotted at a given time of day on the x-axis is the mean across all days in JJA over either 1979–1988 or 2005–2014 at that specific time, averaged over all grid points in the given elevation range. For the precipitation line plot in Fig. 3g, the given value is the accumulation in the 3 h up to the given time on the x-axis, for example, the value plotted at 0000 is the mean accumulation from 2100 to 0000. In these line plots in Fig. 3e–g, significance was established by computing the mean value for each year of downscaling for the given variable and time of day, yielding 10 values for 1979–1988 and 10 values for 2005–2014. The Student’s T-statistic and probability that two sets of 10 values have significantly different means were then computed, assuming that the two sets of values have the same true variance. The markers on the 2005–2014 curves are where there is a probability of <0.05 that 1979–1988 and 2005–2014 have the same true mean.

Reanalysis data

In addition to CFSR, used for initial and boundary conditions of the WRF simulations, the ERA-Interim,30 MERRA 2,31 and JRA-5532 reanalyses were all used to compare to the trends captured by CFSR. To calculate the trends in Fig. 1a and Supplementary Fig. 9, the JJA mean (over all times of day) of the given variable was taken in the given dataset for each available year during the downscaling period (1980–2014 in MERRA 2 and 1979–2014 in the other reanalyses). At each grid point, the same Monte Carlo test was performed as above to establish significance of the trends.