1 Introduction

The Western Himalaya and surrounding region receive a significant fraction of their precipitation during the winter months of December to March (e.g., Singh et al. 1995; Archer and Fowler 2004; Iqbal and IIyas 2013; Zaz et al. 2019) as the subtropical jet moves southwards, bringing extratropical storms known as western disturbances (WDs) to the area. WDs bring rainfall to lower elevations, where it is important for agriculture, and snow to higher elevations, where it is important for replenishing glaciers and deepening the snowpack, both of which are subsequently vital for water security during the spring and summer months. Variability in WD frequency and behaviour can thus affect the ecology, socioeconomics, and population health of the Himalayan region. WDs typically originate over the Mediterranean, Europe, or North Atlantic and travel as perturbations in the subtropical jet until they intensify upon reaching South Asia. WDs are often associated with hydro-meteorological disasters like floods, landslides and avalanches across the Himalaya and surrounding regions. Because WDs bring the majority of winter precipitation to the Western Himalaya (Midhuna et al. 2020), and because they are intimately linked with the behaviour of the subtropical jet (Hunt et al. 2018), there is the opportunity for seasonal predictability through global teleconnections that influence the position and intensity of the jet, such as the North Atlantic Oscillation (NAO) and El Niño Southern Oscillation (ENSO). The NAO is a major mode of winter variability in the Northern Hemisphere, its effect extending from North America to Europe and a large portion of Asia (Hurrell 2003) and accounts for more than 36% of the variance in mean sea level pressure during winter months (Walker and Bliss 1932).

The statistical relationship between the NAO and winter precipitation over the western Himalaya has been explored previously, but authors have disagreed on both the sign and significance of the relationship. Using a 21-year rolling window, Yadav et al. (2009) showed that the NAO is positively correlated with winter precipitation over northwest India, and that this relationship was strongest between 1940 and 1980 (when the mean correlation coefficient was 0.42). Using composite analysis, they showed that positive phases of the NAO during this period resulted in a stronger meridional pressure gradient over Europe and North Africa, resulting in a stronger subtropical jet, and hence, they hypothesised, more intense WDs—a result later corroborated by Attada et al. (2018). They noted a weakening of this relationship over recent decades in favour of a strengthening relationship with ENSO. A similar analysis by Nageswararao et al. (2016a) also found a positive correlation. Midhuna and Dimri (2019) found increased WD activity and colder, wetter, winters occurred during years with a positive AO and linked this to increased strength of westerlies originating over North Atlantic, in turn associated with a stronger upper-level meridional pressure gradient. They suggested that a weaker Hadley cell during these years allowed for greater equatorward incursion of midlatitude air masses. Zaz et al. (2018) found a strong positive correlation (0.68) between the NAO and seasonal mean winter precipitation over Jammu and Kashmir, which they attributed to local meridional displacement of the subtropical westerly jet. Like Nageswararao et al. (2016a), they highlighted an increasing influence of the Arctic Oscillation over recent decades. Filippi et al. (2014) found a positive correlation between the NAO and winter precipitation over the Hindu-Kush-Karakoram in the Western Himalaya but—like Yadav et al. (2009)—demonstrated that the strength of the teleconnection varies on interdecadal timescales and is sensitive to the angle of the jet over the North Atlantic.

Other studies have found a positive correlation through more indirect means. Nageswararao et al. (2018) showed that wheat production (typically grown in the winter months) in northwest India was strongly correlated with NAO and AO activity. Forke et al. (2019) analysed carbonate sediments in the Arabian Sea and found that increased winter runoff was associated with periods of positive ENSO and positive NAO. Greene and Robertson (2017) found evidence of a link between the NAO and decadal variability of winter and spring precipitation in the Upper Indus Basin in CMIP5 models. Hingmire et al. (2019) found that non-foggy days over the Indo-Gangetic Plain—which they associated with the presence of WDs—preferentially occurred during periods of positive NAO, although this result may be complicated by the fact that WDs often are associated with fog (e.g., Sawaisarje et al. 2014; Patil et al. 2020). Syed et al. (2010) used a regional climate model to investigate the effects of the NAO and ENSO on precipitation over North Pakistan, Afghanistan and Tajikistan and found that positive phases of the NAO were associated with anomalously low 500 hPa height over the region and associated moisture intake from Mediterranean, Caspian, and Arabian Seas, leading to excess precipitation.

