Dominant controls of cold-season precipitation variability over the high mountains of Asia

A robust understanding of the sub-seasonal cold season (November–March) precipitation variability over the High Mountains of Asia (HMA) is lacking. Here, we identify dynamic and thermodynamic pathways through which natural modes of climate variability establish their teleconnections over the HMA. First, we identify evaporative sources that contribute to the cold season precipitation over the HMA and surrounding areas. The predominant moisture contribution comes from the mid-latitude regions, including the Mediterranean/Caspian Seas and Mediterranean land. Second, we establish that several tropical and extratropical forcings display a sub-seasonally fluctuating influence on precipitation distribution over the region during the cold season. Many of them varyingly interact, so their impacts cannot be explained independently or at seasonal timescales. Lastly, a single set of evaporative sources is not identifiable as the key determinant in propagating a remote teleconnection because the sources of moisture anomalies depend on the pattern of sub-seasonally varying dynamical forcing in the atmosphere.


INTRODUCTION
The High Mountains of Asia (HMA) include the Himalayas, Karakoram, Hindu Kush in South Asia, and the Pamir, Tian Shan in Central Asia, which are home to over 100 mighty peaks that reach the sky over 7000 meters above sea level. Except for the eastern and central Himalayas, a significant fraction of the annual precipitation over these mountains and surrounding valleys is received during the cold season (November-March), which regulates the largest non-polar glacial mass, also known as the Water Towers of Asia. The meltwater from these mountains serves as one of the significant sources of river recharge for several river basins in Asia, including the Brahmaputra in the eastern Himalayas 1,2 (>65%), the Ganges in the central Himalayas 1-3 (~48%), the Indus in the western Himalayas and Karakoram 4,5 (>60%), and the Amu-Darya and Syr-Darya in Pamir (>70%) 6 . The region encompassing the HMA is quite susceptible to high interannual climate variability 7,8 . To this end, while the runoff through meltwater is a relatively more sustainable source for river flows as glaciers partially buffer extreme dry conditions in year-to-year precipitation variability 9 , this source of freshwater is now considered threatened due to a persistent rise in global temperatures 10 . Consequently, economies in the downstream regions of Asian mountains that heavily depend on irrigated agriculture through meltwater-dependent river flows 11,12 and hydropower generation 9 are now considered vulnerable due to the uncertainty surrounding the future of seasonal snow accumulation and overall glacial mass in the coming decades. Therefore, a precise knowledge of processes dictating climate variability and change over the HMA is imperative for the sustenance and sustainability of regional social-ecological systems.
Most winter precipitation over the HMA comes through the eastward-moving extratropical cyclones, known as the Western Disturbances (WDs) 13,14 . Given that the HMA is in a landlocked region, WDs are the critical mechanism transporting moisture from distant water bodies and terrestrial areas. They interact with the Indian subcontinent's warm air mass and topography to produce precipitation that falls as rain over the lower elevations and valleys and snow over the higher elevations. The trajectories and the strength of WDs, which determine the amount of moisture distributed as precipitation across the mountains in Central and South Asia 14,15 , are influenced by the positioning and strength of the subtropical westerly jet over the region. Several earlier studies establish that natural climate fluctuations, including El Niño-Southern Oscillation (ENSO) and North Atlantic Oscillation (NAO), have a positive influence on the winter precipitation variability over this region [16][17][18][19][20] , linking it to the southward shift and strengthening of the subtropical westerly jet 21 and enhanced moisture transport from surrounding water bodies (Arabian, Mediterranean and Caspian Seas) 22,23 into the region during the positive phase of NAO [24][25][26] and ENSO 20,[23][24][25]27 .
However, recent studies demonstrate a lack of spatiotemporal robustness in these teleconnections 23 . For instance, one study suggests a dipolar pattern of NAO influence with a negative association over the Pamir Mountains and a positive association over the Karakoram and Himalayas 28 . Another study determines that direct ENSO teleconnection over this region lacks subseasonal persistence and suggests that the tropical Indian Ocean plays a crucial role in regulating the positive ENSO forcing via inter-basin interactions 23 . A role of contemporaneous atmospheric wave activity in the Atlantic and Eurasian regions that, in addition to NAO, gives rise to several recurrent atmospheric circulation patterns such as East Atlantic Mode (EAM), East Atlantic/West Russia pattern (EAWR), Scandinavian pattern (SCAND), has also been suggested in precipitation variability over the region comprising of the HMA 29 . Likewise, Siberian High (SH), a semipermanent feature in the Northern Hemisphere, also influences the precipitation distribution over the HMA 30 . However, a clear understanding of the collective influence of these modes of climate variability is currently lacking. In fact, these complementary and/or competing paradigms reflect knowledge gaps that need to be bridged for a more rigorous understanding of the cold season precipitation variability over the HMA.
One of the critical aspects of the cold season climate over the HMA, which has received less attention, is identifying complementary moisture sources that contribute to the regional precipitation. A few studies have taken on this topic using Eulerian water-source tagging capability in a Global Climate Model (GCM) 31 or employing a Lagrangian modeling approach on a few years of reanalysis data 32 . However, method-based uncertainties in moisture source identification, such as those arising from the use of single simulated or reanalyzed data, different water vapor tracking approaches, and relatively short analysis periods warrant further studies to establish a comprehensive understanding. Moreover, these studies have primarily focused on Central Asia and do not fully cover the region encompassing the HMA. More importantly, no study thus far has systematically evaluated thermodynamic pathways through which the abovementioned modes of natural climate influence precipitation variability over the HMA. Due to the lack of consensus among previous studies regarding their influences, a reevaluation of remote teleconnections over this region is also necessary on a sub-seasonal timescale.
Therefore, the primary aims of this study are 1) to establish a more robust understanding of evaporative sources that contribute to the cold season precipitation over the HMA and surroundings areas, 2) to reevaluate remote teleconnections on the subseasonal scale that exert an influence on precipitation variability over the HMA during the cold season, and 3) to explain remote teleconnections in terms of anomalies in atmospheric circulations and moisture transport. For this purpose, this study employs a Lagrangian moisture tracking algorithm to identify primary evaporative sources and examines their sub-seasonal variability over the HMA at the climatological scale in multiple reanalysis datasets. Furthermore, a multi-linear regression analysis is used to elucidate the relative influences of various modes of natural climate variability on the cold season precipitation over the HMA and unravel dynamic and thermodynamic pathways through which their teleconnections are established.

