Famines and likelihood of consecutive megadroughts in India

Consecutive failures in the summer monsoon rainfall led to widespread and severe droughts with profound implications for agricultural activities in India. However, the likelihood of successive megadroughts in India’s past and future climate remain poorly understood. Using Palmer Drought Severity Index (PDSI) from the Monsoon Asia Drought Atlas (MADA), we show that the major famines that affected millions of people during 1200–2018 were linked with summer monsoon droughts. Four megadroughts covering more than 40% of the country occurred for two consecutive summer monsoon seasons during 1200–2018. The most recent and severe megadrought occurred in 2002–2003. Simulations from the Community Earth System Model (CESM) for the last millennium (850–2005) ensemble (LME) show that the likelihood of two and three-year consecutive megadroughts during the summer monsoon is about 0.7 and 0.3 events per 100 years, respectively. Large ensemble simulations from CESM (CESM-LE) show a decline in the frequency of megadroughts in the future. Summer monsoon megadroughts are strongly associated with the warm sea surface temperature (SST) anomalies over the Pacific Ocean in the past and future climate. Substantial warming under the projected future climate can cause megadroughts under near-normal precipitation during the summer monsoon season. Despite the projected decline in the likelihood of the summer monsoon megadroughts under the warming climate, megadroughts in the future can have considerable implications for India’s food production and water availability.


