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.