Moisture sources
Most cold season precipitation over the HMA is distributed in the western Himalayas, Karakoram, Hindu Kush, and Pamir Mountain ranges ( Supplementary Fig. 1). The evaporative source domain of the current study (Fig. 1a) provides up to 97% of the moisture for this precipitation, which is shown as the recycling ratio or the fraction of moisture falling as precipitation in a specific grid cell that originated as evapotranspiration from the entire source region (Fig. 1b). Note that these estimates are based on the averages of the European Center for Medium-Range Weather Forecast (ECMWF), Re-Analysis-Interim (ERAI) 33 , ECMWF Re-Analysis version 5 (ERA5) 34 , and National Aeronautics and Space Administration (NASA) Modern-Era Retrospective Analysis for Research and Application Version 2 (MERRA2) 35 , while their individual magnitudes are shown in Supplementary Fig. 2, which display close correspondence. The average recycling ratios based on regional contributions are shown in Fig. 2 for the most prominent sources. Supplementary Fig. 3 shows the individual contributions from all regions separately for ERAI, ERA5, and MERRA2. Every sub-region provides moisture, but the significant contribution comes from the Mediterranean/Caspian Seas (>25%), Mediterranean land (>24%), Indian Ocean (>14%), local recycling (>11%), and North Atlantic Ocean (>9%). The strongest influences of the Mediterranean/Caspian Seas, Mediterranean land, and North Atlantic Ocean are displayed in the west-northwest, with 42%, 36%, and 14% contribution to the total seasonal precipitation. The input from these three regions gradually decreases to 5% going west to east. On the other hand, the most robust footprint of local recycling is seen over eastern Nepal (up to >50% of the seasonal total) with a gradual decrease in contribution to 25% and less while going west over the Hindu Kush mountains and beyond. The Indian Ocean contribution is the strongest over relatively dry regions in southwestern Pakistan and the Uttar Pradesh region in India, contributing up to 35% of the total seasonal precipitation. The spatial distribution of regional contributions displays close correspondence across three reanalyses, with most of their differences within ±5%. The most noticeable disagreement exists in the case of local recycling, where its relative contribution is >10% lower in MERRA2 than that in ERAI and ERA5 over the Himalayas and Karakoram (Supplementary Fig. 3 Fig. 1 The evaporative source region, seasonal (November-March) recycling ratio, and recycled precipitation. a The evaporative source region with its subdivisions. b Seasonal fraction (expressed in %) of moisture falling as precipitation at each grid cell that originated as evapotranspiration (filled contours) over the evaporative source region (recycling ratio) and corresponding recycled seasonal precipitation (mm/day; line contours). This analysis shows an average of ERA5, MERR2, and ERAI from 1981 to 2018.
Moisture contributions from individual sources exhibit subseasonal variations, while their total contribution remains stable throughout the season, represented as the recycling ratio (0.89-0.96) (Fig. 3). Some disagreements in the sub-seasonal magnitudes of regional contributions also exist across the reanalysis datasets. Over half of moisture comes from the Mediterranean/Caspian Seas and Mediterranean land throughout the cold season. The moisture from the Mediterranean/Caspian Seas reduces to almost half (~20%) of their contribution to the early season (~35%). In contrast, an increase of 15-20% is seen from the Mediterranean land by the end of the season. Up to a 50% decline in contribution from the Indian Ocean is also  witnessed from its peak in January (~14%). The local recycling is also the strongest towards the end of the season (~18%), while moisture contribution from North Atlantic Ocean remains relatively steady after November (~10%). The evaporation from the Mediterranean and adjacent seas exhibits gradual decreases in the cold season with an annual minimum around April-May 36 , which explains the steady decline in moisture contribution from the Mediterranean/Caspian Sea during the cold season. Alternatively, a gradual increase in the contribution from extratropical land regions is consistent with their annual cycle of surface temperatures. Noticeable disagreements among reanalysis are present for moisture contributions from the Mediterranean/ Caspian Seas and local recycling (up to ±4%). Nonetheless, all datasets display a similar pattern of relative contributions and subseasonal variability from each source region (Fig. 3). Given that there is no significant disparity among the reanalysis (Fig. 3, Supplementary Figs. 2, 3), we perform the rest of the analyses using the average of three datasets.

Influence of natural modes on precipitation variability
Several natural modes of climate variability (see Methods) exert influences on year-to-year variations in monthly precipitation amounts during the cold season, as the coefficient of determination (R 2 ) magnitudes and areas with significant F-test statistics predominantly display gradual improvements throughout the season, highlighting the importance of each predictor in elucidating the drivers of precipitation variability over the study region ( Supplementary Fig. 4). However, the explained variance is relatively low in February and, to some extent, in January over most of the study area. In the Multi Linear Regression (MLR) model for total recycled precipitation, the relative comparison of standardized partial regression coefficients (β * ) across predictors and months reveals several important details (Fig. 4). First, the Tropical Western/Eastern Indian Ocean precipitation dipole (TWEIO) is the most prominent and persistent forcing throughout the cold season that exhibits a statistically significant positive influence on precipitation variability in several months, with the strongest teleconnection in December (|β * | > 0.6). March is the only exception when TWEIO's effect over parts of the Tien Shan mountains is significantly negative. SCAND is the other forcing that exhibits relatively persistent (though not always significant) positive impact throughout the season, except for the parts of the Hindu Kush mountains in November and the central Himalayas in January. EAM exerts the strongest, widespread, and statistically significant negative influences in December and January over the study region (|β * | > 0.6). Except for the Himalayas, EAM teleconnection is also negative in February. SH footprints are conspicuously comparable to EAM but relatively subdued during these three months, except in February over the Himalayas. Other modes exhibit sparingly noticeable but spatiotemporally varying influences on the recycled precipitation. Overall, the mix of spatiotemporally varying tropical and extratropical influences on the precipitation distribution over the HMA highlights the inherent complexity of precipitation teleconnections, which perhaps is one of the causes of the limited seasonal predictability skill over this region 37 .

Dynamical and thermodynamical pathways of remote teleconnections
The thermodynamic pathway through which a natural climate mode establishes its remote influence implies identifying moisture sources that provide excess or less moisture during an active phase of that mode. Likewise, the dynamic pathway means identifying atmospheric anomaly patterns that facilitate the propagation of a remote influence and make conditions conducive for significant precipitation variations. Therefore, subseasonal thermodynamic and dynamic pathways through which tropical and extratropical forcings exert influence can be explained by examining monthly anomalies in moisture contributions from respective evaporative sources and overlying circulation anomalies causing anomalous moisture transport. Note that the key feature of atmospheric circulations over the HMA during the cold season is the subtropical westerly jet, which guides the movement of WDs over the region 38,39 . The changes in the strength or positioning of the jet stream and near-surface wind patterns can be manifested using the anomalies in geopotential height at 200 mb and upper (200 mb) and lower (850 mb) atmospheric winds. Likewise, changes in the moisture contribution from source regions can be explained using the anomalies in their contributions to the total recycled precipitation over the study region. Consequently, to provide a physical explanation of spatiotemporally varying influences exerted by each mode of natural climate variability, we construct separate MLR models for 200 mb geopotential height, 200 mb and 850 mb zonal and meridional winds, and contributed precipitation from individual source regions, using the same set of predictors employed in the total recycled precipitation MLR (see Methods  Fig. 3 Climatological sub-seasonal variability of the recycled precipitation from individual source regions. Total recycling ratio (top line; right y-axis; unitless) and percent contribution from individual source regions to the total recycled precipitation (bottom six lines; left y-axis; %) at sub-seasonal (5-day running mean) timescale, averaged over the study region shown as black in Fig. 1a.
The lines represent the average of three reanalyses, and the shaded area represents the range across them. Regional contributions to the total recycled precipitation are only shown for regions with contributions exceeding >6% during the season.
few prominent cases in the main figures (Figs. 5-7) to elucidate the significant findings. Across all the modes and months, TWEIO exhibits the most substantial positive teleconnection in December (Fig. 4). A statistically significant strengthening of the subtropical westerly jet ( Fig. 5a; negative geopotential height anomalies) is witnessed in December with the anomaly center over the study region that stretches in semi northwest-southeast orientation. The anomalous [TWEIO] [SH] [NAO] [EAM] [  atmospheric wave train pattern also intensifies the southeastward moisture flow over the subtropics. As a result, these atmospheric anomalies lead to a statistically significant increase in moisture supply from the mid-latitude source regions, including the Mediterranean/Caspian Seas, Mediterranean land, and northern Eurasia, in addition to enhanced moisture advection from the Indian Ocean ( Supplementary Figs. 6-7). Local recycling also increases due to excessive precipitation (Fig. 5d). The strengthening of the subtropical westerly jet is also present in January. An anomaly center is now located west of the study region, stretching diagonally over the Arabian Peninsula and northern Africa ( . Notably, the enhanced moisture transport from the African land is spatially widespread and broadly statistically significant (Fig. 5n). The most striking aspect to note in these results is that no single set of moisture source regions is identifiable as the key driver of TWEIO influence over the HMA. The pattern of atmospheric anomalies dictates where the excess moisture would be coming from. The atmospheric anomaly pattern varies month to month (Fig. 5, Supplementary  Figs. 6, 7), regulating the moisture anomalies from individual source regions (Fig. 5, Supplementary Fig. 5), which dictate the spatiotemporal variability in the TWEIO influence. EAM displays the most substantial and widespread negative teleconnection in December and January (Fig. 4). In January, two statistically significant high-pressure centers emerge, one over the study region and one over central Europe. The latter one is zonally  Same as in (a-g) but for January. Standardization of regionally recycled precipitation is with respect to the total recycled precipitation. Stippling represents statistical significance at a 90% confidence level.
elongated over the Atlantic sector Fig. 6a). These EAM-induced atmospheric anomalies mutually limit moisture advection over the HMA from all source regions ( Fig. 6b-g, Supplementary Fig. 7), significantly reducing monthly recycled precipitation (Fig. 4). A similar but slightly shifted and less intense pattern of atmospheric anomalies also emerges in December (Supplementary Figs. 6, 7), causing a statistically significant but relatively less intense precipitation reduction (Fig. 4). While the overall structure of EAM-induced zonally elongated high-pressure cell over Europe persists throughout the cold season, the strength and location of atmospheric high over the study region substantially vary during the rest of the season, leading to insignificant or spatially heterogeneous precipitation responses in other months (Fig. 4, Supplementary Figs. 6-8). The pattern of SH teleconnection in January exhibit peculiar similarity with that of EAM (Fig. 4), which is a result of correspondence in the dynamic and thermodynamic anomalies induced by these two modes of climate variability ( Supplementary Figs. 6-9). As previously noted, SCAND is the only other forcing after TWEIO that exhibits relatively persistent (weak positive) teleconnection throughout the season (Fig. 4). The atmospheric anomalies induced by SCAND tend to strengthen westerly flow from the Mediterranean Sea to Russia, which diagonally traverses western parts of the study regions (Fig. 6f, Supplementary Fig. 6). This pattern of circulation anomalies fundamentally persists throughout the season with some variations in magnitude that explain the relative persistence of SCAND teleconnection. Considering the SCAND influence in February, two zonally connected low-pressure cells, one over the east Atlantic and one over the northwest of the HMA (Fig. 6f), increase moisture advection from midlatitude evaporative sources over the Hindu Kush and Pamir Mountains through the eastern flank of the atmospheric low (Fig. 6g-j, Supplementary Fig. 7 6 Standardized partial regression coefficients associated with the EAM and SCAND in the multi-linear regression models. Standardized partial regression coefficients of 200 mb geopotential height, zonal and meridional winds (a), and regionally recycled precipitation (b-g) on the EAM index in January. Six selected evaporative source regions exhibit the most prominent statistically significant relationship. (h-n) Same as in (a-g) but for SCAND in February. Standardization of regionally recycled precipitation is with respect to the total recycled precipitation. Stippling represents statistical significance at a 90% confidence level.
positive teleconnection of SCAND over those regions. Alternatively, moisture intrusion from the subcontinent land and the Indian Ocean is limited over the western flank of the atmospheric low. Given the northwest location of the SCAND-induced atmospheric low during the cold season over the study region, its positive influence is less noticeable over the Himalayas than over the other areas in the northwest (Fig. 4, Supplementary Fig. 10). Likewise, EAWR teleconnection is more significant over the Karakoram, Hindu Kush, and Pamir Mountains than over the rest of the study region (Fig. 4). January is the only exception when it exerts significant negative influence over the Himalayas. The main features of EAWR-induced atmospheric anomalies exhibit a southeast-northwest dipole of high-and low-pressure cells over the study region that display spatiotemporal variability in their patterns and strength ( Supplementary Figs. 6, 7). Interestingly, unlike other natural modes of variability where sources of moisture anomalies may vary sub-seasonally, the EAWR-induced influence primarily propagates over the HMA through anomalies in the climatologically leading moisture source regions (Mediterranean/Caspian Seas, Mediterranean land; Supplementary Fig. 11). The most prominent influences of NAO and ENSO are displayed in November (Fig. 4), which are also spatially counteracting. NAO establishes a trough (ridge) over Central Asia (India) (Fig. 7a) in November, facilitating enhanced (reduced) moisture transport from the Indian Ocean, East Africa, and mid-latitude regions to the north-northwest parts (central Himalayas) (Fig. 7b-g, Supplementary Fig. 7). In contrast, ENSO establishes negative teleconnection by inducing an opposite pattern of atmospheric and moisture anomalies in November (Fig. 7h-n). With few exceptions, the ENSO and NAO influences are relatively unremarkable during the rest of the cold season (Fig. 4, Supplementary Figs. 12, 13). Overall, the remote teleconnections over the HMA are a combination of Standardized partial regression coefficients associated with the NAO and ENSO in the multi-linear regression models. Standardized partial regression coefficients of 200 mb geopotential height, zonal and meridional winds (a), and regionally recycled precipitation (b-g) on the NAO index in November. Six selected evaporative source regions exhibit the most prominent statistically significant relationship. (h-n) Same as in (a-g) but for ENSO in February. Standardization of regionally recycled precipitation is with respect to the total recycled precipitation. Stippling represents statistical significance at a 90% confidence level.
compounding or counteracting influences that are often spatiotemporally superimposed on each other.

DISCUSSION
The interplay of tropical forcing and several internal modes of climate variability in the Northern Hemisphere plays a key role in precipitation variability over the HMA. Throughout the cold season, a statistically significant positive correlation exists between the two tropical forcings, i.e., ENSO and TWEIO ( Table 1). Note that TWEIO reflects the strength of the east-west precipitation anomaly dipole in the tropical Indian Ocean in response to disruption in the Walker circulation in the Indian Ocean sector 23,40,41 . It should not be confused with the Indian Ocean Dipole (IOD), which represents the dipole pattern of sea surface temperature anomalies in the Indian Ocean and is primarily active during the summer and fall seasons. A potential role of anomalous convection over the tropical Indian Ocean in cold season precipitation variability over the parts of HMA has also been highlighted in earlier studies 21,42,43 . More recently, studies have reported a significant relationship between TWEIO and sea surface temperature anomaly in the equatorial central Pacific, suggesting that TWEIO is likely an outcome of ENSO-forced interbasin interactions via atmospheric link 23,40 . Such interannual modality of tropical climate variability induces anomalous diabatic heating in the atmospheres over the equatorial central Pacific and the tropical Indian Ocean, as precipitation distribution deviates from its mean in these parts of oceanic basins. As a result, eastward propagating Rossby waves are generated with origins in the equatorial central Pacific and the tropical Indian Ocean that impact extratropical circulations 40,41 . These waves interact with each other and the internal atmospheric modes of variability such as NAO, EAM, etc 29 . These interactions during the Northern Hemisphere winter are reflected in the statistically significant monthly correlations between tropical forcing (ENSO, TWEIO) and several of these modes used in this study. When monthly sea level pressure anomalies are regressed onto the ENSO and TWEIO indices, the spatial patterns show that ENSO projects on NAO in November, EAM in November-January, and EAWR in December.
Similarly, TWEIO also projects on these modes in several months (Fig. 8). These projections are evident because, in each case, the regressed SLP pattern displays a close resemblance with the corresponding atmospheric mode of variability ( Supplementary  Fig. 14). Note that sea level pressure regressions show a similarity between ENSO and TWEIO for all months, as they exhibit persistent high correlations. These patterns reveal several pathways through which ENSO can influence a distant region. In the case of HMA, ENSO teleconnection is a combination of direct and indirect influences, where indirect effects come through interbasin connection (TWEIO) and interactions of ENSO and TWEIO with internal modes of variability in the Northern Hemisphere, including NAO, EAM, and EAWR. The influence of SH over the HMA also comes through multiple pathways, as it exhibits interplay with other internal atmospheric modes of variability in the Northern Hemisphere. The tropical forcings (ENSO, TWEIO) and SH operate relatively independently, as their indices lack significant correlations (Table 1). However, strong interactions between SH and other internal modes of variability (EAM, EAWR, and SCAND) are displayed through statistically significant correlations in several months. Monthly patterns of sea level pressure regressions onto EAM and EAWR ( Supplementary Fig. 14) reveal a low-pressure pattern over Mongolia and Siberiathe vicinity of SH maximathat explains SH's significant out-of-phase relationships with these modes in several months. The sea level pressure regression onto SH also shows that it negatively projects on EAM and EAWR throughout the cold season (Fig. 8). In both cases, the regressed SLP pattern resembles the negative phase of the corresponding atmospheric mode of variability ( Supplementary Fig. 14). A weak (significant) positive projection of SH on NAO (SCAND) is also displayed in November (December). Given the spectacle of interactions within the extratropics and between the tropics and extratropics, it is evident that the impacts of tropical forcings and internal atmospheric modes of Northern Hemisphere climate on the HMA precipitation cannot be explained in isolation. What makes it even more intriguing is that several of these modes of climate variability can co-exist during the cold-season 29 , and when superimposed, they can either compound or compensate each other's impact.
The lack of spatiotemporal persistence in several teleconnections also indicates that previously established scientific perspectives on the role of natural modes of climate variability, such as ENSO and NAO, in seasonal precipitation variability over the HMA are primarily equivocal. The analyses described in this study demonstrate that an individual mode can have opposing influences within a season and that its teleconnection may be a combination of multiple varying pathways. Moreover, unlike earlier studies that suggest a crucial role of moisture transport from the Arabian, Mediterranean, and Caspian Seas in establishing a teleconnection, our analysis demonstrates that a single set of evaporative sources are not identifiable as the sole determinant in propagating a remote teleconnection over the HMA. Note that an atmospheric anomaly, such as wave activity, alone does not tell us where the excess or deficiency of moisture originates. Just by looking at the pattern of low or high-pressure anomalies, a source of moisture anomaly cannot be established.
Because this study quantifies the amount of moisture coming from various oceanic and terrestrial moisture sources through back trajectory analysis, one can identify the exact origins of moisture anomalies in response to each natural climate fluctuation. The lack of persistence in teleconnections arises due to subseasonally varying dynamical forcing in the atmosphere, which also induces variations in sources of moisture anomalies. While moisture anomalies from the two most prominent sources (Mediterranean/Caspian Seas and Mediterranean land) are vital to establishing these relationships, complementary or contrasting contributions from other moisture sources often play an essential role. For instance, note that moisture from Africa, tropical Atlantic, and northern Eurasia contributes only~5% to the total at the  Supplementary Fig. 3), but in several cases, excessive moisture transport from these regions in response to a remote teleconnection is notable and significant ( Supplementary Figs.  5-13). An important aspect that may not be ignored is the existence of a relatively large unexplained variance in February (Supplementary Fig. 4). While our analyses are limited in what can be inferred from them in this regard, they provide notable clues that may be useful for more targeted analyses in future research endeavors. A large portion of the explained variance during the cold season is because of the tropical forcings. This role of tropical forcings is absent in February, as ENSO and TWEIO display little influence on the precipitation variability over the HMA 23 . Also, note that February is the only month with no significant correlation between tropical forcings and mid-latitude atmospheric modes of the Northern Hemisphere, suggesting the importance of tropical-extratropical interactions in propagating planetary-scale teleconnections. Moreover, our investigation is limited to contemporaneous relationships, but some earlier studies have noted the possibility of lagged teleconnections, such as the significant positive correlation between prior season sea surface temperature anomalies in the western tropical Indian Ocean and the cold season precipitation over northwestern South Asia 43 , which requires further exploration. There can also be a potential role of sub-seasonal variability of moisture sources in varying teleconnections; however, the current analysis is limited in covering this context. Nonetheless, our analyses provide a more comprehensive evaluation of sub-seasonal precipitation variability over the HMA during the cold season. Specifically, we establish a robust understanding of the terrestrial and oceanic moisture sources over the HMA during the cold season and demonstrate a crucial role of moisture recycling over the mid-latitude regions, which is also partly supported by findings in recent studies 29 . Further, we note that several tropical and extratropical modes of variability have a sub-seasonally fluctuating influence on the cold season climate over the HMA. Given that these modes varyingly interact, their impacts should not be studied individually or at seasonal timescales. An individual mode of climate variability can have opposing influences within a season, as its teleconnection is often a combination of multiple varying pathways. The dynamic and thermodynamic pathways dictating the variability of cold-season precipitation over the HMA are admittedly as complex as the region itself. While our analyses partly unravel their complexity, further research and modeling studies are needed to characterize the underlying mechanisms fully. Half of the moisture recycles from terrestrial sources during the cold season (Figs. 2, 3, Supplementary Fig. 3  land-atmosphere interactions that should be systematically explored in future studies.

METHODS Datasets
We use three reanalysis datasets: the ECMWF ERAI 33 , available at 0.75°× 0.75°horizontal grid spacing, the NASA MERRA-2 35 , available at 0.5°× 0.625°horizontal grid spacing, and the ECMWF ERA5 34 , available at 0.25°× 0.25°horizontal grid spacing. The ERA5 and MERRA2 are remapped to the ERAI horizontal grid using bilinear interpolation before their use. We obtained four variables from each reanalysis for moisture back trajectory analyses: daily evapotranspiration, daily precipitation, six-hourly vertically integrated eastward moisture flux, and six-hourly vertically integrated northward moisture flux. For the rest of the analyses, we obtain six monthly variables from each reanalysis: precipitation, sea level pressure, geopotential height, zonal winds, meridional winds, and sea surface temperature. All datasets cover 1980-2018. Note that the use of reanalyzed precipitation in moisture back trajectory analysis is necessary for it to be consistent with the reanalyzed moisture fluxes and evapotranspiration used in the determination of evaporative sources, as observed precipitation is often not well aligned with the reanalysis at the daily timescale. An inaccurate understanding of the relative contribution from various evaporative sources may be established if significant disparities exist between reanalyzed and observed precipitation. The use of multiple reanalyses minimizes the possibility of such inaccuracies. Apart from that, most of the HMA region is infamous for the low density of ground station network 44 and poor satellite-based enhancements in precipitation estimates 45 . Therefore, there is no way to conclude that gridded daily observations are superior to the reanalyzed precipitation over the HMA. In fact, the three datasets that provide observed daily precipitation for this region have less agreement among themselves relative to the consensus among the reanalysis (Supplementary Fig. 1).

Moisture back trajectory analyses
We use a modified version of the Dynamic Recycling Model (DRM)a Lagrangian moisture tracking algorithmto determine the moisture sources contributing to the cold season precipitation over the HMA. Previously, we have employed it in several studies over South Asia 46 , Mediterranean 47 , and Arabian Peninsula. The model was initially developed by Dominguez et al. (2006) 48 , which was later modified to remove reliance on parameterization [48][49][50] . The version used in the current study is detailed in Mei et al. (2015) 46 . It assumes a well-mixed atmosphere and uses vertically integrated zonal and meridional moisture fluxes. This assumption is reasonable over regions with low vertical shear, such as the HMA during the cold season. Note that data assimilation in reanalysis may lead to a violation of the moisture conservation assumption. However, Dominguez et al. (2006) 48  For this study, we make use of the evaporative source region between 60°W-90°E along the longitudes and 0°N-60°N along the latitudes, which is divided into nine exclusive but complementary moisture sources (Fig. 1a). The sub-division broadly separates known oceanic basins and large water bodies to ensure relatively accurate attribution of precipitation over the study region to remote sources. Most subdivisions of terrestrial areas are based on the homogeneity of precipitation seasonality within a sub-region during the November to March season. The nine subdivisions are defined as:1) Study regionland in Central and South Asia that receives a significant amount of cold-season precipitation, 2) Indian Oceangrids between 30°E-90°E and 0°N-30°N that include the Indian Ocean, Arabian Sea, Bay of Bengal, Persian Gulf and Red Sea, 3) tropical Atlanticthe Atlantic Ocean between 0°N-20°N and 60°W-30°E, 4) Mediterranean/ Caspian Seaswater grids that cover the Mediterranean Sea, Black Sea, Caspian Sea, and the Aral Sea, 5) Mediterranean landland between 25°N-40°N and 20°W-60°E, 6) Africa-land between 0°N-25°N and 20°W-60°E, 7) Subcontinent landland grid points in India and Pakistan that are not part of the study region, 8) northern Eurasialand grids up to 60°N, covering Europe (excluding Mediterranean land region), Russia and China, and 9) North Atlanticocean grid points between 20°N-60°N and 60°W-20°E.
We describe moisture source analyses using recycling ratio and recycled precipitation. The recycling ratio over a grid point is the fraction of moisture falling as precipitation that originated as evapotranspiration from the entire source region. Its average over the whole evaporative domain represents the overall recycling ratio. In other words, the multiplicative product of recycling ratio and precipitation provides recycled precipitation, i.e., the contribution of evapotranspiration within a region to precipitation in the same region. When the evaporative source domain is divided into several parts to understand the terrestrial and oceanic moisture sources over a sub-region, the evaporative contribution of grid points within that sub-region to the precipitation in the same region is referred to as local recycling. The total recycled precipitation over the sub-region is now a combination of local recycling and moisture advection from other parts of the evaporative source domain. Therefore, the recycling ratio over a subregion is the aggregate of local recycling and recycled moisture from different sub-regions within the evaporative source domain.