In contrast, some studies have also asserted a negative correlation between the NAO and winter precipitation over the western Himalaya. Devi et al. (2020) constructed a simple seasonal forecast model to predict winter precipitation over the Western Himalaya and found a strong negative correlation with the NAO (− 0.57). Basu et al. (2017) attributed a weak 2015 winter monsoon to reduced WD activity caused by positive phases of ENSO and the NAO, stating that the positive NAO set up an unfavourable standing wave in the subtropical jet with a blocking anticyclonic anomaly over the subcontinent. Using plankton-based proxies, Munz et al. (2017) showed that intense winter monsoons since 1750 coincided with strong negative phases of the NAO. Similarly, Giosan et al. (2018) used plankton-based proxies to show that strong winter monsoons during the mid-Holocene (4500 to 3000 years ago) were contemporaneous with a more negative NAO and southward migration of the subtropical jet, bringing increased WD activity to South Asia.

Not only do these studies disagree on the sign and significance of the relationship, but they are mostly statistical in nature. Our study proposes to resolve this by establishing causality between the NAO and downstream winter precipitation over the Western Himalaya, through the medium of WDs. We start by showing how the subtropical jet is affected by different NAO phases, then use wavelet analysis to link this relationship to WD behaviour. We then use composite analysis to show how this North Atlantic variability projects onto WD behaviour, moisture flux, and precipitation. Finally, we develop a simple model, linking North Atlantic variability with characteristics of the subtropical jet, and evaluate the ability of the model to forecast precipitation over the Western Himalaya at a lead time of one month. The results of this study may help improve seasonal forecasts for WD-related winter precipitation which will in turn help in developing policies and mitigation strategies related to socio-economic and natural hazards associated with these disturbances.

2 Data and methods

2.1 ERA5

To quantify moisture transport and the synoptic-scale structure of the troposphere over the Western Himalaya, we use data from the ECMWF ERA5 reanalysis (Hersbach et al. 2020). Data are available globally, at hourly resolution from 1950 onwards, on a 0.25° × 0.25° grid. Data are available over 37 pressure levels from 1000 to 0.01 hPa, as well as at selected heights above the surface. Data are assimilated into the forecasting system from a large variety of sources, including satellites, automatic weather stations, and radiosondes. We also use ERA5 as an alternative source of precipitation data for verification. Data can be freely downloaded from https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview.

2.2 APHRODITE

Precipitation data in this study come from the APHRODITE-2 project (Yatagai et al. 2012, 2020) which is based on a dense network of in-situ gauge data, gridded using a kriging method to 0.25° × 0.25° at daily resolution. The dataset runs from 1951 to 2015. Over India, the rainfall amount recorded for a given day is the accumulation from 0300 UTC (0830 LT) the previous day through 0300 UTC of the current day. This dataset was chosen for its longevity and resolution. Given we are not interested in rainfall over the ocean, satellite data are not needed. Despite the relatively sparse gauge network over the Himalaya and Hindu Kush, APHRODITE was shown to outperform other datasets, including reanalyses and satellite products in this area (Baudouin et al. 2020). Data can be freely downloaded from http://aphrodite.st.hirosaki-u.ac.jp/products.html.

2.3 Other precipitation datasets

Accurately measuring precipitation in northern India can be challenging, on account of the complex orography and typically sparse population. We therefore occasionally need to verify results using a variety of precipitation datasets from different sources. These are described briefly below, with the exception of ERA5, which is discussed in Sect. 2.1.

IMDAA (Indian Monsoon Data Assimilation and Analysis; Rani et al. 2021) is a high-resolution reanalysis focused over the India region, available at 0.12° resolution an hourly frequency between 1979 and 2018. Developed jointly by the UK Met Office and the Indian National Centre for Medium-Range Weather Forecasting, it uses the Met Office Unified Model and assimilates a large volume of in-situ, buoy, aircraft, and satellite observations over India. Preliminary analyses suggest it captures interannual variability well but has a wet bias over the Western Himalaya (Sharma et al. 2022). Data can be freely downloaded from https://rds.ncmrwf.gov.in/.

IMD (India Meteorological Department gridded precipitation dataset; Pai et al. 2015) is a gauge-only gridded dataset with 0.25° resolution produced by the IMD, available at daily frequency from 1901 onwards. This dataset provides a valuable alternative source of in-situ precipitation observations but comes with the typical caveats of low gauge density in the Western Himalaya and under-representation of snowfall due to blowing (Dahri et al. 2018). Data can be freely downloaded from https://cdsp.imdpune.gov.in/home_gridded_data.php.