INTRODUCTION
Indian summer monsoon (June-September) accounts for more than 80% of the total annual rainfall 1 . Rainfall during the summer monsoon plays a vital role in the Indian economy 2,3 . The Indian summer monsoon has a year-to-year variability of about 10% of its long-term mean 1,2 . A rainfall deficit of more than 10% is considered a drought, while rainfall of more than 10% of the long-term mean is termed surplus 2 . A deficit of 19-21% in rainfall in 2002 caused a drought that resulted in an estimated loss of billions of dollars, which was about 1% in the Gross Domestic Product (GDP) of India 4 . From 1901 through 2010, 17% of the total years experienced drought during the summer monsoon season. Characteristics and the occurrence of droughts during the summer monsoon season are well studied 1,[4][5][6][7][8] . About 18 meteorological droughts occurred in India during the 1870-2018 period 7 . The severe meteorological and hydrological droughts occurred in 1876, 1899, 1918, 1965, and 2000 7 . While droughts due to summer monsoon failure have implications, consecutive droughts disrupt India's socio-economic well-being. For instance, the drought of 2015-2018 caused groundwater depletion and affected about one-fourth of the total Indian population 7,8 . Similarly, the drought of 1875-1876 caused the great famine and affected millions of people in India 9, 10 .
The year-to-year variability of the Indian summer monsoon is influenced by El Nino Southern Oscillation (ENSO [11][12][13][14][15][16]. Warm sea surface temperature (SST) anomalies over the Pacific Ocean are negatively related to the Indian summer monsoon rainfall (ISMR). The negative relationship between summer monsoon rainfall and SST anomalies is the strongest over the eastern Pacific region 15 . The physical mechanism of the linkage between ENSO and the Indian summer monsoon is associated with atmospheric circulation and ENSO-induced SST variability 17 . The coupled atmospheric and oceanic variability modify Walker circulation and connect the rising motion over the central-eastern Pacific during a warm phase of ENSO with reduced rainfall over India 17 . While ENSO is one of the significant drivers of the year-to-year variability of the Indian summer monsoon, other forcings are related to SST patterns over the Indian Ocean, Atlantic Ocean, and other parts of the Pacific Ocean affect ISMR. For instance, Indian Ocean Dipole (IOD) is a dominant mode of coupled ocean-atmospheric variability, influencing summer monsoon [18][19][20][21][22] . A positive phase of IOD is linked with the above-normal summer monsoon rainfall [23][24][25] . Moreover, IOD modulates the teleconnection between ENSO and ISMR 25,26 . At the decadal time scale, Pacific decadal oscillation (PDO) influences the summer monsoon through a modification of the Walker circulation and adjustment in the Indian ocean Hadley cell 27,28 . Similarly, at a multidecadal time scale, Atlantic multidecadal oscillation (AMO) influences the summer monsoon 29,30 through modifying the tropospheric temperature gradient 17 .
India has a long history of famines [31][32][33] . Famines caused significant mortality due to lack of food availability, resulting from crop failures and other factors 9 . The history of famines in India reported in the scientific literature goes back to 503-443 BC 31 . Famines occurred in a large number in India during the pre-British era. While the number of local/regional famines that affected a small part of the country is considerably higher 31 , India witnessed twenty major famines from 1200 to 1945. Out of these major famines, nine famines were widespread and caused the loss of millions of lives in India. The summer monsoon failures and droughts 9,31,33 caused a majority of the famines. Agricultural drought due to soil moisture deficit and crop failure resulted in food insecurity, which primarily caused most famines in India during the British Era 9,10 . For instance, out of six major famines (1873-1874, 1876, 1877, 1896-1897, 1899, and 1943) during 1870-2018, soil moisture drought 9 was the main reason for the five (1873- 74, 1876, 1877, 1896-97, and 1899). On the other hand, summer monsoon failure and drought were not associated with the Bengal famine of 1943; instead, it occurred due to policy failures and conflict 9,34 . Therefore, understanding major droughts and the causes of their occurrence in the long-term historical record is crucial.
A relatively short observational record of a century and half period does not provide crucial information on the megadroughts that caused famines before the 1870s in India. Long-term paleoclimate reconstructions based on the historical documents and tree ring datasets [35][36][37][38][39] have been used to examine droughts and famine across the globe [35][36][37][38][39] . Paleoclimate records over the Indian monsoon region showed that the Indian summer monsoon's year-to-year variability is higher than the estimates based on the observations 40 . Seasonally resolved gridded reconstruction of monsoon failures based on the network of tree ring chronologies can reveal the occurrence of megadroughts over the past millennium 35 . As the past droughts were the leading cause of most famines that India witnessed 9,31 , it is crucial to examine the likelihood of the megadroughts and their drivers in the past and future. We identify the two critical research gaps: (i) several famines that occurred before the 1870s and affected millions of people in India; however, their causes and linkage with the monsoon failures are not well understood and (ii) the role of substantial projected warming on the occurrence of megadroughts in India remains unexplored. In consecutive years, a widespread failure of the summer monsoon can cause megadroughts with profound implications for agriculture and water availability. However, estimation of megadrought using a shortterm record is challenging due to its low probability of occurrence. To address these gaps, we use Palmer Drought Severity Index (PDSI) based on the Monsoon Asia Drought Atlas (MADA, 1200-2012) and Climate Research Unit (CRU, 2013-2018) to examine the likelihood of megadroughts that occurred for the 2 and 3 consecutive years during 1200-2018. Also, we use long-term simulations from the Community Earth System Model (CESM) to estimate the likelihood of megadroughts during the summer monsoon season in India in the last millennium (850-2005) and the projected future climate (2021-2100).

Consecutive megadroughts in India, 1200-2018
We used the PDSI for June-August (JJA) period from MADA 35,36 (MADA-PDSI hereafter), which is available at 1°spatial resolution for the monsoon Asia domain (see Methods for more details). MADA-PDSI, available for 1200-2012, was extended till 2018 using gridded observations from Climate Research Units (CRU). Since MADA-PDSI for 1990-2012 is based on the CRU datasets, we used CRU observations to construct PDSI for the 2013-2018 period. We identified consecutive megadroughts during two or three summer monsoon seasons using the long-term (1200-2018) record of PDSI. As consecutive summer monsoon failures have considerable implications for water availability and agriculture 10,41 , megadroughts were identified if these affected more than 40% of the country, and all-India averaged PDSI was less than −1.5 for 2/3 consecutive years. We selected only widespread droughts to reduce uncertainty in MADA-PDSI 35,42,43 , which might be associated with localized droughts. We identified six megadroughts that occurred consecutively for 2 or 3 years during the entire record of 1200-2018 (Supplementary Table 1 Table 1). On the other hand, only two (1431-1433 and 1721-1723) megadroughts remained consecutively for 3 years. The mean area for five out of six megadroughts was more than 50%, while the maximum area affected by megadrought ranged from 48 to 72% of the country (Supplementary Table 1). We ranked these megadroughts based on the overall severity score that combines intensity, duration, and area under drought 7 . Based on the overall severity score, 1721-1723 was the most severe drought during the monsoon season in India. The peak intensity of the megadroughts occurred in 1721, 1431, 2002, 1257, 1424, and 1703 with the affected area of 72%, 52.5%, 65%, 51.5%, 52%, and 48% of the country, respectively (Fig. 1). 2002-2003 was the only megadrought that occurred during the recent decades, which affected more than 65% of the country (Fig. 1). The 2002 drought was primarily caused by a rainfall deficit (~51%) in the month of July 44,45 .

Reconstruction of famines caused by droughts
The spatiotemporal information from the MADA on the summer monsoon failures provides a long-term context of droughts that occurred during the past over Monsoon Asia 46 . Famines in India were instigated by droughts and caused the death of millions of people 9 . Instrumental meteorological record in India does not go beyond 1870 9,10 , which hampers our understanding of the underlined causes of the famines before 1870. Nine major famines caused by droughts were identified from the database developed by Walford (1878) 31 and other sources 47 48 . As droughts and the summer monsoon failures 9 caused most famines in the post-1870 period, we used MADA-PDSI to identify the droughts and affected regions during each famine (Fig. 2).
India experienced a famine that was extended to the entire country during 1344-1345 31,49 . The famine of 1344-1345 severely affected the Deccan part of India 31 . The famine was so acute that the Mughal emperor could not procure necessary items for his household 31,49 . MADA-PDSI captured the 1345 drought, which affected about 43% of the country (Supplementary Table 2, Fig. 2). The subsequent major famine occurred in 1631. Walford (1878) 31 reported that a famine was caused by drought and war in 1631. Moreover, the great drought that occurred in 1631 affected a large part of India and Asia 31 . The 1630-1632 famine is also known as the Deccan famine, which affected Deccan Plateau, Khandesh, and Gujarat 50 . The famine resulted from staple crops' failure for 3 consecutive years and led to intense hunger, diseases, and migration. 1630-1632 was one of India's worst famines during the Mughal Empire 51,52 . We find that MADA-PDSI captured the drought during 1630-1632 peaked in 1632 (Fig. 2). The 1630-1632 drought affected about half of the country and was centered in India's central and northwestern parts (Fig. 2). The other major famine in India occurred during 1702-1704 31 . The famine affected western India's parts and was mainly centered in Thar and Parkar (now in Pakistan) districts 31 . The famine was caused by drought killed more than two million people 51 . We find that 1702-1703 is also captured using MADA-PDSI as the megadrought of the 2 consecutive years (Fig. 2 India experienced another great famine during 1769-1773 31,53 , also known as the Great Bengal Famine of 1770. The total reported deaths due to the famine ranged between 2 and 10 million 53 , with an approximate reduction of 7-33% in the total population of Bengal 54 . Based on the MADA-PDSI, we find that the 1769-1771 drought severely affected northern India and the Gangetic Plain, which peaked in 1771, affecting about half of the country (Supplementary Table 2; Fig. 2). The subsequent major famine in India occurred during 1802-1804 31 . The drought of 1802-1804 peaked in 1804 and affected about 55% of the country, which was well captured by MADA-PDSI (Fig. 2, Supplementary Table 2). India witnessed another notable famine in 1837-1838, which affected northwest provinces 31 . The 1837-38 famine is also known as the Agra famine of 1837-1838, which involved a large part of Uttar Pradesh. We find that the MADA-PDSI captured the drought during 1837, which affected about 35% of the country (Fig. 2). India witnessed a major famine in 1866, which is also known as the Orissa famine of 1866. The famine affected Bengal, Orissa, and Bihar's lower provinces, which caused the death of one and half million people 31 . The famine was caused by the failure of the summer monsoon rainfall during the late season 48 . The famine had an unprecedented influence on the demography as about one-third of the total population of Orissa died 55 . MADA-PDSI captured the drought reasonably well (Fig. 2, Supplementary Table 2). Madras's Presidency witnessed one of the most extended famines in 1877 that caused the death of more than 10 million people. MADA-PDSI captured the drought in Bombay (now in Maharashtra), but the drought intensity is weaker in south India than reported in the previous studies 9, 10 . As MADA-PDSI covers only JJA period, it does not account for the summer monsoon variability during the late monsoon season (September). Also, MADA-PDSI does not account for the rainfall variability during the northeast monsoon (October-December), affecting southern India 56 . The last famine of the 19th century occurred in 1899-1900, which affected central and western parts of India. The famine caused deaths of about 1-4.5 million people and affected central India, Hyderabad, and Rajputana agency 47 . MADA-PDSI captured drought in some famine-affected regions (Fig. 2), which covered about 40% of the country in 1900 (Supplementary Table 2, Fig. 2). Overall, the famines caused by the summer monsoon droughts in the previously reported studies are captured by MADA-PDSI, which provides us with a basis to examine the occurrence and likelihood of the consecutive megadroughts in India.
The leading mode of variability of droughts Next, we compared MADA-PDSI based on paleoclimate records against the PDSI estimated using gridded observations from CRU for a common period of 1901-1989 ( Supplementary Fig. 1). PDSI from MADA v2 for the 1200-1989 period is based on the paleoclimate reconstruction 35,36 while the estimates after 1989 are based on the CRU observations 57 . The reconstructed PDSI from the paleoclimate records compares well (correlation = 0.63, pvalue < 0.05) with the CRU-based PDSI ( Supplementary Fig. 1). All the major droughts during the 1901-1989 period were well captured by the paleoclimate-based reconstructed PDSI (Supplementary Fig. 1b). Comparison of PDSI based on paleoclimate record and CRU data further demonstrates that reconstruction based on a network of tree-ring chronologies provides valuable information on droughts during the summer monsoon over India 35 . Identifying droughts that caused famines provides confidence to use the MADA-PDSI to examine the megadroughts' variability and occurrence in India.
We used empirical orthogonal function (EOF) analysis to examine the dominant mode of variability in MADA-PDSI for the entire record of 1200-2018 (Fig. 3). The principal component of the leading mode (PC-1) of MADA-PDSI is well correlated (r = 0.78, p-value < 0.05) with PDSI. Also, all the six megadroughts were well captured by the PC-1 (Fig. 3a). The leading mode of variability in MADA-PDSI explains about 22.7% of the total variance (Fig. 3b). The spatial pattern of EOF-1 of MADA-PDSI shows positive loadings in most of India except for the southern peninsula, northeastern regions, and Jammu and Kashmir (Fig. 3b). Since the negative values of PC-1 co-occur with the megadroughts (PDSI), leading EOF mode (EOF-1) of MADA-PDSI indicates that most of the megadroughts were centered in the central and northern parts of the country (Fig. 3b). We constructed the composite of MADA-PDSI for the six megadroughts (Supplementary Table 1), which exhibits a similar spatial pattern obtained from EOF-1 (Fig.  3c). We also identified all the droughts during the summer monsoon (JJA) with PDSI less than −1.5 and covered 40% or more of the country during 1200-2018. A total of fifty such droughts occurred in the record of 819 years (1200-2018). MADA-PDSI composite of all the 50 droughts displayed a somewhat identical spatial pattern that was reflected in EOF-1 and the 6 megadroughts. Therefore, the leading mode of drought variability obtained from the MADA PDSI shows that drought-affected northern and central India more than the other parts of the country (Fig. 3).
We used maximum covariance analysis (MCA 58,59 ) to obtain the leading mode of coupled variability between the summer monsoon precipitation and SST anomalies (Fig. 4). SST anomalies over the Pacific Ocean are negatively associated with the summer (JJA) monsoon precipitation over India during the observed record of 1901-2018, for which observed SST is available (Fig. 4a, b). The leading mode of coupled variability obtained from the MCA explains about 70% of total squared covariance and corresponding principal components (PC-1) are strongly correlated (r = 0.53, p-value < 0.05). The long-term record (1200-2018) of MADA-PDSI showed that the megadrought of 2002-2003 was the most severe  Table S1, Fig. 4c). We, therefore, examined the spatial variability in MADA-PDSI during the summer monsoon (JJA) of 2002 and associated SST anomaly (Fig. 4d, Supplementary Fig. 2). A majority of north India was affected by drought during 2002 (Fig. 4c). For instance, drought was severe over the Indo-Gangetic Plain and northwestern India. Based on the Hadley Center's (HadSST) observations, SST exhibits warmer anomaly over the Pacific and Indian Ocean (Fig.  4d, Supplementary Fig. 2). Warmer SST anomaly over the Pacific during the summer monsoon of 2002 indicates the presence of El Nino 4,44 . El Nino and warmer SST over the central equatorial Pacific are more strongly related to India's monsoon failures than the warmer SST anomaly over the eastern equatorial Pacific 60 . In addition to El-Nino, IOD also influences the summer monsoon rainfall, a dominant mode of coupled variability between the ocean and atmosphere [18][19][20][21][22] . The positive dipole mode index (DMI) is linked with above normal ISMR [23][24][25] and modulates the teleconnection between El Nino and summer monsoon rainfall 25,26 . For instance, concurrent El Nino and positive DMI reduces the effect of El Nino on the Indian summer monsoon and strengthen the monsoon rainfall 25 . We find a negative DMI and warm SST anomaly over the Indian Ocean during the summer monsoon season of 2002 (Fig. 4d, Supplementary Fig. 2). El Nino and negative DMI strengthen the intensity of drought during the year 2002. Moreover, warmer SST anomalies in the Indian Ocean are linked with the drying over the Indo-Gangetic Plain 13,61 . Land-sea thermal gradient has been declining over South Asia due to higher warming over the Indian Ocean than lesser warming over the land 61 . Due to a decrease in land-ocean thermal gradient, Hadley circulation weakens, which results in reduced rainfall over parts of northern India 61 . Moreover, warming in the Indian Ocean SST is linked with the troposphere ridge over north India, resulting in a decline in rainfall 62 . PDO and AMO also influence summer monsoon rainfall over India directly or indirectly 27,63 . We find a negative phase of PDO during the onset of the summer monsoon drought of 2002, which turned to a positive phase during late 2002 and 2003. Droughts during the summer monsoon season can intensity during the El-Nino overlapped with a positive phase PDO 27 . We find a negative phase of AMO during the summer monsoon of 2002, which suppresses rainfall over India 29,30 . Therefore, SST warming in both the Pacific and Indian Oceans affected the 2002 drought in India during the summer monsoon.
We obtained total column rainwater (TCRW) and wind (u and v components) from ERA5 reanalysis to examine the climatological and anomalous features during the 2002 observed drought (Supplementary Fig. 3). Climatological mean patterns show high TCRW over India and the central Pacific Ocean during the summer monsoon. A strong westerly flow signifies the moisture transport from tropospheric winds from the Indian Ocean to the Indian landmass ( Supplementary Fig. 3a). During the summer monsoon of 2002, the weakening of westerly winds can be noted over the Indian mainland along with the dry anomaly of TCRW (Supplementary Fig. 3b). Therefore, the decline of westerly moisture transport is associated with the 2002 drought over India, as reported in the previous studies 44 . The 2002 drought was caused by a rainfall deficit of 21.5% during the monsoon, while 56% rainfall deficit occurred during July 64 . Consistent with the SST anomaly patterns over the Pacific and Indian oceans during the summer monsoon of 2002, a negative relationship between MADA-PDSI and SST anomalies (JJA) exhibited the role of ENSO on the summer monsoon failure and droughts over India 13,61 (Fig. 4a,  b). For instance, MADA-PDSI is significantly (P-value < 0.05) correlated with the SST anomalies over the Pacific and Indian Oceans (Fig. 4b). Also, Overall, the long-term MADA-PDSI and observed SST exhibited the role of ENSO on droughts over India as reported in the previous studies 13,60,61 .

Consecutive megadroughts in the past and future
The occurrence of the consecutive megadroughts during the summer monsoon season over India has profound implications 1,7,8 . For instance, consecutive monsoon failures over the Indo-Gangetic Plain during 2014-2015 caused water scarcity and groundwater depletion 8 . Moreover, the drought during 2016-2018 in South India due to the northeastern monsoon's failure resulted in a water crisis over a large region 36 . As CESM-Last Millennium Ensemble (LME) [CESM-LME hereafter] captures the coupled variability between summer monsoon rainfall and SST reasonably well 65 (Supplementary Fig. 4), we estimated the occurrence of the consecutive megadroughts using the 12 ensemble members from the CESM-LME. Also, CESM-LME captures the critical feature of the Indian summer monsoon and seasonal cycle of rainfall well ( Supplementary Fig. 5). First, we estimated PDSI for each ensemble member of the CESM-LME, drought characteristics were calculated. Out of the four megadroughts that occurred consecutively for 2 years based on the MADA-PDSI (Supplementary Table 1), the drought of 2002-2003 had the highest overall severity score (274.6), while the drought of 1702-1703 had the lowest overall severity score (160.6). We considered the overall severity score of the 1702-1703 drought as a benchmark (with the lowest overall severity score) to estimate the likelihood (number of megadroughts/total number of years) CESM-LME during 850-2005. The likelihood was estimated using the number of megadroughts divided by the total number of years. We first estimated the number of megadroughts in all the ensemble members of CESM-LME (12*1156), divided by the number of years. The likelihood of such (~1702-1703) megadroughts (affected area more than 40% and PDSI less than −1.5) during the summer monsoon (JJA) in two consecutive years is less than one (~0.73) event in 100 years. On the other hand, the likelihood of the megadroughts like 2002-2003 is even lesser (~0.34 events in 100 years). Long-term PDSI from the CESM-LME for 1156 years (850-2005) from each ensemble member enables us to estimate the likelihood of megadroughts with a low probability of occurrence.
The composite MCA patterns of coupled variability between the summer monsoon precipitation and SST from all 12 ensemble members of CESM-LME show a negative relationship between summer monsoon rainfall and SST over the Pacific Ocean (Fig. 5a,  b). The relationship exhibits clear patterns associated with the ENSO signature in CESM-LME simulations during the weak summer monsoon, consistent with the observed period (Fig. 4). To further confirm the role of SST variability over the Pacific Ocean, we constructed atmospheric and oceanic anomaly composites for all the megadroughts that occurred during the two consecutive years in CESM-LME (Fig. 5c, Supplementary Fig.  6). A total of 101 megadroughts of 2 successive years were identified in all the 12 ensemble members of CESM-LME. Consistent with the observed megadroughts (Fig. 4), the positive SST anomaly over the Pacific Ocean is strongly associated with the consecutive 2-year megadrought during the summer monsoon season in CESM-LME (Fig. 5d). A majority (~60%) of megadroughts (61 out of 101) were associated with El Nino, positive IOD, and negative phase of PDO during the last millennium based on the CESM-LME simulations (Fig. 6a), which is consistent with the SST anomalies present during the observed megadrought of 2002-2003. Out of 61 megadroughts that were associated with El Nino during the last millennium, around 70% were also linked with positive DMI, negative PDO, and positive AMO (Fig. 6c). Analysis of climatological mean total precipitable water (TMQ) and wind field at 850 hPa showed the presence of higher TMQ over the tropical region and strong westerly wind flow during the summer monsoon (JJA) [ Supplementary Fig. 6a]. Anomaly composite of TMQ and tropospheric winds showed dry TMQ anomalies and weaker westerly winds during the megadroughts of the two consecutive years ( Supplementary Fig. 6b). The anomalous anticyclonic pattern observed for the megadroughts in CESM-LME is consistent with the observational analysis based on ERA5 reanalysis (Supplementary Figs. 5 and 6). Similar climatological and anomalous patterns of TMQ and wind at 850 hPa were obtained for the megadroughts that occurred for the 3 consecutive years (Supplementary Fig. 6c). Megadroughts were widespread, with a mean coverage area of about 72%. Only two (1431-1433 and 1721-1723) megadroughts based on MADA-PDSI occurred for the 3 consecutive years (Supplementary Table  1). We considered the overall severity score (~311) of the 1431-1433 drought as a benchmark (the lowest severity score) to estimate the likelihood using CESM-LME. The probability of the three consecutive years megadroughts during the summer monsoon season is about 0.3 events per 100 years based on CESM-LME. Consistent with the megadroughts of two successive years, the role of positive SST anomalies on 3 consecutive years megadroughts over India was also exhibited ( Supplementary  Fig. 6c).
Like CESM-LME, CESM-Large Ensemble (CESM-LE hereafter) also captures the summer monsoon rainfall variability reasonably well during the historical climate 59,66 . We used CESM-LE simulations from 40 different ensemble members for the 2021-2100 period to estimate the likelihood of megadroughts in the future climate. We calculated PDSI for each ensemble member of CESM-LE to identify megadroughts and their characteristics (Fig. 7a, b). The likelihood  Supplementary Fig. 7), which is consistent with observations and CESM-LME simulations. We find that the future megadroughts during the summer monsoon over India will also be associated with El Nino, positive IOD, and negative phase of PDO (Fig. 6c, d). Moreover, the climatological and anomaly patterns of total precipitable water and wind at 850 hPa showed dry anomalies for TMQ and weaker westerly flow during the megadroughts that are projected to occur in the future (Supplementary Fig. 7). However, we note that the TMQ anomalies are much stronger for the projected future climate than for the last millennium, attributable to the rising trend in total precipitable water under the warming climate ( Supplementary  Fig. 7). The likelihood of the megadroughts during the summer monsoon declines in the future (Fig. 7), which can be attributed to an increase in the summer monsoon precipitation under the warming climate 66 . Moreover, MADA-PDSI is based on the June-August period, which does not consider the summer monsoon's variability during the late season. For instance, a considerable decline in India's summer monsoon season precipitation is reported during the late monsoon season in CESM-LE 59 , affecting the severity and frequency of droughts in India.