Multi-linear regression analyses
We use seven internal and/or inter-annual modes of natural climate variability in a MLR to investigate the variability of cold season recycled precipitation. These modes are defined as 1) ENSOthe standardized Niño 3.4 index, which is the area-averaged detrended sea surface temperature anomaly over the equatorial Pacific Ocean (5°S-5°N, 170°-120°W), 2) TWEIO 23,40the standardized differences of detrended area-averaged precipitation over the tropical western (40-80°E; 10°S-10°N) and eastern (90-140°E; 10°S-10°N) Indian Ocean, 3) SHthe standardized areaaveraged detrended sea level pressure over 40-65°N; 80-120°E, 4-7) the standardized indices of NAO, EAM, EAWR and SCAND from the Climate Prediction Center (https://www.cpc.ncep.noaa.gov/data/teledoc/ telecontents.shtml), which are based on the rotated principal component analyses of mean standardized 500 mb geopotential height anomalies over 20-90°N latitudinal band. The spatial pattern of the last four indices is shown in Supplementary Fig. 14 using the monthly sea level pressure regression onto each of these indices. All variables used in these analyses are based on the average of ERAI, ERA5, and MERRA2. We employ Variance Inflation Factor (VIF) analysis to determine the multicollinearity in the model and note that the MLR model does not suffer from this issue as VIF corresponding to all independent variables remains less than 1.5. The importance of each of the seven predictors is determined by comparing the R 2 and F-test statistics (Supplementary Fig. 3). The same seven predictors construct separate MLR models for 200 mb geopotential height, 200 mb and 850 mb zonal and meridional winds, and nine precipitation contributions from individual evaporative source regions. The standardized regression coefficient (β * ) for each independent variable in an MLR model is found by multiplying the unstandardized regression coefficient with the ratio of the standard deviation of the independent variable and dependent variable. In the case of those MLR models where contributed precipitation from one of the nine evaporative source regions is the dependent variable, β * for each independent variable is found by multiplying unstandardized coefficients by the ratio of the standard deviation of the independent variable and total recycled precipitation. Using total recycled precipitation instead of regionally contributed precipitation for standardization ensures that β * for each independent variable is comparable across all MLR models in question.

CODE AVAILABILITY
The scripts developed to analyze these datasets can be made available on request from the corresponding author.