GPM IMERG (Global Precipitation Mission—Integrated Multi-satellitE Retrievals for GPM; Huffman et al. 2015) has global coverage at a half-hourly, 0.1° resolution, starting June 2000 and continuing to the present day. Over the tropics, GPM IMERG primarily ingests retrievals from (for 2000–2014) the now-defunct Tropical Rainfall Measuring Mission (Kummerow et al. 1998; 2000) 13.8 ~ GHz precipitation radar and microwave imager (Kozu et al. 2001) and (for 2014-) the Global Precipitation Measurement (GPM; Hou et al. 2014) Ka/Ku-band dual-frequency precipitation radar. When an overpass is not available, precipitation is estimated by calibrating infrared measurements from geostationary satellites. While GPM IMERG performs well when compared against gauge-based products, performance falls at higher elevations or when quantifying extreme precipitation events (Prakash et al. 2018), although it performs well over the Western Himalaya (Baudouin et al. 2020). Data can be freely downloaded from https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGHH_06/summary?keywords=%22IMERG%20final%22.

2.4 Western disturbance database

For WD track and intensity information, we use the University of Reading ERA5-T42 track database, computed by applying the algorithm developed in Hunt et al. (2018). The mean relative vorticity in the 450–300 hPa is computed, and then spectrally truncated at T42 to remove small-scale features and other short-wavelength noise. Regions of positive vorticity are then identified, and their centroids linked subject to a set of physical constraints. The database is then filtered, removing tracks which do not pass-through north India (20°N–36.5°N, 60°E–80°E), last less than 48 h, or dissipate at a point westward of their genesis. These data can be downloaded from https://gws-access.jasmin.ac.uk/public/incompass/kieran/track_data/WD/ and are described in greater detail in Sharma et al. (2022).

2.5 NAO and ENSO indices

For NAO data, we use the monthly index published by the NWS CPC, available from 1950 onwards. The index is computed by applying a rotated principal component analysis (Barnston and Livezey 1987) to standardised 500 hPa geopotential height anomalies between 20° and 90°N. The NAO has highest variance in the winter months. Throughout, we define NAO+ as months where the NAO index is greater than 1 (giving 160 months in total, 49 in winter), and NAO− as months where the index is less than -1 (giving 142 months in total, 58 in winter). NAO data can be freely downloaded from https://www.cpc.ncep.noaa.gov/products/precip/CWlink/ pna/nao.shtml.

For ENSO, we use the CPC Oceanic Niño Index (ONI), which can be downloaded from https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php. This index is computed as the mean SST anomaly in the box [5°S–5°N, 120–170°W], computed against a 30-year baseline that is updated every 5 years.

2.6 Wavelet spectra