Influence of climate warming on megadroughts
Droughts during the summer monsoon are mainly driven by the precipitation deficit, which can be exacerbated by warming 7,66 .
We examined the fraction of droughts caused by precipitation deficit, warming, and both by precipitation deficit and warming using the CESM-LME and CESM-LE simulations during the past and future climate, respectively (Fig. 8). We estimated the precipitation and temperature anomalies for all the megadroughts that occurred for two consecutive years for each ensemble member of CESM-LME and CESM-LE to identify the role of temperature in the droughts in the past and future. Droughts with precipitation deficits (negative precipitation anomaly) and negative temperature anomalies were considered precipitation deficit droughts (Fig.  8). During the summer monsoon and positive temperature anomalies, drought with precipitation deficits were considered precipitation deficits and temperature-driven droughts. On the other hand, droughts with no precipitation deficit and positive temperature anomaly were identified as temperature-driven droughts. Since precipitation and temperature anomalies during the summer monsoon are negatively correlated 66 , precipitation deficit can cause warming due to land-atmospheric feedback 59,65,67,68 . However, we did not quantify the role of land-atmospheric coupling on warming and its impact on drought. Therefore, precipitation deficit and temperature-driven droughts may be more dominated by precipitation anomalies during the historical climate. Anomalies of precipitation and temperature during the summer monsoon season were estimated against the reference of 1951-2000 for both CESM-LME and CESM-LE simulations. In the last millennium (CESM-LME), the megadroughts were driven by either only precipitation deficit or a combination of precipitation deficit and warm temperature (Fig. 8a). For instance, a majority (~78%) of megadroughts during the last millennium were caused by both precipitation deficit and warming during the summer monsoon. Around 22% of the total megadroughts in the last millennium were caused by precipitation deficit only during the summer monsoon season (Fig. 8a). However, megadroughts in the future are likely to be caused either by precipitation deficit and warming or solely due to a substantial rise in temperature (Fig. 8b). For instance, about 55% of future megadroughts are projected to be driven by both precipitation deficit and warming, while 45% of future megadroughts is driven only by due to rising in temperature (Fig. 8b). The substantial rise in air temperature and near-normal summer monsoon precipitation can also cause drought due to higher evapotranspiration 69 and increased atmospheric water demands. Moreover, only precipitation deficit megadroughts and a combination of precipitation deficit and warm temperature megadroughts showed similar drought severity and intensity during the last millennium (Fig. 8c, d). However, precipitation deficit and warm temperature-driven megadroughts are projected to be more severe and intense than temperaturedriven megadroughts under the projected future climate (Fig. 8e,  f). Overall, we find that the likelihood of megadroughts consecutively for two years during the monsoon season is relatively low in future climate, which can be influenced by the significant warming that is projected. Therefore, the occurrence of megadroughts cannot be neglected under the warming climate.