Wavelet power spectra are computed using the PyCWT library (https://pycwt.readthedocs.io/), which is an implementation of Torrence and Compo (1998). For our wavelet analysis, we use a derivative-of-Gaussian continuous wavelet. This only has a real part, and thus computed power identifies the individual peaks and troughs of the underlying data.

2.7 Empirical orthogonal functions

To explore the variability of the jets over the North Atlantic and their relationship with the NAO, we use empirical orthogonal function (EOF) analysis, a method by which a timeseries is decomposed into orthogonal (i.e. uncorrelated with each other) basis functions that seek to maximise explained variance. In this study, we use the eofs package in Python (Dawson 2016).

3 Results

3.1 The spectral relationship between the NAO and WDs

It is known that the NAO can significantly affect the location and intensity of the North Atlantic jet stream (Woollings and Blackburn 2012). As the North Atlantic jet and the subtropical westerly jet (SWJ) over Eurasia are connected (Strong and Davis 2008), this implies the possibility of a physical pathway between the NAO and the Western Himalaya. Figure 1 shows the climatology of global subtropical jets during winter months and how that climatology changes during NAO+ and NAO− months.

Fig. 1
figure 1

Zonal wind speed at 200 hPa and location of the Northern Hemisphere westerly jet axis for a the Dec–Mar climatology and b the difference between NAO+ and NAO−. The jet axis is defined as a meridional local maximum in 200 hPa U wherever U > 30 m s−1. In b the jet axes are shown for both composite NAO+ (red) and NAO− (blue). Subplots on the right show the behaviour over South Asia

Full size image

The climatological subtropical jet has two distinct branches: a short one over North America and the North Atlantic, associated with North Atlantic storm tracks, and a longer one covering North Africa and South Asia, in which WDs are embedded. During NAO+ months, the North Atlantic jet moves northwards, responding to the strong meridional pressure gradient there. During NAO− months, the North Atlantic jet moves southward, almost merging with the subtropical jet over Africa. Despite this, the subtropical jet is more intense and less stable during NAO+ , likely owing to an increased upper-level meridional temperature gradient in the midlatitudes and subsequently stronger eddy forcing (e.g., Li et al. 2021), which plays a significant role in WD genesis. There is also a slight meridional shift in the jet axis as it passes over India, where it is about 200 km further north during NAO+ than NAO−.

It is known that both the latitude of the subtropical jet (Hunt et al. 2017) and its intensity (Hunt et al. 2018) exert controls on WD behaviour over South Asia. In addition, the greater vertical wind shear exerted by a stronger jet is associated with more intense updrafts in WDs (Baudouin et al. 2021). As the NAO affects both of these jet characteristics, we can reasonably hypothesise that it will also modulate downstream WD behaviour. We test this hypothesis in two ways, using both wavelet and compositing techniques. Figure 2 shows the frequency spectra for both the NAO index (a,b) and WD frequency (c,d) between 1950 and 2020, computed using a derivative-of-Gaussian wavelet. Because this wavelet only has a real part, the spectra highlight the locations of peaks and troughs (both have high power) across the range of frequencies. Black contours that mark the locations of significant peaks and troughs of the NAO are superimposed on the WD frequency spectrum so that potential regions of overlap can be identified.

Fig. 2
figure 2

Frequency spectra of a, b the NAO and c, d WD frequency over South Asia, computed using a derivative of Gaussian wavelet with an angular frequency of one month. Both the NAO index and WD frequency are standardised before computation, so that the spectral power is normed. Black line contours on the power spectra a, c indicate where the NAO spectral power is significantly greater than zero at a 90% confidence level. Hatching on the power spectra a, c indicates the cone of influence, where edge effects are significant. e Shows the p value of the correlation between the NAO and WD spectral power timeseries at each frequency. The relationship is plotted for a selection of lead times (0 to 10 months) between the two timeseries, such that WD frequency lags the NAO. The grey dashed line indicates a 5σ (p < 3 × 10–7) threshold

Full size image

Spectral power in the NAO index is largely concentrated into three bands: annual, 2–4 year, and 8–16 year. The sub-annual bands also contain some periods of significantly high power, and there may be power at much lower frequencies not detectable in the 70-year dataset used here (e.g. Trouet et al. 2009). The WD frequency spectrum is dominated by annual variability; this is expected since WDs are almost entirely absent during the summer monsoon (June to September) but occur at rates approaching ten per month during the winter. This annual cycle is also controlled by the subtropical jet, which migrates northward over China every summer in response to seasonal shifts in the Hadley Cell (Schiemann et al. 2009) and is not related to the NAO.

The two other bandwidths with significant NAO activity, 2–4-year and 12–16-year, show evidence of relationship with WD frequency. Peaks and troughs within the 2–3 year band align closely with those in the WD frequency spectrum, especially between 1950 and 2000 and most strongly between 1970 and 1990, which may partially explain why previous authors have noted declining correlations between the NAO and synoptic winter activity over the Western Himalaya since the 1980s. Interestingly, there are two periods where NAO power in this band weakens significantly, 1963–1966 and 2003–2007, during which the WD frequency power in this band also weakens significantly. There is a slowly varying lag relationship between the two: sometimes the extrema nearly coincide, whereas sometimes the extrema in the NAO spectrum lead the extrema in the WD spectrum by as much as six months. Since 1980, there has also been a more marked overlap in the extrema of the 8–16-year band of both datasets, although the power of both the NAO and WD frequency is weaker here than in the 2–4-year band.

We can test the potential predictive power of these relationships by quantifying the strength of correlation between the two timeseries as a function of lead time and wavelet frequency, as shown in Fig. 2e. We extract the timeseries of spectral power at each frequency for both the NAO and WD frequency and compute the p value of Pearson’s correlation coefficient. This is then repeated for lead times (NAO leading WDs) of up to ten months. The significance of these correlations exceeds a 5σ threshold at both 4-year and 10–12-year bands, both related to the high activity bands we saw earlier. However, in the 10–12-year band the relationship weakens as a function of lead time, implying that the NAO would not be a useful predictor of decadal variability of WDs, and that instead they share a common forcing at these timescales. However, in the 4-year band, the relationship between the NAO and WD frequency is negligible at zero lead time (i.e., when they are taken simultaneously) but increases as the NAO lead time is increased, reaching a maximum at about eight months. This implies the existence of a robust lead-lag relationship between the two that could be used to improve predictability of interannual WD variability. We also see a number of significant correlations at higher frequencies (i.e., shorter periods), typically peaking in strength at lead times between 0 and 2 months. Given the very large power of WD frequency in annual and sub-annual bands, quantifying and understanding these relationships is important.

The results of our wavelet analysis lend credibility to the idea that the NAO can modulate WD frequency through its interactions with the subtropical jet, and that this relationship offers an opportunity for seasonal weather predictions of the South Asian winter.

3.2 The role of western disturbances

Before we investigate predictability, we first quantify the role of WDs in supporting the teleconnection between the NAO and Western Himalayan winter precipitation. This is important because WDs are the likely conduit, as their frequency and intensity are sensitive to changes in the subtropical jet discussed earlier and they are also the primary source of winter precipitation to the Western Himalaya. As such, investigating WDs allows us to quantify the effect of the NAO on the local dynamics responsible for winter precipitation over North India.

We start with composite spatial maps of WD frequency and intensity and contrast them between NAO+ and NAO− periods. Although this method does not capture lead-lag relationships between the two, their extrema are sufficiently simultaneous that this should not be a problem. These composites are shown in Fig. 3. The highest values of WD frequency are found over the Western Himalaya, where the maximum value is about 4.3 season−1 (100 km)−2. Over the box [68–78°E, 27–38°N; marked in Figs. 3 and 4] which contains the Western Himalaya, the highest climatological WD frequency, and the highest mean winter precipitation, the mean WD frequency is about 2.6 season−1 (100 km)−2. During NAO+, the peak value is 31% higher than during NAO− (Fig. 3b), and the average frequency over the box is 20% higher. Following the slight northward shift in the subtropical jet during NAO+ (see Fig. 1), there is a marginal poleward shift in the location of highest WD frequency. There is also a significant change to the mean intensity: the average 350 hPa relative vorticity at the centre of WDs passing through the box is 7.7 × 10–5 s−1; this is 7% higher during NAO+ than NAO−.

Fig. 3
figure 3

The role of the NAO in modulating WD frequency and intensity. a Climatological frequency in Dec–Mar; b composite difference in frequency between NAO+ and NAO− months; c climatological intensity in Dec–Mar; and d composite difference in intensity between NAO+ and NAO− months. The region with the highest climatological winter precipitation and WD frequency surrounding the Western Himalaya is marked with a black rectangle in (a). Stippling in b and d indicate where the NAO+ and NAO− populations differ significantly at a 95% confidence level, according to Welch’s t test

Full size image
Fig. 4
figure 4

Effect of the NAO on Dec–Mar precipitation and vertically integrated moisture flux over South Asia. a Dec–Mar climatology, b composite difference between NAO+ and NAO− months and c the ratio of the difference and the climatology. The black rectangle in a marks the same region as in Fig. 3. Stippling in c indicates where a Welch’s t test indicates where the distributions of NAO+ and NAO− precipitation are significantly different at the 95% confidence level

Full size image

Both WD frequency and intensity exert strong controls on precipitation over the Western Himalaya and surrounding region by modulating southwesterly moisture flux from the Arabian Sea (Hunt and Dimri 2021). Figure 4 shows the climatology of winter precipitation and vertically integrated moisture flux over South Asia and how these fields change between NAO+ and NAO−. Mean precipitation within the black rectangle over the APHRODITE period (1950–2015) is 0.97 mm d−1 with a peak value of 4.22 mm d−1. During NAO+, these are 45% and 57% higher respectively than in NAO−. Much of this change can be explained by increased westerly moisture flux entering the box and impinging on the orography. The mean value is 47.8 kg m−1 s−1, which is 39.9% larger during NAO+ than NAO−. This results in a significant increase in precipitation across much of north India and Pakistan. Precipitation over some areas of Jammu and Kashmir is as much as 1.5 mm d−1 higher (about 50% of the climatological mean) during NAO+ than NAO−, and states with climatologically low winter precipitation—Rajasthan and Gujarat—see almost twice as much precipitation.

3.3 Developing a model of the NAO–precipitation teleconnection

Having established a relationship between the NAO and WDs that offers potential predictability for the latter, we now seek to extend this to winter precipitation over the Western Himalaya, from which we will explore potential teleconnection mechanisms.

3.3.1 A case study: Gulmarg, Jammu and Kashmir

We start with a simple model, linking the winter (NDJF) precipitation, measured by a gauge at Gulmarg meteorological station, with seasonal mean values of the NAO and ONI. Gulmarg is a popular skiing destination, located in the northwestern corner of the Kashmir valley in the Pir Panjal Range, at an elevation of 2940 m. Firstly, we fit a simple linear model of the form \(p = a_0 + a_1 {\text{NAO}} + a_2 {\text{ONI}}\), where the coefficients, \(a_n\), are found through a least-squares optimisation. Doing this (Fig. 5a), we find \(a_0 =\) 4.14 mm d−1, \(a_1 =\) 0.97 mm d−1, and \(a_2 =\) 0.30 mm d−1. The correlation coefficient between modelled and observed rainfall is 0.35, significant at a 95% confidence interval, with underestimated variance being responsible for much of the error. The model is only weakly improved by the inclusion of the ONI—the correlation between the seasonal mean observed precipitation and the NAO alone is 0.31. It is worth noting that the majority of winter precipitation over Gulmarg is snow, resulting in large uncertainties in gauge measurements (Strangeways 2004).

Fig. 5
figure 5

A simple teleconnection model for Nov-Mar precipitation over Gulmarg in Kashmir. a Gauge observations from Gulmarg station are given by blue bars; the mean value of the NAO index is given in red; the mean value of the SOI is given in orange, and black markers denote the predicted precipitation. b The same model is fit to other sources of precipitation: two reanalyses (IMDAA and ERA5), two gridded gauge datasets (IMD and APHRODITE) and a satellite dataset (GPM IMERG). The correlation coefficient between each pair of observed and modelled precipitation is given in the legend

Full size image

We repeat the process for a variety of precipitation sources in Fig. 5b, in each case, fitting the linear NAO/ONI model to simulate the seasonal mean precipitation. In each case, the precipitation is taken from the grid square containing Gulmarg (0.25° in width and length for ERA5, IMD, and APHRODITE; 0.12° for IMDAA; 0.1° for GPM IMERG) with no additional downscaling or post-processing. Using this larger area results in a smaller interannual variance, and in each case the correlation between the modelled precipitation and “observed” precipitation is better than for the individual gauge. This tells us two things: firstly, that it may be possible to use the NAO or related variables to construct a simple predictive model for winter precipitation over north India; and secondly, that such a model could benefit from considering precipitation averaged over a large area. We now discuss the construction of such a model with the aim of providing seasonal forecasts at a one-month lead time, starting with a more detailed analysis of variability over the North Atlantic.

3.3.2 Flavours of the NAO and a predictive model

To understand this teleconnection, and build our model, we need to start at the source—quantifying the effect of the NAO on jet variability over the North Atlantic and investigating how information about that variability is transported downstream through the subtropical jet. As Wang et al. (2012) demonstrated, there is more to NAO variability than can be represented through just its static dipole index. Following Filippi et al. (2014), we quantify this variability using a simple EOF analysis, taking daily 200 hPa winds from ERA-Interim during the winter months (December to March) over the whole dataset length (1979–2018) over the domain [80–0°W, 10–70°N] and conduct EOF analysis (Sect. 2.6) to extract the primary modes of variability.

The first four EOF loading patterns are given in Fig. 6. The first EOF explains 20.4% of the variance and is strongly linked to the NAO: the correlation between the monthly means of its principal component (PC) timeseries and the NAO is over 0.8. Large positive values of PC1 correspond to NAO+ and a split jet over the North Atlantic, whereas large negative values correspond to NAO− and a single jet. EOF2 explains a similar amount of variance and is also associated with the latitude of the jet over the North Atlantic. The two patterns are in quadrature and can therefore be taken as a pair that explains latitudinal variability of the jet. EOF3 and EOF4 similarly form a pair that explains tilting of the jet. For example, positive values of PC3 tend to produce a southwesterly jet that passes over the British Isles, reminiscent of Scandinavian blocking. Combined, these four EOFs provide a more complete quantification of North Atlantic variability than the NAO index alone, reflecting different ‘flavours’ of the NAO.

Fig. 6
figure 6

The first four EOFs of winter (Dec–Mar) 200-hPa zonal wind over the North Atlantic, computed using daily means over the domain shown. Fraction of variance explained is given in parentheses in the title of each subfigure

Full size image

Both the latitude and tilt of the jet over the North Atlantic are likely to play a role in the downstream stability of the subtropical jet as it passes over Eurasia. To quantify this relationship, and its potential role in the predictability of Himalayan precipitation, we construct bivariate distributions (Fig. 7) to investigate the relationship between each of the first four principal components of the EOFs described in Fig. 6 and four characteristics of the subtropical jet which we will now describe briefly. Firstly, the core of the subtropical jet is identified at daily frequency in the domain [0–80°E, 15–55°N] by finding the strongest zonal wind speed (and its latitude) at each longitude. ‘Jet strength’ is the monthly mean value of these core 200 hPa zonal wind speeds. Defining the jet as the region 750 km either side of the core at each longitude, we then compute ‘jet shear’ as the monthly mean difference of 200 hPa and 500 hPa zonal wind 500 hPa in the jet region. We define ‘jet waviness’ is the monthly standard deviation of meridional wind in the jet region. Finally, ‘jet latitude’ is computed as the mean latitude of the core between 60°E and 80°E (i.e., as it passes over South Asia). Together these comprise a toolkit that allows us to assess the strength, variability, and instability of the jet from the point of view of a disturbance passing along it. We also include mean precipitation over the state of Uttarakhand (situated along the Western Himalaya) for the following month, to assess whether these relationships could have predictive power.

Fig. 7
figure 7

Bivariate distributions of the principal components associated with the leading four EOFs of North Atlantic winter jet variability (y-axis) and characteristics of the subtropical jet and mean precipitation over the western Himalayan states one month later (x-axis). Distributions are coloured grey if the correlation is not significantly different from zero at the 90% confidence level. If the correlation is significant and positive, the distribution is coloured red; if negative, blue. Darker shades indicate higher frequencies. The correlation coefficient between each pair of timeseries is given in the top right corner of the respective subplot

Full size image

PC1, which is very closely related to the NAO index, is positively correlated with jet latitude, corroborating our composite analysis in Fig. 1. It is also positively correlated with waviness but has a negative correlation with shear. At a one-month lead time, it is positively correlated with precipitation over the Western Himalayan states. PC2, which adjusts the latitude of the jet over the North Atlantic, seems to play a smaller role in the behaviour of the subtropical jet over Eurasia: it is negatively correlated with its strength and weakly (though still significantly) positively correlated with its waviness. There is no significant correlation between PC2 and Uttarakhand precipitation a month later. PC3, associated with tilting of the jet over the North Atlantic, is positively correlated with the strength, shear, and waviness of the downstream subtropical jet, as well as with Uttarakhand precipitation a month later. PC4 is not significantly correlated with any of the selected characteristics of the subtropical jet, nor the precipitation.

In order to be able to propose a mechanism linking NAO flavours to Himalayan winter precipitation, we also need to understand how the jet characteristics above project onto the precipitation. At a one-month lead time (bottom row of Fig. 7), we find a significant negative correlation with jet latitude, and a significant positive correlation with jet shear (i.e., a ‘deeper’ jet). Precipitation responds much more strongly to these jet characteristics at zero lead time (not shown): the correlation with jet strength is 0.49, with jet latitude is − 0.35, and with jet shear is 0.19; but much more weakly to the four PCs. This supports our hypothesis that the subtropical jet responds slowly to changes in the NAO, resulting in a delayed response of the precipitation over the Western Himalaya. With positive values of PC1, it is the strong relationship with jet waviness that overcomes the weakening effects of reduced shear and increased (i.e., further north) latitude to provide increased Himalayan precipitation. With positive values of PC3, increased jet strength, shear, and waviness come together to increase Himalayan precipitation.

Now we know that these EOFs/PCs have potential predictive power for Uttarakhand precipitation, and a reasonable intuition for why—let us develop a simple statistical model to see if we can achieve (and subsequently quantify) similar results for other states of North and Central India. We do this by using a least-squares method to fit a linear model with an offset, i.e.:

$$p\left( {t + 1\, {\text{month}}} \right) = c + a_1 {\text{PC}}1\left( t \right) + a_2 {\text{PC}}2\left( t \right) + a_3 {\text{PC}}3\left( t \right) + a_4 {\text{PC}}4\left( t \right) + a_5 {\text{ONI}}\left( t \right) ,$$

where \(p\) is monthly state-averaged precipitation and \(c\) and \(a_n\) are the sought coefficients. For completeness, we also include the ONI term. Values of these coefficients for states where the correlation between modelled and observed precipitation is significant at the 99% confidence level are given in Table 1. We also use a bootstrap test—refitting the models several thousand times, each after shuffling the predictors along the time axis—to confirm that the models are not overfitted (not shown). We see that in each of these five states, the model captures a reasonable amount of the precipitation variability a month in advance—indeed the \(r\)-values actually fall if we use precipitation from the same month (not shown)—supporting our hypothesis that these PCs have some predictive power. \(a_1\) and \(a_3\) are typically the largest coefficients, suggesting, as before, that winter precipitation is most sensitive to changes in PC1 and PC3. We also note the consistently low and statistically insignificant values of \(a_4\) and \(a_5\), which points to neither PC4 nor the SOI being useful predictors at this lead time. Following our earlier discussion that such models may perform better predicting precipitation over larger regions, we also include the Western Himalayan states—i.e., the mean precipitation over Himachal Pradesh, Punjab, Jammu and Kashmir, Uttarakhand, Haryana, and Uttar Pradesh. This performs well (\(r\) = 0.323), but not as well as over Punjab or Himachal Pradesh alone.

Table 1 Coefficients c and an, following a least-squares fit to monthly mean precipitation over each state (at a one-month lead time)
Full size table

4 Conclusions

Winter (Dec–Mar) precipitation is vital to the Western Himalaya and surrounding region. At low elevations, it provides the rainfall needed to grow rabi crops; at higher elevations, it provides snowfall, which, when it melts later in the year, is a crucial component of the region’s water security. As such, identifying potential sources of predictability for seasonal and subseasonal precipitation is important for improving both preparation and mitigation. Many previous studies have established a statistical relationship between Western Himalaya winter precipitation and the NAO, but these disagree on the significance, duration, and even sign of the relationship, with none establishing a causal mechanism.

In this study, we used seventy years of reanalysis, precipitation, and NAO data to confirm that there is a strong positive correlation between the NAO and winter precipitation over the Western Himalaya, and to demonstrate the mechanism through which this relationship is realised. In months where the NAO index is greater than 1 (NAO+), the subtropical jet is significantly more intense than in months where the NAO index is less than − 1 (NAO−), likely because of a stronger upper-level meridional temperature gradient and greater eddy forcing upstream. This results in the production of more western disturbances (WDs), the winter storms that bring the majority of winter precipitation to the Western Himalaya. On average, WDs striking the Himalaya are 20% more common during NAO+ than during NAO−, rising to an increase of 31% in areas where WDs are most frequent. The more intense subtropical jet is also associated with an environment of greater baroclinic instability, meaning that WDs are about 7% more intense during NAO+ than during NAO−. Combined, these increases in WD frequency and intensity result in 40% moisture entering the region around the Western Himalaya, resulting in 45% more winter precipitation during NAO+ than during NAO−.

We used spectral analysis to explore the co-variability of the NAO and WD frequency and showed that they were strongly and significantly correlated on 4-year and 10–12-year timescales. By lagging the WD timeseries, we showed that the NAO had little potential predictability for WD frequency on the 10–12-year band—where the two timeseries instead probably share a common forcing—but good potential predictability up to several months ahead on the 4-year band. This implies that the NAO at least partially forces WD frequency on interannual timescales, and that the delay of a few months could be used to develop a simple seasonal forecast.

A simple model using just the NAO and ONI (an index quantifying the sign and strength of ENSO) was able to replicate precipitation at a single test gauge—Gulmarg in Kashmir—reasonably well but was unable to capture much of the interannual variability associated with such a small catchment. Verification of the model using other datasets over larger scales, where the interannual variability was lower, performed better.

To develop a more advanced, predictive model, we extracted additional modes of variability associated with the NAO by applying EOF analysis to 200 hPa winds over the North Atlantic. Four EOFs and their associated PC timeseries were extracted: the first pair describing the latitude of the jet over the North Atlantic, explain 37.4% of its variability; the second pair, describing the tilt of the jet, explain 18.1% of its variability. The PC timeseries associated with the first EOF (PC1) is significantly correlated with downstream waviness in the subtropical jet, whereas PC3 was significantly correlated with downstream jet strength. A linear model combining these PCs was used to predict mean monthly precipitation over five states of North India at a one-month lead time. The predicted precipitation was significantly correlated with observations in each case, suggesting that the model, and hence North Atlantic variability, has useful predictive skill for Western Himalayan precipitation on seasonal timescales. The model performs slightly better—albeit at a shorter lead time—than the one proposed by Nageswararao et al. (2016b), who used oceanic heat fluxes between August and November to predict seasonal precipitation over North India the following winter. Unlike some previous studies (e.g., Dimri 2013; Cannon et al. 2016; Kamil et al. 2019), we found the relationship between ENSO and winter precipitation over the Western Himalaya to be quite weak. This relationship is known to be volatile, and so our result could be due to the choice of region, baseline period, or even precipitation dataset. Either way, further research on this teleconnection is clearly required.

Further work is also needed to address the precise set of mechanisms by which the NAO modulates the subtropical jet in more detail, to quantify the lag between the NAO and Western Himalaya precipitation, and to quantify the impact of the NAO on extreme precipitation events in the region. Future research may also want to explore the relationship between the NAO and the northeast monsoon over Sri Lanka and the southern states of India.