DISCUSSION
Summer monsoon rainfall is critical for the food and water security of the large population in India. Understanding the drivers and likelihood of megadroughts is essential for the future food and water security of one of the world's most populated regions. However, the short instrumental record of the Indian summer monsoon underestimates the magnitude of rainfall variability 37 .
Recently gridded rainfall and temperature observations were extended to 1870, which were used to reconstruct the soil moisture droughts during the famines of the British Era 9 . We used long-term paleoclimate reconstructions for monsoon Asia, simulations from the CESM ensembles for the last millennium, and projected future climate to examine drivers of the megadroughts. As MADA-PDSI is based on paleoclimate reconstructions, we considered the widespread and severe droughts in our analysis. We found that the reported famines due to monsoon failures were well captured in MADA-PDSI. Failure of the summer monsoon for consecutive years has detrimental implications for water availability, groundwater recharge, and food security. Recent droughts that occurred for the two successive years (2014)(2015) in north India and South India (2016-2018) had lasting implications 8,56 . We used a longterm PDSI based on MADA and CRU for the 1200-2018 period to identify the megadroughts that occurred for the two and three consecutive years. The likelihood of megadroughts for 2 and 3 successive years during the monsoon season is about 0.5 and 0.25 events per 100 years based on MADA-PDSI for 1200-2018. However, the likelihood of megadroughts was slightly increased in the longer-term record from the CESM-LME simulations. For instance, we find that the likelihood of 2-and 3-years megadroughts are about 0.73 and 0.3 events per 100 years in CESM-LME. The likelihood of 2-and 3-year megadroughts is projected to decline in the future climate, attributed to an increase in summer monsoon precipitation 66,70 .
Summer monsoon failure and droughts played a vital role during the six significant famines captured using the MADA-PDSI. While most droughts in the past were caused due to the combination of precipitation deficit and warming, substantial warming that is projected in the future can cause droughts with near-normal precipitation during the summer monsoon. Given the dependence of a large population on agriculture and depletion of groundwater in several parts of the country 71,72 , the occurrence of megadroughts in the future can have large implications despite its lower likelihood. While the megadroughts in the past caused food insecurity and famines in India 9 , famines are unlikely in the future despite the occurrence of megadroughts. India has made tremendous progress in improving crop yields due to enhanced agricultural practices and irrigation. Moreover, there has been a considerable growth in transportation and public distribution systems in the post-independence in India. For instance, despite the several intense and widespread droughts after the independence in India, the country did not experience famines 9 . The projected warming can result in a decline in the crop yields 73,74 , which the hot and dry extremes 66 can further worsen. However, the risk of famines in the future remains low due to flood storage, irrigation, public distribution system, rural employment, and transportation network 75,76 .
Based on our findings, we conclude the following: • Based on CESM-LME simulations, the likelihood of 2 and 3-year consecutive megadroughts during the summer monsoon is about 0.7 and 0.3 events per 100 years, respectively. A majority of the past and future megadroughts are linked with El Nino. The frequency of megadroughts in the future is projected to decline due to increased monsoon season precipitation. However, substantial warming that is projected can influence the future megadroughts in India, which will have implications for food and water security.

Observed data (1200-2018)
We used the combined long-term Monsoon Asia Drought Atlas (MADA 35 (Dai-PDSI) was developed using the precipitation and temperature from the CRU. Therefore, we used the PDSI from the CRU (https://crudata. uea.ac.uk/cru/data/drought/) at 0.5°spatial resolution for 1901-2018 80,81 , which was further regridded at 1°spatial resolution to validate the MADA-PDSI data. Using PDSI (MADA-PDSI from 1200-2012 and CRU-PDSI from 2013-2018), we evaluated the summer monsoon drought severity and occurrence over India for 1200-2018. We used observed precipitation and sea surface temperature (SST) for the period 1901-2018. The monthly gridded precipitation was obtained from CRU (version 4.03 57 ) for 1901-2018 at 0.5°spatial resolution over India. CRU-precipitation is developed using many observed datasets using the angular-distance weighting interpolation at a monthly time scale 57 . Gridded products from the CRU have been widely used for regional and global drought and hydroclimatic extreme studies 82,83 .
We obtained monthly SST from the Hadley Centre's (HadSST, https:// www.metoffice.gov.uk/hadobs/hadisst/) from 1901 to 2018 to evaluate the linkage between drought (PDSI and precipitation anomalies) and SST 84 . The HadSST data were reconstructed using the two-stage reduced-space optimal interpolation of observation data at 1°spatial resolution 84 . We estimated oceanic indices, including the ENSO, DMI, AMO, and PDO using HadSST for the observed period. ENSO was estimated using the areaweighted average SST anomalies over the Niño 3.4 region (5°N-5°S, 120-170°W; Supplementary Fig. 8). Indian Summer Monsoon droughts are highly correlated with the SST anomalies over the Nino-3.4 region compared to the other Nino regions 7 . DMI was calculated using the difference between the area-averaged SST anomalies of western (10°S-10°N , 50°E-70°E; Supplementary Fig. 8) and eastern (10°S-0°N, 90°E-110°E) Indian Ocean regions 85 . Similarly, area-averaged SST anomaly in the North V. Mishra and S. Aadhar Atlantic region (0°N-60°N, 0-80°W; Supplementary Fig. 8) represents the AMO 86,87 while PDO was estimated using the leading pattern of SST anomalies over Northern Pacific Ocean (north of 20°N) 27,88 ]. We used EOF analysis to identify the leading principle component over the Northern Pacific Ocean (20°N-60°N, 120°E-100°W; Supplementary Fig. 8), which was used for estimating PDO. Moreover, we obtained TCRW, and wind (u and v components) at 850hPa from the European Centre for Medium-Range Weather Forecasts Reanalysis version 5 (ERA-5 89 ) for the period 1979-2018 to evaluate available moisture content and atmospheric condition during observed megadroughts.

Data from CESM
Community Earth System Model realistically simulates the coupling between ISMR and SST variability over the Pacific Ocean 66 . Moreover, CESM-1 simulations reproduce seasonality and spatial variability in summer monsoon rainfall reasonably well 46 . Most global climate models show bias in summer monsoon precipitation and do not capture the ENSO-monsoon coupling well 70,90 . We, therefore, selected the CESM-1 for the present work. Moreover, simulations from many ensemble members are available for the past, present, and future periods, making CESM suitable for the analysis to estimate the likelihood of megadroughts in India. We used simulations from the CESM-LME 91 project (https://www. cesm.ucar.edu/projects/community-projects/LME/) to evaluate the likelihood of the occurrence of the consecutive droughts over India during the period 850-2005. We used 12 runs (out of 13; the first run does not have the data for the entire period of 850-2005) of CESM-LME for the 850-2005 period to estimate monthly PDSI. Monthly precipitation, surface air temperature (maximum and minimum), specific humidity, wind speed, latent heat flux, sensible heat flux, and surface pressure were used to estimate India's PDSI. PDSI was calculated using precipitation, atmospheric water demand (potential evapotranspiration: PET), and available water content (AWC) of soil 92 . Atmospheric water demand (PET) was estimated using the Penman-Monteith method 93,94 , while AWC was obtained from the Harmonized World Soil Database (HWSD 95 ). We analyzed the consecutive droughts and their characteristics over India in the last millennium using PDSI from the 12 ensemble members from CESM-LME.
We used data from the 40 ensemble members from the Community Earth System Model-Large Ensemble (CESM-LE 96 ) project (http://www. cesm.ucar.edu/projects/community-projects/LENS/) for the period 2021-2100 to examine the drought occurrence in the future climate. We estimated the atmospheric water demand (PET) using the modified Penman-Monteith method 97 for the future (2021-2100) using monthly surface air temperature, wind speed, specific humidity, surface pressure, and latent and sensible heat flux. We used the modified Penman-Monteith method since the changes in surface resistance due to the increased warming climate were not considered in the Penman-Monteith method, affecting the drought severity in the future climate 69 . Using the PET, precipitation, and AWC, we estimated PDSI for 2021-2100 using the 40 runs of CESM-LE to understand the future drought severity. The CESM-LE simulations have been widely used to analyze drought and extreme climate studies under the future climate 66,69,98 . We used monthly SST, total precipitable water (TMQ), and wind at 850hPa from CESM-LME and CESM-LE for 850-2005 and 2021-2100, respectively. As SST and TMQ from CESM-LE for the future period (2021-2100) has a strong warming trend, we used Ensemble Empirical Mode Decomposition (EEMD 99 ) to remove secular trend 100 from TMQ and SST for the period 2021-2100. Moreover, we examined the year-to-year variability between the summer monsoon (JJA) precipitation and SST 13 . The MCA 58 method was used to evaluate the coupled variability between precipitation and SST. We identify the first leading mode of variation between rainfall and temperature using the MCA, which shows the substantial coupled variability between two fields (precipitation and SST). Similar to observation, we estimated oceanic indices including ENSO (Nino 3.4), DMI, AMO, and PDO from CESM-LME and CESM-LE simulations for the past and future climate. Moreover, TMQ and wind at 850hPa were used to evaluate available moisture content and atmospheric condition during megadroughts in LENS-LME and LENS-LE.

DATA AVAILABILITY
All the datasets used in this study can be obtained from the corresponding author.

CODE AVAILABILITY
Data processing scripts can be obtained from the corresponding author.