Droughts and Mega-Droughts

ABSTRACT

Drought is the deadliest natural disaster on Earth due to its long duration, wide spatial coverage and direct connection with the food supply for human beings. Drought develops slowly and thus is called a silent killer. This paper reviews the history of research on droughts and mega-droughts. Interannual droughts are in many places driven by the El Nino-Southern Oscillation, multi-decadal mega-droughts are often driven by the Atlantic Multi-Decadal Oscillation, while centennial mega-droughts are often driven by the Global Inter-Centennial Oscillation. More generally, droughts are affected by multiple factors including global sea surface temperature anomalies, local land-atmosphere feedbacks, internal atmospheric variability and external forcings from outside the climate system. Possible future research directions are also suggested.

RÉSUMÉ

[Traduit par la redaction] La sécheresse est la catastrophe naturelle la plus meurtrière sur Terre en raison de sa longue durée, de sa large couverture spatiale et de son lien direct avec l’approvisionnement alimentaire des êtres humains. La sécheresse se développe lentement et est donc appelée un tueur silencieux. Le présent article porte sur l’histoire de la recherche sur les sécheresses et les méga-sécheresses. Les sécheresses interannuelles sont dans de nombreux endroits dirigées par l’oscillation El Niño-oscillation australe, les méga-sécheresses multi-décennales sont souvent dirigées par l’oscillation atlantique multi-décennale, tandis que les méga-sécheresses centennales sont souvent dirigées par l’oscillation globale inter-centennale. Plus généralement, les sécheresses sont affectées par de multiples facteurs, notamment les anomalies de la température de surface de la mer, les rétroactions terre-atmosphère locales, la variabilité atmosphérique interne et les forçages externes au système climatique. Des directions de recherche futures possibles sont aussi suggérées.

Keywords:

• drought
• mega-drought
• El Nino-Southern Oscillation
• Atlantic Multi-Decadal Oscillation

1 Introduction

Drought, sometimes called the silent killer, is probably the deadliest natural disaster on Earth. Research shows that droughts have strongly affected human history (Fig. 1). A series of megadroughts between 135,000 and 75,000 years ago may have been responsible for the first migrations of early humans out of Africa, who were the ancestors of all non-Africans today around the world (Scholz et al., 2007). A megadrought hitting the Middle East and parts of Europe 4,500 years ago was a significant reason for the collapse of the Old Kingdom in Egypt around 4200 cal BP (Barta & Bezdek, 2008; Marriner et al., 2012). A century-long mega-drought in Europe between 200 and 300 AD preceded the collapse of Western Roman Empire during the Migration Period (Fig. 1a; Buntgen et al., 2011). A halving of annual rainfall in Central America 1200 years ago caused a war with neighbouring nations over dwindling water resources leading to the demise of the Mayan civilization (Evans et al., 2018; Kennett et al., 2012). Long-term drought records in China correspond well with the collapse of the Tang, Yuan, and Ming Dynasties, all times that were characterized by famine and popular unrest (Fig. 1b; Zhang et al., 2008). The deadliest drought in human history was the 1875–1878 Great Drought and resultant Global Famine, which caused a death toll of more than 50 million people in China, India, Indonesia, Australia, Brazil, Colombia, Venezuela and South Africa (Davis, 2001; Singh et al., 2018).

Fig. 1 (A) Reconstructed European AMJ precipitation totals (top) and JJA temperature anomalies (bottom) with respect to the 1901–2000 period. Bold lines are 60-year low pass filtered data (from Buntgen et al., 2011). (B) The WX42B d18O record from China. Chinese dynasties are indicated (from Zhang et al., 2008).

World population has almost doubled in the past four decades from 4 billion in 1975–7.9 billion in 2020 (World Bank, 2021), and the largest growth happened in the poorest developing countries in Africa and Central Asia (see Fig. 2 for the global map of population density), making those regions more vulnerable to droughts and associated failure of agriculture. The deadliest and costliest droughts in the past 120 years (1900–2020) are shown in Fig. 2a. The 1899–1900 drought and resultant famine affected 59.5 million people with mortality estimates of 1–4.5 million (Mishra et al., 2019). The 1984 drought in Africa affected 250 million people in 22 countries and caused a death toll of one million people (Druyan, 1989). The 1988 North American drought caused $60 billion in damage (E. R. Cook et al., 2007; about $120 billion in today’s dollar), which was similar to the damage of Hurricane Katrina and wiped out more than 25% of the U.S. GDP growth in that year. Figure 2b shows the normalized precipitation anomalies for global continent mean from the Global Precipitation Climatology Centre (GPCC; Schneider et al., 2016) dataset, global mean from the Global Precipitation Climatology Project (GPCP; Adler, 2003) dataset, Sahel rainfall from Nicholson et al. (2012), Indian monsoon rainfall from the Indian Institute of Tropical Meteorology (Mooley et al., 1981), north China mean from the GPCC dataset, and central/western U.S. mean from the GPCC dataset. The precipitation anomalies are smoothed with 11-year running mean to display the low-frequency variability at decadal and longer timescales. The decadal precipitation anomalies in Sahel, India, China and the global continent mean show coherent variations. The decadal precipitation anomalies in central/western U.S. tend to be in opposite sign with the other regions, which has been multiplied by −1 in Fig. 2b. The deadliest and costliest droughts often occurred in the multi-decadal dry periods. This coherent variation of precipitation anomalies around the globe has been speculated to be driven by the coherent variation of SST anomalies in different ocean basins (Fig. 2c; Folland et al., 1986; Sutton & Hodson, 2005; Zhang & Delworth, 2006). However, the modern instrument records are too short, and long-term paleoclimate datasets are needed to derive statistically robust relationships.

Fig. 2 (A) The deadliest and costliest droughts in the past 120 years (1900–2020). Shading shows the population density in 2015. Damage has been adjusted to 2020 US dollars. (B) Global and regional mean precipitation anomalies (see text for details). Note the central/western U.S. anomalies have been multiplied by −1. (C) Global and basin-averaged SST anomalies from the ERSST dataset (detrended). The anomalies are smoothed with 11-year running mean to display the low-frequency variability at decadal and longer timescales.

Drought can be defined as the shortage of water. Since water exists in all components of the global water cycle including the atmosphere, land, ocean and sea ice, shortage of water can happen in any of the components due to repartitioning of the total amount of water among the different components, or due to decrease of the total amount of water on Earth. Drought studies generally focus on the shortage of water over land, which directly affects human society (see reviews by Mishra & Singh, 2010; Peters-Lidard et al., 2018; Schubert et al., 2016; Wilhite, 2000, 2016). The land surface water budget has five key components: precipitation, evapotranspiration, runoff (also called streamflow), infiltration and soil moisture storage (Fig. 3a). There are four types of droughts generally considered in drought research: meteorological, hydrological, agricultural, and socioeconomic (Eslamian et al., 2017; Keyantash & Dracup, 2002; Wilhite & Glantz, 1985). Meteorological drought results from a shortage of precipitation, while hydrological drought describes a deficiency in the volume of the water supply at the ground, which includes streamflow, reservoir storage, and/or groundwater heights. Agricultural drought relates to a shortage of available water for plant growth, and is assessed as insufficient soil moisture to replace evapotranspirative losses. The socioeconomic drought is a consequence of the other drought types when societal demand exceeds natural supply. The four types of droughts are connected to one another. Generally, meteorological drought propagates to hydrological drought, then to agricultural drought, and finally to socioeconomic drought. When meteorological drought propagates to hydrological drought, other terms of the surface water budget also contribute. When hydrological drought propagates to agricultural drought, human activities such as irrigation and crop management also contribute. When agricultural drought propagates to socioeconomic drought, political and economic conditions also contribute, and it was often the case in history that corrupted governments badly enhanced the drought disasters. Scientifically, all the four types of droughts are ultimately driven by shortage of precipitation, which is the leading forcing of the surface water budget (Qian et al., 2006, 2007), although evapotranspiration change can also contribute significantly in some regions (B. I. Cook, et al., 2014; B. I. Cook, Ault et al., 2015; Schubert et al., 2016; Teuling et al., 2013). In the following we will focus on the meteorological droughts.

Fig. 3 (A) Schematic depiction of land-surface water budget and physical mechanisms for droughts and mega-droughts. (B) Global map showing the season with the largest climatological precipitation at each grid from the GPCC dataset. Magenta lines show the monsoon regions with summer-winter precipitation difference larger than 300 mm.

Many drought indices have been developed for the meteorological droughts (Eslamian et al., 2017; Heim, 2002; Keyantash & Dracup, 2002; Mishra & Singh, 2010; Oladipo, 1985), such as the discrete and cumulative precipitation anomalies, the rainfall deciles, the rainfall anomaly index (van Rooy, 1965), the Drought Area Index (Bhalme & Mooley, 1980), and the Standardized Precipitation Index (McKee et al., 1993). The Palmer Drought Severity Index (PDSI; Palmer, 1965), which has been used by many studies on climate change and paleoclimate droughts, was considered as a meteorological drought index in earlier reviews (e.g. Keyantash & Dracup, 2002; Mishra & Singh, 2010; Zargar et al., 2011), but as an agricultural drought index in later reviews (e.g. Eslamian et al., 2017; Yu et al., 2019). Palmer (1965) used the surface precipitation and the surface air temperature to drive an idealized two-layer soil model, and calculated other terms of the land surface water budget. The PDSI has well-known limitations, such as region-dependent empirical parameters preventing spatial comparison of drought severity (Alley, 1984; Guttman et al., 1992), overly simplified soil model (Xu et al., 2012), overestimation of drought duration (Yu et al., 2019), and most importantly, overly large sensitivity to temperature, which tends to exaggerate drought severity in a warming climate (Sheffield et al., 2012). In this review, keeping in mind the above limitations, we nevertheless discuss drought studies which used the PDSI, simply because the limited number of existing reconstructions of paleoclimate drought atlas for different continents mostly used the PDSI rather than precipitation. Evaluations of the drought indices revealed that rainfall deciles and rainfall anomaly are among the best indices for meteorological drought (Keyantash & Dracup, 2002; Oladipo, 1985).

Since meteorological drought is shortage of precipitation, studies of the meteorological drought are intimately tied to precipitation variability over the continents. Precipitation on Earth can be generated by different physical mechanisms. Over the tropical continents and the warm-season extratropical continents, precipitation is often generated by deep convection. Over the cold-season extratropical continents, precipitation is often generated by fronts and orographic effects (see review by Houze, 2014; Lin et al., 2022). Figure 3b displays global map showing the season with the largest climatological precipitation at each grid point based on the GPCC global land precipitation dataset covering the period 1901 till present. Red colour (northern summer, June–July–August) and brown colour (northern autumn, September–October–November) are for the northern hemisphere warm season, while blue colour (northern winter, December–January–February) and green colour (northern spring, March–April–May) are for the northern hemisphere cold season. Magenta lines show the monsoon regions with summer-winter precipitation difference larger than 300 mm. Over most of the global continents, especially over the monsoon regions, the largest climatological precipitation falls during the local warm season. However, over the Middle East, Central Asia, the western United States, the southeastern United States, the south coast of Australia and Chile, the largest climatological precipitation falls during the local cold season. Shortages of precipitation during those seasons have been the priority of most drought studies and predictions, although deficits in precipitation in the other (transition) seasons would also contribute to the development of droughts.

Predicting droughts is one of the most difficult problems in climate science, because the occurrence of droughts involves various modes of climate variability, climate change, SST-driven teleconnections, atmospheric internal dynamics, land-atmosphere feedbacks and external forcings (Fig. 3a). State-of-the-art global climate models are used for drought forecasts. Currently, 3-month official forecasts and 12-month experimental outlooks are being issued by the leading global climate forecasting centres, such as the U. S. National Center for Environmental Prediction (NCEP; https://www.cpc.ncep.noaa.gov/products/predictions/90day/), the European Center for Medium-range Weather Forecast (ECMWF; https://apps.ecmwf.int/webapps/opencharts/products/seasonal_system5_standard_rain), the International Research Institute (IRI; https://iri.columbia.edu/our-expertise/climate/forecasts/seasonal-climate-forecasts/), and the Australian Bureau of Meteorology (BOM; http://www.bom.gov.au/climate/outlooks/#/rainfall/median/seasonal/0). For the 3-month forecasts, the skills are generally higher in the tropics, especially over Australia, Africa and South America in the local warm seasons, but lower in the extratropics over North America and Europe (Goddard et al., 2003, their Table 1). In fact, droughts are connected to climate variability at all time-scales. The seasonal monsoons produce most of the land-surface precipitation in the tropics and subtropics (Fig. 3b), and droughts in those regions are often associated with the failure of monsoon rains. The global monsoon constitutes the rising branch of the Hadley Circulation, and thus some droughts are associated with variations in the seasonal zonal-mean atmospheric circulation. Droughts are modulated by the El Nino-Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), the Atlantic Multi-decadal Oscillation (AMO), and the Global Inter-Centennial Oscillation (GICO), among other modes of climate variability. Droughts are also affected by local feedback processes such as the circulation-convection feedback and land-atmosphere feedback. Moreover, external forcings such as human-released greenhouse gases, solar forcing, volcano forcing and tidal gravitational force can impact droughts.

This paper reviews the observations and theories of meteorological droughts. Based on the time-scales, we divide droughts into two groups: (1) the interannual to decadal droughts which last for less than 20 years (hereafter called droughts), and (2) the mega-droughts which last longer than 20 years. The latter are further divided into multi-decadal and centennial mega-droughts. We use 20-years as a somewhat arbitrary separation between droughts and mega-droughts, though this is supported somewhat by observations/reanalyses showing that the global climate system has a natural separation of variance between timescales longer than 40 years and shorter than 40 years (see Fig. 15). With the oscillations of precipitation that are longer than 40 years, the wet/dry periods last for more than 20 years. We do not discuss the short-term “flash droughts” at subseasonal/seasonal time scales, such as the droughts/heat waves in Northern Eurasia (e.g. Schubert et al., 2014), which generally do not cause as much damage as the long-term droughts. In the following, we review separately the history of research on droughts and mega-droughts.

2 Droughts

a Overview

Rainfall began to be measured in the ancient times by different civilizations, with records of the simplest rain gauges going back as far as the sixth century BC in ancient Greece (Gallicchio, 2017). A bowl with its mouth as wide as 457 mm was used as a rain gauge in India in the fourth century BC (Frisinger, 1977; Shamasastry, 1915). In China, rain gauges were installed across the country around 1247 and methods were developed to convert the point measurements to area averages (Needham, 1959). In Korea, rain gauges were distributed to the country around 1441 and were used continuously until 1907 (Wada, 1911). In Italy, Benedetto Castelli, a long-term friend of Galileo, used a glass to measure rainfall around 1639 (Strangeways, 2010).

The modern instruments for measuring rainfall were invented in the seventeenth century (Strangeways, 2010). Wren & Hooke developed a tipping bucket (Grew, 1681), and Towneley (1694) published the first systematic measurement of rainfall. Starting from 1861, the British Rainfall Organization and UK Met Office published 126 volumes of “British Rainfall”, with about 500 stations in 1861 (Burt, 2010). Stations were also set up in the late nineteenth century in the United States and the British colonies, such as India, South Africa, and Australia, and then spread to other countries in the early twentieth century. Today, numerous surface stations exist around the world (Kidd et al., 2017), while satellites with different sensors provide global coverage especially over the ocean (Kidd & Levizzani, 2011). The longest instrument global precipitation dataset goes back to 1900 (Sun et al., 2018).

Research on droughts began in the late nineteenth century for various regions around the world, such as India (Blanford, 1884), South Africa (Hutchins, 1889), North America (Guinn, 1890; Henry, 1930a, 1930b), Australia (Berney, 1906), and Brazil (Freise, 1938; Mossman, 1919). Brooks and Glasspoole (1922) presented an excellent pioneering study on the relationships between severe droughts in the British Isles and large-scale atmospheric circulation. They recognized that the establishment and maintenance of local anticyclonic conditions is closely associated with dry weather, and aimed to discover how the local anomalies of pressure distribution are related to others in the general circulation in different parts of the world. Using 72 years (1850–1921) of data, they identified nine severe drought periods, which can be categorized into two types. Type A is characterized by high pressure anomalies to the north of the British Isles. The Icelandic low is shifted toward west Greenland and the Azores high is weakly developed. In this case, the drought is said to be due to dry winds rather than to anticyclonic conditions. Type B is characterized by high pressure anomalies right over the British Isles, and low pressure anomalies over the Arctic Ocean and at the Azores (Fig. 4a), which is similar to the North Atlantic Oscillation (NAO) pattern discovered by Exner (1913).

Fig. 4 (a) Mean pressure anomaly of seven droughts in British Isles (from Brooks & Glasspoole, 1922). (b) Change in number of rain days per year associated with decrease of 10% in forest area (from Brooks, 1928).

Brooks (1928) published a pioneering study on land-atmosphere interaction and effects of deforestation and afforestation. There was a hot debate at that time on the influence of forests on rainfall. Brooks analysed the vegetation data and rainfall data from a deforestation/afforestation region (Fig. 4), and found that the replacement of forests by field crops or grass land probably tends to increase the general rainfall of the district slightly, while their replacement by bare ground probably tends to decrease the general rainfall slightly. Forests increase the local rainfall by one to two per cent, which is due to the increase of surface roughness. The most noticeable effect of afforestation of grass or crop land is increase of the local run-off by ten to twenty per cent, which is partly due to the decreased loss of water in evaporation and transpiration, and partly due to the slight increase of actual rainfall. Afforestation of bare ground would decrease the run-off, but this effect would be masked by the much greater regularity of the flow. Brooks suggested that forests, drawing their water supply from deep sources, possibly exercise a slight mitigating effect on droughts.

Since then, because of the strong hydrological, agricultural and socio-economic consequences of droughts, numerous case studies have been conducted on major drought events in different regions of the globe, such as North America (Chang & Wallace, 1987; Chen & Newman, 1998; Hoerling et al., 2014; Hoerling & Kumar, 2003; Klein, 1952; Liu et al., 1998; Namias, 1955, 1960, 1966, 1978, 1982, 1991; Palmer & Brankovic, 1989; Schubert et al., 2004b; Seager, 2007; Seager, , Kushnir, et al., 2008; Seager, Burgman, et al., 2008; Trenberth et al., 1988; Trenberth & Guillemot, 1996; Wang et al., 2014), Europe (Green, 1977; Hurrell, 1995; Kingston et al., 2013, 2015; Obukhov et al., 1984; Polozova & Grigoryeva, 1984; Ratcliffe, 1978; Schubert et al., 2011; Trenberth & Fasullo, 2012), Africa and Middle East (Barlow et al., 2016; Bryson, 1973; Folland et al., 1986; Fontaine et al., 1995; Giannini et al., 2003; Grist & Nicholson, 2001; Held et al., 2005; Lu & Delworth, 2005; Lyon, 2014; Nicholson, 1980, 1981; Rodriguez-Fonseca et al., 2011, 2015; Rowell, 2001; Shinoda, 1990; Yang et al., 2014; Winstanley, 1973a, 1973b), Latin America (Chu, 1983; Chung, 1982; Coelho et al., 2016; Enfield, 1996; Enfield & Alfaro, 1999; Grimm, 2011; Grimm & Zilli, 2009; Magana et al., 1999; Markham & McLain, 1977; Moron et al., 1995; Moura & Shukla, 1981; Namias, 1972; Yoon & Zeng, 2010), East Asia (Guo, 1992; Park & Schubert, 1997; Simmonds et al., 1999; Xu & Zhang, 1983; Ye, 1981; Zhang & Zhou, 2015; Zhou & Yu, 2005), and India (Berlage, 1957; Bhalme et al., 1983; Bhalme & Mooley, 1980; Goswami, 1998; Kanamitsu & Krishnamurti, 1978; Keshavamurty & Awade, 1974; Krishnan et al., 2006, 2009; Mishra et al., 2019; Raman & Rao, 1981; Ramaswamy, 1962, 1976; Sikka, 1980; Walker, 1923, 1924, 1928, 1937; Walker & Bliss, 1930, 1932). The above case studies provided detailed analysis of the physical processes associated with the development of individual drought events, which has led to interesting new ideas and new hypothesis on the underlying physical mechanisms. As just one example, the 1988 North American drought, which caused $60 billion damage, has been examined extensively. In one of the earliest studies of that drought, Trenberth et al. (1988) found strong anticyclonic conditions and a northward displaced jet stream in the upper atmosphere over North America throughout this period were part of a distinctive wavetrain that emanated from the tropical Pacific. They further found that anomalies associated with the displaced ITCZ can force an anomalous wavetrain across North America similar to that observed. They proposed that the tropical Pacific SST anomalies are a notable (but not the sole) external forcing of the patterns. However, the generality of such mechanisms needs to be confirmed by looking more broadly across an ensemble of drought events, and carrying out additional modelling studies. In fact, Schubert et al. (2009) provide evidence from simulations with an ensemble of atmospheric general circulation models (AGCMs) forced by idealized SST, that the precipitation deficits associated with droughts such as the one that occurred over the U.S. Great Plains during 1988, tend to be associated with cold tropical Pacific SSTs and warm tropical Atlantic SSTs, with the opposite (warm tropical Pacific SSTs and cold tropical Atlantic SSTs) true for wet conditions (see also next section).

With the advances in global precipitation datasets, SST analysis and atmospheric reanalysis, climatology and composites of droughts have been presented for different continents, such as North America (Andreadis et al., 2005; Karl & Koscielny, 1982; Seager & Hoerling, 2014), Europe (Briffa et al., 1994; Schubert et al., 2014; Spinoni et al., 2015), Africa (Nicholson et al., 2012; Rodriguez-Fonseca et al., 2011), South America (Grimm, 2011; Grimm & Zilli, 2009; Marengo & Espinoza, 2016), Asia (Mishra et al., 2019; Schubert et al., 2014; Wang et al., 2011), and the whole globe (Dai, 2011; Dai et al., 1998; Sheffield et al., 2009). These composite studies provided statistically robust results on the physical processes and physical mechanisms, which are reviewed below. Long-term ensemble climate model experiments further help to decompose the physical processes and feedbacks, which cannot be easily isolated in observational analysis.

The proposed physical mechanisms for droughts include (1) SST-driven teleconnections, (2) atmospheric internal dynamics, and (3) land-atmosphere feedback. Generally, droughts in the tropics and subtropics are better correlated with SST anomalies, while droughts in the mid-latitudes and high latitudes are more related to atmospheric internal dynamics with weaker but nevertheless important ties to remote SST anomalies. Land-atmosphere feedbacks often act in tandem with the other mechanisms leading to longer and/or more extreme events. Each of these mechanisms will be discussed in more detail below.

b SST-Driven Teleconnections

The first mechanism is SST-driven teleconnections. Following the methods of Schubert et al. (2016), Fig. 5 shows the correlations of the 2–10-year filtered monthly global SST anomalies in the past 116 years (1901–2016) from the ERSST dataset (Huang et al., 2017) with the 2–10-year filtered area-averaged precipitation anomalies from the GPCC dataset (Schneider et al., 2016) for (a) Asia, (b) Oceania, (c) northern South America, (d) Africa, (e) North America, (f) Middle East/southern Europe, (g) southern South America, and (h) northern Europe. Black stars denote the grids with correlation coefficient above the 95% confidence level. Magenta lines denote the region for precipitation averaging. The interannual precipitation anomalies over Asia are significantly correlated with a negative ENSO SST pattern in the Pacific Ocean, and a negative Indian Ocean Dipole (IOD; Saji et al., 1999; Webster et al., 1999) pattern in the tropical Indian Ocean (Fig. 5a). The precipitation anomalies over Oceania are associated with a negative ENSO SST pattern in the Pacific Ocean, basin-wide cooling in the tropical Indian Ocean, cooling symmetric to the equator in the tropical Atlantic Ocean, and warming in the South Atlantic Ocean (Fig. 5b). The precipitation anomalies over northern South America are associated with very similar global SST patterns as Oceania (Fig. 5c). The precipitation anomalies over Africa are associated with negative SST anomalies in the tropical central/eastern Pacific and the tropical Indian Ocean (Fig. 5d). The precipitation anomalies over North America are associated with a positive ENSO SST pattern in the Pacific Ocean, warming in the Atlantic warm pool (Wang & Enfield, 2001; C. Wang, Enfield, et al., 2006), and basin-wide warming in the tropical Indian Ocean (Fig. 5e). The precipitation anomalies over Middle East and southern Europe are associated with very similar global SST patterns as North America but with stronger correlations (Fig. 5f). The precipitation anomalies over southern South America are also associated with very similar SST patterns as North America in the Pacific and Indian Oceans with stronger correlations, but different SST patterns in the Atlantic Ocean with warming in the South Atlantic Ocean (Fig. 5g). The precipitation anomalies over northern Europe are associated with a dipole SST pattern in the North Atlantic Ocean and warming in the northwest Pacific Ocean (Fig. 5h). Overall, the wet years over Asia, Oceania, northern South America and Africa are associated with La Nina in the Pacific Ocean and general cooling in the tropical Indian Ocean, while the wet years over North America, Middle East, southern Europe and southern South America are associated with El Nino in the Pacific Ocean and general warming of the Indian Ocean. Note that basin-wide warming (cooling) of the tropical Indian Ocean occurs in the summer following the peak of El Nino (La Nina).

Fig. 5 Correlation of the 2–10-year filtered monthly global SST anomalies in the past 116 years (1901–2016) from the ERSST dataset with the 2–10-year filtered monthly area-averaged precipitation anomalies from the GPCC dataset for (a) Asia, (b) Oceania, (c) northern South America, (d) Africa, (e) North America, (f) Middle East/southern Europe, (g) southern South America, and (h) northern Europe. Black stars denote the grids with correlation coefficient above the 95% confidence level. Magenta lines denote the region for precipitation averaging.

An interesting question is what the typical spatial scale of droughts is. When presenting the precipitation or drought index time series, the choice of the region to average is often arbitrary or based on continents or administrative regions. An objective method for identifying the spatial scale of droughts is the EOF analysis, which can break the continent-averaged precipitation into the leading modes of droughts. Rotation of EOF analysis could further alleviate the artificial effect of mathematical orthogonality. EOF analysis and rotated EOF analysis have been applied to Africa (Rodriguez-Fonseca et al., 2011), South America (Grimm, 2011; Grimm & Zilli, 2009), North America (Karl & Koscielny, 1982; Seager & Hoerling, 2014), and Europe (Briffa et al., 1994; Schubert et al., 2014). The spatial scales of the interannual droughts are generally larger than 2000 km (see Figs 6 and 7).

Fig. 6 (a) The leading modes of interannual precipitation anomalies in Western African monsoon and their correlations with SST anomalies (from Rodriguez-Fonseca et al., 2011). (b) The leading modes of interannual precipitation anomalies in South America and their correlations with SST anomalies (from Grimm & Zilli, 2009).

Fig. 7 (a) The leading EOF patterns of precipitation in North America (upper) and the correlation with SST anomalies (lower) (from Seager & Hoerling, 2014). (b) The temporal correlation between the monthly JJA 250-mb meridional wind and leading RPCs of combined 2-m temperature and precipitation in Europe (left) and the corresponding precipitation patterns (right) (from Schubert et al., 2014).

The SST anomalies play a key role in driving the regional droughts over the tropical and subtropical continents, such as Africa, South America, Oceania and Asia, but the SST-driven teleconnections become weaker over the mid-latitude and high-latitude continents, such as North America and Europe. When calculating the correlations of the leading modes with global SST anomalies, mode 1 of the interannual precipitation anomalies in west Africa is highly correlated with the ENSO, while mode 2 is correlated with the tropical Atlantic SST anomalies (Fig. 6a; Rodriguez-Fonseca et al., 2011). For South America, mode 1 and mode 3 are highly correlated with the ENSO, while mode 2 is related to the South Atlantic SST anomalies (Fig. 6b; Grimm & Zilli, 2009). For North America, an unrotated EOF analysis shows that mode 1 is strongly correlated with the ENSO, but mode 2 and mode 3 are not forced by SST anomalies (Fig. 7a; Seager & Hoerling, 2014). For northern Eurasia, none of the leading summer rotated EOF modes is significantly correlated with SST anomalies (Fig. 7b; Schubert et al., 2014).

Theoretical studies showed that teleconnections can be forced by SST anomalies in both the tropics and the extratropics (though less so). Land-sea temperature contrast was first proposed to drive the monsoons in the tropics (Halley, 1686; Holton, 1992; Hoskins & Rodwell, 1995; Rodwell & Hoskins, 1995; Wallace & Hobbs, 1977; Webster, 1987; Webster et al., 1998). The tropical atmospheric responses to tropical SST anomalies can be understood using the Gill model (Gill, 1980; Zhang & Krishnamurti, 1996). Diabatic heating symmetric to the equator generates a Kelvin wave response to the east and a Rossby wave response to the west, while heating anti-symmetric about the equator produces a mixed Rossby-gravity wave and a Rossby wave close to the heating. Zhang and Krishnamurti (1996) showed that a linear sum of the symmetric and anti-symmetric heating responses can explain quite well the observed monthly tropical circulation. The extratropical atmospheric response to SST forcing is characterized by propagation of Rossby waves (Bjerknes, 1937; Bolin, 1950; Charney & Eliassen, 1949; Held et al., 1985; Hoskins & Karoly, 1981; Rossby, 1939; Smagorinsky, 1953; see reviews by Held, 2019; Trenberth et al., 1998). The Rossby waves tend to propagate through the westerly waveguides and thus are very sensitive to seasonal mean state of westerly flow (Ambrizzi et al., 1995; Ambrizzi & Hoskins, 1997; Hoskins & Ambrizzi, 1993). The global distribution of Rossby waveguides and preferred propagation paths have been shown by Ambrizzi et al. (1995) for JJA, and by Hoskins and Ambrizzi (1993) for DJF. The Rossby waveguides are basically the jet streams, with the warm season hemisphere having one subtropical jet while the cold season hemisphere having a subtropical jet and a subpolar jet. The preferred propagation paths of Rossby waves connect the tropics with the extratropics. For example, in DJF, SST anomalies in the tropical Pacific force a wave train similar to the Pacific-North American Pattern (Hoskins & Karoly, 1981), while SST anomalies in Indian Ocean drive a pattern heavily projected onto the North Atlantic Oscillation pattern (Hoerling et al., 2004). SST anomalies in the North Atlantic force a dipole between Greenland and Europe (Kushnir & Held, 1996), whereas SST anomalies in the North Pacific drive a dipole between Aleutian low and Hawaiian high (Kushnir et al., 2002), which is quite similar to the North Atlantic case. As shown in Fig. 3b, over most of the mid-latitude and high-latitude continents, the largest climatological precipitation occurs in the warm season when background westerly flow is weaker, leading to a weaker teleconnection driven by SST forcing, which may explain the weaker SST correlations for droughts over northern Europe (Figs 5h and 7b). However, over some of the subtropical continents such as the western United States, southeastern United States, Middle East and southern Europe, the largest climatological precipitation occurs in the cold season when background westerly flow (the subtropical jet) is stronger, leading to strong teleconnections driven by SST forcing (Fig. 5e, f).

The dominant interannual mode of the global climate system is the ENSO, whose discovery was stimulated by the 1875–1878 Great Drought and resultant Global Famine, which killed more than 50 million people around the world (Blanford, 1884; Hildebrandsson, 1897; Lockyer & Lockyer, 1902; Walker, 1923). Figure 8a summarizes the observed teleconnections of ENSO to the extratropics, tropical Indian and Atlantic Oceans, tropospheric zonal mean circulation and stratospheric circulation. El Niño drives the positive barotropic PNA pattern to northern hemisphere extratropics leading to a negative NAO (Horel & Wallace, 1981; Li & Lau, 2012), the positive barotropic PSA pattern to southern hemisphere extratropics affecting the Antarctic circumpolar wave (Karoly, 1989; Peterson & White, 1998), and the positive PJ/AEP pattern to east Asia which is a highly populated region (Nitta, 1986). Within the tropics, El Niño suppressed the Walker circulation and drives teleconnections to tropical Indian and Atlantic Oceans through atmospheric bridge, ocean dynamics and local air–sea feedbacks, leading to SST warming 3–6 months after the eastern Pacific warming (Klein et al., 1999). The Indian Ocean warming in turn serves as a heat capacitor to enhance the PJ/EAP pattern (Yang et al., 2007). In addition to the basin-wide warming, ENSO is also linked to the dipole modes such as IOD (Ashok et al., 2003) and AMM (Enfield & Mayer, 1997). El Niño enhances the tropospheric Hadley circulation and Ferrel circulation in certain seasons (Feng & Li, 2013; Oort & Yienger, 1996), and the stratospheric Brewer-Dobson circulation which leads to colder temperature and reduced ozone in tropical lower stratosphere, but warmer temperature and enhanced ozone in polar vortex region (Randel et al., 2009; Randel & Cobb, 1994). El Niño weakens the polar vortices in both hemispheres, while the northern polar vortex further affects climate in Europe by downward control (Domeisen et al., 2019; van Loon & Labitzke, 1987).

Fig. 8 (A) Schematic depiction of global teleconnections of ENSO. Walker circulation is represented by dashed cells. The PNA and PSA patterns are illustrated by centres of geopotential height anomalies and grey arching arrow towards North and South America, respectively (red – positive geopotential anomaly; blue – negative geopotential anomaly). The East Asia–Pacific (EAP) pattern or Japan – Pacific (JP) pattern is indicated by orange (positive) and light blue (negative) centres of sea level pressure. The Northwest Pacific Anticyclone (NWP AC) is represented by the orange circle. The ‘capacitor effect’ is represented by purple dashed arrow. Shading in the oceans: Sea surface temperature anomalies (red – warm, blue – cool). Shading over land: precipitation anomalies (brown – dry, green – wet). A positive Indian Ocean Dipole (IOD+) event can co-occur with El Niño in austral spring. An Indian Ocean Basin-wide warming (IOBW) and a tropical North Atlantic (tNA) warming can develop during and after El Niño’s peak. The tropical Atlantic can exhibit a positive Atlantic Meridional Mode pattern (AMM+). The zonal mean circulation and El Niño teleconnections via stratospheric pathway are illustrated by arrows on the left side of the schematic. (B) Schematic depiction of the observed physical mechanisms leading to the switch between El Niño and La Niña (from Lin & Qian, 2019b).

What drives the ENSO? Many ENSO theories have been developed, all involving ocean-atmosphere feedbacks. The most important feedback is the SST gradient-trade wind feedback discovered by Bjerknes (1969), with the SST gradient along the equator driving the Walker circulation, and the Walker circulation, in turn, enhancing the SST gradient by forcing equatorial ocean upwelling in the eastern Pacific cold tongue. In addition to the Bjerknes feedback, other feedback mechanisms also play important roles in the Walker Circulation and the ENSO (Lin, 2007), such as the SST-latent heat flux feedback (Wallace, 1992) and the SST-shortwave flux feedback (Ramanathan & Collins, 1991). The switch between El Nino and La Nina has been linked to the ocean equatorial waves (e.g. Picaut et al., 1997; Suarez & Schopf, 1988; Weisberg & Wang, 1997), or recharge of ocean heat content (Jin, 1997), or stochastic forcings (e.g. Lau, 1985; Penland & Sardeshmuhk, 1995). Recently, Lin and Qian (2019b) found that the switch between El Niño and La Niña is caused by a subsurface ocean wave propagating from the western Pacific to the central and eastern Pacific and then triggering development of SST anomaly. Further analysis of Earth's angular momentum budget and NASA's Apollo Lunar Laser Ranging Experiment suggests that the subsurface wave is likely driven by lunar tidal gravitational force (Fig. 8b).

ENSO SST anomalies often force strong global impacts in precipitation. Global climate models forced by the observed SST anomalies generally could reproduce the observed precipitation anomalies over the tropical and subtropical continents, but tend to underestimate the precipitation anomalies over the mid-latitude and high-latitude continents, although inter-model differences exist. Schubert et al. (2016) reviewed the SST-forced teleconnections associated with droughts in eight regions around the world (Fig. 9a). The synthesis is primarily based on regionally focused articles submitted to the Global Drought Information System (GDIS) collection together with new results from a suite of atmospheric general circulation model experiments intended to integrate those studies into a coherent view of drought worldwide. On interannual time scales, the preeminence of ENSO as a driver of meteorological drought throughout much of the Americas, eastern Asia, Australia, and the Maritime Continent is well established, whereas in other regions (e.g. Europe, Africa, and India), the response to ENSO is more ephemeral or nonexistent. Northern Eurasia, central Europe, and central and eastern Canada stand out as regions with few SST forced impacts on precipitation on interannual time scales. Many AGCMs could simulate well the precipitation anomalies over different continents. However, some big drought events remain difficult to simulate, such as the 1930s Dust Bowl in North America (Seager, Kushnir, et al., 2008). Moreover, intercomparisons of AGCMs’ AMIP experiments showed that inter-model differences exist. Tippett (2006) compared seven AGCM simulations of Sahel precipitation, and found that the models tend to simulate better decadal variability than interannual variability of Sahel precipitation. Mohino et al. (2011) evaluated the AMIP simulations of five AGCMs and found that not all models were able to reproduce the observed negative link between rainfall over the Sahel and El Nino-like anomalies in the Pacific. When the AGCMs are coupled to the OGCMs, the model biases become larger. Ukkola et al. (2018) evaluated the historical simulations (1950–2014) for 20 coupled models from the Coupled Model Intercomparison Project (CMIP5), and found that at the global scale, the coupled models agree well on most precipitation drought metrics, but 1980systematically underestimated precipitation drought intensity compared to observations.

Fig. 9 (a) Background map shows the measure of impact of the prescribed SSTs on the precipitation. The inset maps show the correlations between the area-averaged precipitation with SST. All results are for 1979–2011 based on 60 ensemble members: 12 AMIP simulations for five models (from Schubert et al. 2016). (b) Global impacts of the ENSO lifecycle (from Lin & Qian, 2019a).

The observed impacts of ENSO on global continental precipitation were first studied in the s by Ropelewski and Halpert (1986, 1987, 1989) using the World Monthly Surface Station Climatology (WMSSC) dataset compiled by the National Center for Atmospheric Research (NCAR), which were summarized later by the review of Trenberth et al. (1998). Recently, Lin and Qian (2019a) demonstrate a new picture of the global impacts of ENSO throughout its whole lifecycle based on the latest satellite, in situ and reanalysis datasets (Fig. 9b). The ENSO lifecycle significantly affects precipitation over the tropical and subtropical continents, and also affects precipitation over north Russia and north Canada. In general, the warm phase of ENSO is associated with dry anomalies over the tropical continents such as the west Africa, southeast Asia, Australia and tropical South America, wet anomalies over the subtropical continents such as the Middle East, central Asia, subtropical North America and subtropical South America, and dry anomalies over high-latitude continents such as north Russia and north Canada. The global impacts of the cold phase generally mirror those of the warm phase. This tri-belt seasaw among tropical continents, subtropical continents and high-latitude continents is likely linked to the ENSO-driven response of the global zonal mean atmospheric circulation, since the tropical belt is the ascending branch of the Hadley Cell while the subtropical belt is the descending branch of the Hadley cell and Ferrel cell, and the high-latitude belt is the ascending branch of the Ferrel cell and polar cell.

c Atmospheric Internal Dynamics

The second mechanism is atmospheric internal dynamics including monsoon dynamics (Chao & Chen, 2001a; Ramage, 1971; Riehl, 1954, 1979; Schneider & Bordoni, 2008; Sikka & Gadgil, 1980; Simpson, 1921; Trenberth et al., 2000), ITCZ dynamics (Bacmeister et al., 2006; Chao, 2000; Chao & Chen, 2001b, 2004; Charney, 1969, 1971; Holton et al., 1971; Lindzen, 1974; Schneider et al., 2014), Hadley circulation dynamics (Held & Hou, 1980; Hess et al., 1993; Numaguti, 1993; Plumb & Hou, 1992; Raymond, 2000; Schneider, 1977; Schneider & Lindzen, 1977), and extratropical Rossby wave dynamics (Schubert et al., 2011, 2014). Simpson (1921) first pointed out that land-sea temperature contrast cannot explain some observed key features of the monsoon: India has a much higher surface air temperature in May before the monsoon starts than in July when monsoon is at its peak, and the average temperature is much greater in dry years than in wet years. The northwest India, which is the hottest part of the country, does not get any precipitation during monsoon season. Trenberth et al. (2000) conducted a complex EOF analysis of the global divergent circulation and found a global monsoon mode which accounts for 60% of the variance. The global monsoon mode dominates the Hadley circulation, the Pacific Walker circulation, the Americas-Atlantic Walker circulation, and the Walker-type circulation over Asia, Middle East, North Africa, and the Indian Ocean. Now it is well recognized that monsoons are related to shift of global ITCZ and Hadley circulation.

For the extratropics, Schubert et al. (2011, 2014) for example found that summertime stationary Rossby waves play a key role in the development of droughts, flooding and heat waves, such as the 1988 US drought, the 1993 US flooding, the 2003 European heat waves, and the 2010 Russian heat waves. They show up as the leading modes of rotated EOF of monthly North Hemisphere surface temperature and precipitation (Fig. 7b). At intraseasonal time scales, the leading modes account for more than 30% of the precipitation variability and more than 60% of the surface temperature variability over the U.S. northern Great Plains, parts of Canada, Europe, and Russia (Schubert et al., 2011). The Rossby waves are guided by the mean-state jet streams. Model experiments revealed that the Rossby wave forcing is dominated by submonthly vorticity transients.

d Land-Atmosphere Feedback

The third mechanism is land-atmosphere feedback. Land-atmosphere feedback studies were pioneered by Brooks (1928), Newell (1971) and Charney (1975), Charney et al. (1975). The proposed physical processes include feedbacks through albedo (Charney, 1975; Charney et al., 1975, 1977; Sud & Fennessy, 1982), surface roughness (Sud & Smith, 1985), soil moisture (Guo et al., 2006; Koster et al., 2003, 2004, 2006; Oglesby & Erickson III, 1989; Shukla & Mintz, 1982; Sud & Fennessy, 1984; Walker & Rowntree, 1977), and vegetation (B. I. Cook et al., 2009; Dirmeyer, 1994; Meng et al., 2014; Pielke et al., 1998; W. Pitman et al., 2009, 2012; Schubert et al., 2004a; Zeng et al., 1999; Zeng & Neelin, 2000; Wang, Anderson, et al., 2006). Charney (1975) conducted GCM experiments with different prescribed surface albedo in the Sahara region, and found that an increased albedo causes decreased precipitation. He proposed that decreased plant cover during drought will lead to increased albedo and then be reinforced by decreased precipitation. This positive land-atmosphere feedback could initiate or perpetuate a drought. Koster et al. (2004) used a dozen climate models to identify the key regions on Earth where precipitation is affected by soil moisture anomalies during Northern Hemisphere summer (Fig. 10a). They found that the major sensitive regions are the central Great Plains of North America, the Sahel, equatorial Africa, and India. There are secondary regions over South America, central Asia, and China. The major sensitive regions are in the transition zones between wet and dry climates. They proposed that precipitation is not sensitive to soil moisture either in wet regions where evaporation is controlled by radiation, or in dry regions where evaporation is too small to affect precipitation. In the transition zones, however, evaporation is still sensitive to soil moisture but high enough to affect boundary-layer moisture, and boundary-layer moisture can trigger moist convection, which connect soil moisture to precipitation. Zeng et al. (1999) developed a theoretical model of Sahel to show that interactive vegetation and associated land-surface feedback enhance precipitation variability at interannual and interdecadal time scales. They found that vegetation changes associated with changes in rainfall during the multidecadal drying period from the 1950s to the 1980s.

Fig. 10 (a) The land-atmosphere coupling strength diagnostic for boreal summer averaged across the 12 models participating in GLACE. Insets are areally averaged coupling strengths for the 12 individual models over the outlined, representative hotspot regions (from Koster et al., 2004). (b) Relationship of NLH to SWet in the 16 ensemble members of nine GCMs at the grid box encompassing the ARM Central Facility. Solid blue line is fit through the means of 20 bins of equal number of points. Red points show the ensemble member used as basis for fixed SWet integrations (from Dirmeyer et al., 2006).

Simulations of land-atmosphere feedbacks are strongly affected by the land surface schemes in climate models (Deardorff, 1978; Dickinson et al., 1986; Manabe, 1969; Sellers et al., 1986; see reviews by Avissar & Verstraete, 1990; Betts et al., 1996; Peters-Lidard et al., 2018; Pitman, 2003; Randall et al., 2019). It is important to evaluate the model simulations using observations (Dirmeyer, 2011; Dirmeyer et al., 2006, 2016, 2018). Dirmeyer et al. (2006) examined the local covariability of key atmospheric and land surface variables both in models and in observations from the Atmospheric Radiation Measurement (ARM) stations (Fig. 10b). Most models cannot reproduce well the observed relationships between surface and atmospheric state variables and fluxes, suggesting that these models do not capture correctly the land-atmosphere feedbacks. They found systematic biases in surface air temperature and humidity which may contribute to the unrealistic surface flux sensitivities. The multimodel ensemble mean generally performs better than most or all of the individual models. The models also show a large scatter in regional precipitation characteristics, but the multimodel ensemble mean is remarkably accurate.

Physics-based evaluations of climate models are important for understanding the contributions of the three physical mechanisms to droughts. Modelling studies accessing the relative contributions of the different mechanisms often compare only with observed precipitation, which is the end product of all the physical processes. When model precipitation matches well with the observed precipitation, it does not guarantee that the physical processes are correct. For example, the correct model precipitation may come from an overestimation of SST response offset by an underestimation of land-atmosphere feedbacks. Therefore, it is important to evaluate step-by-step the whole chain of physical processes, as well as the strength of the feedbacks.

3 Mega-droughts

a Overview

Due to their longer durations, the multi-decadal- to centennial-scale mega-droughts have more profound impacts on human society and human history. Multi-decadal- to centennial-scale mega-droughts were linked to the collapse of the Old Kingdom in Egypt around 4200 cal BP (Fig. 1b, Barta & Bezdek, 2008; Marriner et al., 2012), the collapse of Western Roman Empire in Europe during the Migration Period in 300–400 AD (Fig. 1a; Buntgen et al., 2011), the collapse of the Tang, Yuan, and Ming Dynasties in China (Fig. 1b; Zhang et al., 2008), and the development and disintegration of Maya political systems in North America (Evans et al., 2018; Kennett et al., 2012). Rapid progress is being made in recent years in research on mega-droughts thanks to the advances in high-resolution paleoclimate records, paleoclimate field reconstructions, and paleoclimate reanalysis.

High-resolution paleoclimate data with annual to decadal resolutions are needed for studying the mega-droughts. Tree rings provide nice records of precipitation for the common era, while speleothem and ice core provide good histories for Holocene and late Quaternary. Quantitative tree ring studies began in the late nineteenth century and early twentieth century (Abbot, 1936; Douglass, 1909, 1914; Huber, 1935; Kapteyn, 1914; Keen, 1937; Schulman, 1947; Shvedov, 1892). Shvedov (1892) presented a pioneering precipitation analysis using two trees in Russia, and found a 3–9 year cycle in precipitation, making it possible to correctly predict the droughts in 1882 and 1891. Douglass (1914) also developed a method of estimating rainfall from tree ring data. Starting from the 1950s, radioisotopic dating and stable isotope analysis were applied to speleothems (Broecker et al., 1960; Orr, 1952, 1953) and ice cores (Coachman et al., 1958; Dansgaard et al., 1969; Langway Jr. et al., 1965; Scholander et al., 1962).

The past three decades (1990–2020) witnessed a major increase in high-resolution paleoclimate datasets. Numerous high-quality high-resolution datasets have been developed for Greenland (Alley, 2000; NGRIP members, 2004; Vinther et al., 2008, 2010), Antarctica (Abram et al., 2013; Steig et al., 2013), Africa (Scroxton et al., 2017; Shanahan et al., 2009; Thompson et al., 2002; Voarintsoa et al., 2017), Asia (Eichler et al., 2009; Griffiths et al., 2016; Kaspari et al., 2007; Sinha, Berkelhammer, et al., 2011; Sinha, Stott, et al., 2011; Tan et al., 2003, 2009; Thompson et al., 1997, 2000; Treydte et al., 2006; Wang et al., 2005), Australia (Denniston et al., 2017; Saunders et al., 2013), Europe (Baker et al., 2011; Edwards et al., 2017; Fohlmeister et al., 2013; Isaksson et al., 2005; Kinnard et al., 2011; Voronin et al., 2008), North America (Anderson, 2011; Asmerom et al., 2013; Bird et al., 2017; Buhay et al., 2009; Edwards et al., 2008; Kennett et al., 2012; Kinnard et al., 2011; Lachniet et al., 2017; McCabe-Glynn et al., 2013; Medina-Elizalde et al., 2010; Steinman et al., 2012), South America (Bird et al., 2011; Elbert et al., 2013; Kanner et al., 2013; Kellerhals et al., 2010; Novello et al., 2016; Reuter et al., 2009; Thompson et al., 2013; Wortham et al., 2017), Indian Ocean (Doose-Rolinski et al., 2001; Munz et al., 2015), tropical western Pacific (Maupin et al., 2014; Partin et al., 2013), Northwestern Pacific (Isono et al., 2009), North Atlantic (Cronin et al., 2010; Cunningham et al., 2013; Eirìksson et al., 2006; Keigwin, 1996; Reynolds et al., 2016; Sicre et al., 2008; van Hengstum et al., 2015), and tropical Atlantic (Black et al., 2007; Kilbourne et al., 2008; McGregor et al., 2015). The datasets have been summarized by a series of reviews (Jones et al., 2009; Ljungqvist, 2009; Ljungqvist et al., 2012, 2016; PAGES 2k Consortium, 2017).

Using the high-resolution datasets, paleoclimate field reconstructions and paleoclimate reanalysis have been developed to produce spatial distribution of temperature and hydroclimate (Cook & Jacoby, 1977; Fritts, 1965; Stockton & Meko, 1975). Temperature reconstructions have been generated for Europe (Briffa et al., 1990; Ljungqvist, Seim, et al., 2019; Luterbacher et al., 2016), Asia (Cook, Krusic et al., 2013; Zhang et al., 2018), North America (Anchukaitis et al., 2013; Trouet et al., 2013), Australia (Gergis et al., 2016), Arctic (Helama et al., 2012; Werner et al., 2018), and northern hemisphere (Briffa, 2000; Büntgen et al., 2021; Cook, Esper et al., 2004; D'Arrigo et al., 2006; Esper et al., 2002, 2012; Guillet et al., 2017; Jones et al., 1998, 2009; Ljungqvist, 2010; Mann et al., 1999; Moberg et al., 2005; Schneider et al., 2015; Stoffel et al., 2015; Wilson et al., 2016). Drought and precipitation have been reconstructed for North America (B. I. Cook et al., 2018; E. R. Cook et al., 1999, 2007, 2020; Cook, Woodhouse et al. 2004; E. R. Cook, Seager, et al., 2010; D'Arrigo & Jacoby, 1991; Fritts, 1965; Routson et al., 2011; Schulman, 1947; Stahle et al., 1985, 2000; Stockton & Meko, 1975; Weakly, 1965), Asia (Borgaonkar et al., 2010; Buckley et al., 2010; E. R. Cook, Anchukaitis, et al., 2010; Sano et al., 2009; Sinha et al., 2007; Sinha, Berkelhammer, et al., 2011; Sinha, Stott, et al., 2011; Shi et al., 2017, 2018), Australia (McGowan et al., 2012; Palmer et al., 2015), and Europe (Cook, Seager et al., 2015). The recent PAGES 2K consortium is the largest international effort in the history of paleoclimatology to reconstruct the climate history of the common era (Anchukaitis & McKay, 2014; PAGES 2k Consortium, 2017), which has led to the latest continent-scale temperature reconstructions (PAGES 2k Consortium, 2013) and the first generation of paleoclimate reanalysis such as the Last Millennium Reanalysis (LMR; Dee et al., 2016; Hakim et al., 2016; Tardif et al., 2019), Paleo Hydrodynamics Data Assimilation (PHYDA; Steiger et al., 2018), and comparison of multiple methods (Neukom et al., 2019; PAGES 2k Consortium, 2019). Hakim et al. (2016) developed the LMR reanalysis using an “offline” data assimilation method, where static ensemble samples are extracted from existing CMIP5 climate-model simulations to serve as the prior estimate of climate variables. They used linear, univariate forward proxy system models that map climate variables to proxy measurements. They verified the reanalysis with modern instrumental datasets and found that the greatest skill is in the tropics. Tardif et al. (2019) developed LMR version 2 based on the CCSM4 climate model simulations, an updated proxy database and improved seasonal regression-based forward models, which significantly improved surface air temperature and 500 mb geopotential height. Steiger et al. (2018) developed the PHYDA reanalysis using the CESM climate model simulations and a data assimilation method especially aiming at reconstructions of hydroclimate. Neukom et al. (2019) further derived global surface temperature fields for the past 2000 years using six different reconstruction methods.

Significant multi-decadal- to centennial-scale variability of temperature has been found widely in the high-resolution paleoclimate reconstructions and reanalyses. Significant multi-decadal- to centennial-scale variability of temperature was noted in individual paleoclimate records as early as in the 1960s (Dansgaard et al., 1969; Johnsen et al., 1970; Lamb, 1965), and later found in continental mean temperature reconstructions (Gergis et al., 2016; Ljungqvist, Seim, et al., 2019, Ljungqvist, Zhang, et al., 2019; PAGES 2k Consortium, 2013; Werner et al., 2018; Zhang et al., 2018), northern hemisphere mean temperature reconstructions (Cook, Esper et al., 2004; D'Arrigo et al., 2006; Esper et al., 2002, 2012; Guillet et al., 2017; King et al., 2021; Ljungqvist, 2010; Moberg et al., 2005; Schneider et al., 2015; Stoffel et al., 2015; Wilson et al., 2016) and global mean temperature reconstructions (PAGES 2k Consortium, 2019). For example, PAGES 2k Consortium (2019) used seven different statistical and data assimilation methods to reconstruct global mean surface temperature for the past 2000 years based on all available high-resolution paleoclimate records around the world. The seven products consistently show strong multi-decadal- to centennial-scale variability that are coherent with one another (Fig. 11a).

Fig. 11 (A) Global mean surface temperature for the past 2000 years using 7 different reconstruction methods. Data are 30–200-year bandpass-filtered (from PAGES 2k Consortium, 2019). (B) Normalized area-averaged PDSI (z-scores) for Europe (OWDA), North America (NADA), and Asia (MADA) and northern hemisphere (NHDA). Data are low pass–filtered to emphasize variability that was >30 years in duration (from Cook, Seager et al., 2015). (C) Lake level of Lake Bosumtwi in West Africa measured by two paleoclimate proxies: δ¹⁸O (red) and Si intensity (black) (from Shanahan et al., 2009).

Multi-decadal- to centennial-scale mega-droughts have also been found in the high-resolution reconstructions of hydroclimate in the past 2000 years. Centennial-scale mega-droughts were first discovered by Stine (1994). He studied tree ring records from California and Argentine Patagonia, and found that California's Sierra Nevada had extremely severe droughts for more than 200 years between 862 and 1074 (Stine #1 mega-drought) and for more than 140 years between 1122 and 1209 (Stine #2 mega-drought). During these droughts, runoff from the Sierra was much weaker than during any of the modern droughts in that region in the past 140 years. He also found similar drought conditions in Argentine Patagonia. Multi-decadal mega-droughts were first discovered by Stahle et al. (2000). They analysed high-resolution tree-ring data and found significant multi-decadal mega-droughts in different regions of North America in the pre-industrial sixteenth century, as well as in the twentieth century. Drought conditions were prevalent throughout the sixteenth century over the Sierra Madre Occidental, Rocky Mountains, and Mississippi Valley, and also persisted from the 1540s to 1580s in Mexico, from the 1550s to 1590s over the Southwest, and from the 1570s to 1600s over Wyoming and Montana. The droughts of the sixteenth century exceeded the severity, persistence, and coverage of any drought witnessed in the United States during the period of instrumental meteorological or hydrological observation, which is consistent with the results of Woodhouse and Overpeck (1998). Later many studies found strong multi-decadal- to centennial-scale variability of hydroclimate in Africa (Scroxton et al., 2017; Shanahan et al., 2009; Voarintsoa et al., 2017; Wassenburg et al., 2013), Asia (Buckley et al., 2010; Cook, Anchukaitis et al., 2010; Hao et al., 2016; Sankar et al., 2016; Sheppard et al., 2004; Shi et al., 2017, 2018; Sinha, Berkelhammer, et al., 2011; Sinha, Stott, et al., 2011; Tan et al., 2009; Treydte et al., 2006; Wang et al., 2005; Wu et al., 2009; Zhang et al., 2008), Australia (Palmer et al., 2015), Europe (Buntgen et al., 2011; Cook, Seager et al., 2015; Helama et al., 2009; Pauling et al., 2006), North America (Bird et al., 2017; Cook, Woodhouse et al., 2004; Cook et al., 2016; Evans et al., 2018; Fawcett et al., 2011; Kennett et al., 2012; Lachniet et al., 2017; McCabe-Glynn et al., 2013; Medina-Elizalde et al., 2010; O'Mara et al., 2019; Steinman et al., 2012; Thirumalai et al., 2018), South America (Campos et al., 2019; Kanner et al., 2013; Novello et al., 2016), the tropics (Griffiths et al., 2016; Lechleitner et al., 2017), the Northern Hemisphere (Asmerom et al., 2013; Ljungqvist et al., 2016), and the Southern Hemisphere (Scroxton et al., 2017). Cook, Seager et al. (2015) compared the continental mean PDSI for Europe from the Old World Drought Atlas (OWDA), for North America from the North American Drought Atlas (NADA), for Asia from the Monsoon Asia Drought Atlas (MADA), and for the northern hemisphere from the Northern Hemisphere Drought Atlas (NHDA) (Fig. 11b). The continental mean PDSI time series have been normalized and converted to Z-scores, and are low pass–filtered to emphasize variability that was >30 years in duration. There is not any statistically significant trend in the regional mean PDSI for Europe, North America, Asia or the northern hemisphere, and the PDSI values in the industrial era are not significantly different from the pre-industrial values. On the other hand, All the regional mean PDSI time series show clear multi-decadal to centennial variability. This is also the case for the continental mean PDSI for Oceania (Palmer et al., 2015), which does not show any statistically significant trend, but shows substantial multi-decadal to centennial-scale variability. For Africa, Shanahan et al. (2009) combined geomorphic, isotopic, and geochemical evidence from the sediments of Lake Bosumtwi, Ghana, to reconstruct natural variability in the African monsoon over the past three millennia (Fig. 11c). They found that intervals of severe drought lasting for periods ranging from decades to centuries are characteristic of the monsoon and are linked to natural variations in Atlantic temperatures. Thus the severe drought of recent decades is not anomalous in the context of the past three millennia, indicating that the monsoon is capable of longer and more severe future droughts.

The spatial structure of the mega-droughts has been analysed using the reconstructed drought atlas for Europe (Cook, Seager et al., 2015), North America (Cook, Woodhouse et al., 2004; Stahle et al., 2016), Asia (Cook, Anchukaitis et al., 2010; Hao et al., 2016; Shi et al., 2017, 2018), Australia (Palmer et al., 2015), South America (Christie et al., 2011), and northern hemisphere (Ljungqvist et al., 2016). Cook, Seager et al. (2015) found that Europe was significantly drier in the Medieval Climate Anomaly period than in the Little Ice Age period and the modern period. Cook et al. (2016) presented the horizontal maps of North America for the famous centennial-scale Stine #1 mega-drought between 862 and 1074, and Stine #2 mega-drought between 1122 and 1209 (Fig. 12a), both of which were in fact pan-continent drought occurring during the Medieval Climate Anomaly and are much stronger than the modern droughts. They also showed the horizontal maps for the major multi-decadal mega-droughts between 966 and 985, 1271 and 1297, 1323 and 1350, 1435 and 1478 and 1568 and 1591, all of which centred in the central United States. The 1271–1297 mega-drought coincided with the abandonment of the Ancestral Puebloan native American settlement, while the 1323–1350 mega-drought coincided with the abandonment of the Cahokia native American settlement. Shi et al. (2018) calculated the leading EOFs for Asian summer precipitation over the past 554 years, and found that the EOF1 represents a centennial mode with a pan-continent drought pattern, while the EOF2 represents an interdecadal mode with a tri-pole pattern among Pakistan/India, South China and southeast Siberia (Fig. 12b). Ljungqvist et al. (2016) presented a reconstruction of Northern Hemisphere hydroclimate variability for the past 1200 years based on a network of 196 millennium-long proxy records. The century-mean drought maps show that the northern hemisphere was drier in the ninth century and twelfth century than in the twentieth century. They found that the strong seesaw patterns of moisture observed in modern instrumental records over the western USA, Mediterranean, and China have operated consistently in the past 1200 years. For the tropics and the Southern Hemisphere, Scroxton et al. (2017) compared the centennial-scale variability of hydroclimate in different continents and found that they show good coherency. They suggested that the cause of this coherency could be related to symmetrical changes in continental sensible heating, or to a low frequency sea surface temperature climate mode.

Fig. 12 (A) Horizontal maps of PDSI for the Stine (1994) #1 and #2 centennial-scale mega-droughts in North America (from Cook et al., 2016). (B) The first two EOFs of Asian summer precipitation for the past 554 years (from Shi et al., 2018).

Numerous observational studies found that multi-decadal- to centennial-scale mega-droughts are associated with ocean temperature anomalies and suggest that they are driven by ocean forcing. Mega-droughts in Africa are linked to natural variations in Atlantic ocean temperatures (Shanahan et al., 2009), the NAO (Wassenburg et al., 2013), latitudinal shifts of ITCZ (Voarintsoa et al., 2017), and oscillations of the Indian Ocean dipole (Scroxton et al., 2017). Mega-droughts in Asia are linked to the North Atlantic ice-rifting events – an index of the AMOC (Wang et al., 2005), the AMO (Sankar et al., 2016), and ENSO-like SST pattern in tropical Pacific (Griffiths et al., 2016). Mega-droughts in Australia are linked to the IPO (Palmer et al., 2015). Mega-droughts in Europe coincided with an ENSO-like pattern (Helama et al., 2009). Mega-droughts in North America are linked to the PDO (O'Mara et al., 2019; Steinman et al., 2012), SST anomalies in Kuroshio Extension (McCabe-Glynn et al., 2013), and the AMOC (Thirumalai et al., 2018). Mega-droughts in South America are linked to the zonal SST gradient in the Pacific Ocean (Kanner et al., 2013).

Many studies have found that the spectra of hydroclimate and temperature variability in instrumental and paleoclimate observations follow power-law scaling, in which the amplitude of variability increases at longer time scales (Ault et al., 2013, 2014; Huybers & Curry, 2006; Laepple & Huybers, 2014a, 2014b; Loope et al., 2020; Parsons et al., 2017; Pelletier & Turcotte, 1997). This scaling behaviour is often characterized by the exponent β in the power-law function relating climate variance to timescale. For example, Ault et al. (2013) calculated the spectra of multiple PDSI datasets in western North America, and the composite multi-proxy slope is shown in Fig. 13a. Seftigen et al. (2017) calculated the spectra of the standardized precipitation evapotranspiration index (SPEI) in Scandinavia from the ScandH17 reconstruction (Fig. 13b), while Parsons et al. (2017) calculated the β value for 33 coral datasets in tropical Indian Ocean and tropical Pacific Ocean (Fig. 13c). Reconstructed thermohaline circulation from paleoclimate records also shows significant centennial-scale variability (Ayache et al., 2018).

Fig. 13 (A) North American spectral means for PDSI computed from CMIP5 Last Millennium (LM) simulations comparing with the multi-proxy slope (from Ault et al., 2013). (B) Spectra of the standardized precipitation evapotranspiration index (SPEI) in Scandinavia from CMIP5 LM simulations and CRU/ScandH17 observations. Gray shading shows the 95% confidence interval of the ScandH17 observation (from Seftigen et al., 2017). (C) Ocean temperature spectral slope β at the same locations calculated from coral records, CMIP5 LM simulations, and CMIP5 historical runs (1850–2015 CE) (from Parsons et al., 2017).

Global climate models driven by external forcings have significant difficulties in simulating the observed amplitude, timing and teleconnections of multi-decadal- to centennial-scale mega-droughts, as well as the observed multi-decadal- to centennial variability of ocean temperature. External forcings do have significant multi-decadal- to centennial-scale variability (Otto-Bliesner et al., 2016). Volcano forcing has both multi-decadal variability and centennial-scale variability, while solar forcing has a centennial-scale variability, although the greenhouse gases and land use change do not show multi-decadal- to centennial-scale variability. However, the volcano forcing and solar forcing often excite overly weak multi-decadal- to centennial-scale variability of hydroclimate in climate models (Ault et al., 2013; Loope et al., 2020; Seftigen et al., 2017). The model spectra of hydroclimate indices are generally red noise spectra with much smaller variance than observations at the multi-decadal to centennial time scales (Fig. 12a, b). Global climate models also have difficulties in reproducing the timing and duration of the observed mega-droughts in Africa (Klein et al., 2016), Asia (Fan et al., 2009), and North America (Coats et al., 2013; Coats, Smerdon, Cook, et al., 2015; Coats, Smerdon, Seager, et al., 2015; PAGES Hydro2k Consortium 2017; Parsons et al., 2017). For example, Hydro2k Consortium (2017) compared the JJA PDSI over the American Southwest from the CMIP5/PMIP3 models with the observations. The observed five largest megadroughts all occurred during the Medieval Climate Anomaly (the gray bars), with the first four events not coincided with volcano eruptions. The climate models fail to capture the first four events. Moreover, global climate models produce overly weak teleconnections between mega-droughts and ocean temperature (Coats et al., 2013; Coats, Smerdon, Cook, et al., 2016; Coats, Smerdon, Karnauskas, et al., 2016; Stevenson et al., 2015). The model mega-droughts are stochastically generated by atmospheric internal variability and local land-atmosphere feedback. The models difficulties in simulating the observed amplitude, timing and SST-teleconnections of mega-droughts are likely associated with the models’ overly weak multi-decadal- to centennial-scale variability in ocean temperature (Laepple & Huybers, 2014a; Parsons et al., 2017). Laepple and Huybers (2014) compared 20 CMIP5 last millennium simulations with 33 coral, foraminifer, and alkenone records around the world oceans. They found that simulate SST variability that is systematically smaller than instrumental and proxy-based estimates. Discrepancies in variability are largest at low latitudes and increase with timescale, reaching two orders of magnitude for tropical variability at millennial timescales. Parsons et al. (2017) compared with observations the β value in CMIP5 last millennium and history experiments for 33 coral sites in tropical Indian Ocean and tropical Pacific Ocean (Fig. 12c). The CMIP5 models significantly underestimate the β values in tropical Indian Ocean and tropical central/eastern Pacific Ocean. The overly weak ocean variability makes the continental mega-droughts overly weak and stochastically generated.

Unforced climate model experiments showed that the Atlantic Meridional Overturning Circulation (AMOC), which is part of the global deep ocean circulation – the thermohaline circulation, has intrinsic oscillations at multi-decadal to centennial timescales. Mysak et al. (1993) used a 2D thermohaline ocean circulation model of the Atlantic Ocean, and found that 200–300-year period oscillations exist for a wide range of vertical and horizontal diffusivities. These oscillations are associated with propagating salinity anomalies accompanying the thermohaline circulation. Park and Latif (2008) conducted a multimillenial control experiment of the Kiel Climate Model and found significant multi-decadal and centennial variability of the AMOC (Fig. 14a). The multi-decadal variability is similar to the observed AMO, which originates in the North Atlantic and goes to 3000 meters deep (Fig. 14b), while the multicentennial variability is driven in the Southern Ocean and extends to the ocean floor (Fig. 14c). The longer timescale of the multicentennial variability is likely associated with its larger spatial scale and involvement of deeper ocean volume. Both modes are associated with significant changes of sea ice extent. Delworth and Zeng (2012) conducted a 4000-year control experiment of the GFDL CM2.1 climate model, and found a distinct multicentennial oscillation of surface air temperature over the entire Northern Hemisphere extratropics, but with strongest signals in the Atlantic and Arctic regions. This variability is driven by interhemispheric oscillations in AMOC oceanic heat transport, which is associated with propagation of salinity anomalies along the direction of AMOC. The salinity anomalies provide positive feedback to the AMOC. In short, the overly weak multi-decadal to centennial variability in forced climate model simulations suggest that volcano forcing and solar forcing are inefficient in exciting multi-decadal- to centennial-scale oscillations of the AMOC and the global thermohaline circulation.

Fig. 14 (a) Spectrum of meridional overturning circulation strength. (b) Multi-decadal composite of meridional overturning circulation. (c) Multi-centennial composite of meridional overturning circulation. (from Park & Latif, 2008).

The AMOC and the global thermohaline circulation were discovered by Wust (1935, 1936) and Clowes and Deacon (1935), and its global structure was summarized by Gordon (1986, 1991), Broecker (1987, 1992), Rahmstorf (2006) and Talley (2013). Theoretical studies found that the AMOC and the global thermohaline circulation are driven by two physical mechanisms (see review by Kuhlbrodt et al., 2007): (1) the mixing-driven mechanism (Bjerknes, 1916; Huang, 1999; Jeffreys, 1925; Munk, 1966; Munk & Wunsch, 1998; Oort et al., 1994; Sandstrom, 1908, 1916), and (2) the wind-driven mechanism (Sloyan & Rintoul, 2001; Toggweiler & Samuels, 1993, 1995, 1998; Webb & Suginohara, 2001). Mixing is caused by tides, winds, and waves. Observations in recent decades revealed that tidal mixing is the leading component of ocean mixing, especially in the deep ocean (Egbert & Ray, 2000, 2003; Garrett, 2003; Koch-Larrouy et al., 2010, 2015; Ledwell et al., 2000; Loder & Garrett, 1978; Polzin et al., 1997; Ray & Susanto, 2016; Rudnick et al., 2003; Zaron & Egbert, 2006). Unfortunately, tidal mixing has not been well represented in most ocean models, which generally use an overly large diffusion to drive the AMOC and the global thermohaline circulation (see review by Kuhlbrodt et al., 2007). Right now the global ocean modelling community is installing tidal mixing into the global ocean models to replace the overly large diffusion.

The first generation of paleoclimate reanalysis makes it possible to study the global coherent structure and physical mechanisms of the multi-decadal to centennial mega-droughts. Steiger et al. (2019) used the PHYDA reanalysis, in conjunction with radiative forcing estimates, to demonstrate that megadroughts in the American Southwest were driven by unusually frequent and cold central tropical Pacific SST excursions in conjunction with anomalously warm Atlantic SSTs and a locally positive radiative forcing. Figure 15a shows the maximum entropy spectrum of global mean surface temperature from the LMR reanalysis for the past 2000 years. The spectrum demonstrates strong inter-centennial variability between 100 and 400 years and multi-decadal variability between 40 and 100 years. Similar results are found for the Nino3.4 SST anomaly (Fig. 15b) and north Atlantic SST anomaly (AMO index as defined by Enfield et al., 2001, Fig. 15c) from the LMR reanalysis. An independent deep ocean AMOC index from a paleoclimate ocean sediment record (Mjell et al., 2016) is also analysed, which also shows strong inter-centennial variability and multi-decadal variability, although the inter-centennial spectral peaks are wider than the AMO index from the LMR reanalysis (Fig. 15d). For the external forcings, the IPCC AR6 solar forcing shows significant centennial-scale variability with a spectral peak around 200 years (Fig. 15e), while the IPCC AR6 volcano forcing shows significant multi-decadal and centennial-scale variability (Fig. 15f). Figure 15g shows the spectrum for the tidal gravitation force calculated from the NASA JPL DE431 Planetary and Lunar Ephemerides (https://ssd.jpl.nasa.gov/planets/eph_export.html), which is the total tidal gravitation force above the equator from the Moon, the Sun, Mercury, Venus, Mars, Jupiter and Saturn. The tidal gravitation forces from Uranus and Neptune are negligible since they are too far away from the Earth. The tidal gravitation force demonstrates both multi-decadal and centennial-scale variability (Fig. 15g). In the following, we will examine the multi-decadal mega-droughts and the centennial mega-droughts separately.

Fig. 15 Maximum entropy spectra of (a) the global mean surface temperature, (b) the Nino3.4 SST, and (c) the North Atlantic Ocean SST from the LMR reanalysis, and (d) the AMOC index from Mjell et al. (2016), (e) the solar forcing from IPCC AR6, (f) the volcano forcing from IPCC AR6, and (g) the tidal gravitational force from NASA JPL. All data are for the past 2000 years.

b Multi-Decadal Mega-Droughts

As shown in Fig. 2, the deadliest and costliest droughts in the past 140 years, such as the 1930s Dust Bowl and the 1980s Sahel drought, often occurred during the peak of the multi-decadal dry periods. It has long been speculated that these big droughts and the multi-decadal dry periods are driven by the global SST anomalies associated with the AMO (Folland et al., 1986; Sutton & Hodson, 2005; Zhang & Delworth, 2006), but the modern instrument records are too short to study the multi-decadal mega-droughts. Paleoclimate records showed that the droughts in the modern instrument period are not unprecedent, and there were much stronger droughts in the past 2000 years (Fig. 11). The first generation of paleoclimate reanalysis makes it possible to derive statistically robust results. Following the method of Schubert et al. (2016), we examine the global SST patterns associated with the multi-decadal precipitation anomalies over each continent in the past 2000 years. Figure 16 illustrates the correlation of 40–100-year filtered LMR global SST anomalies with 40–100-year filtered area-averaged precipitation anomalies in the past 2000 years for (a) Asia, (b) Oceania, (c) Africa, (d) North America, (e) Europe, and (f) South America. Magenta lines denote the region for precipitation averaging. Black stars denote the grids with correlation coefficient above the 95% confidence level. The multi-decadal precipitation anomalies in Asia are associated with a negative ENSO-like SST pattern in the Pacific Ocean, a negative IOD pattern in the tropical Indian Ocean, negative SST anomalies in the Southern Indian Ocean, and positive SST anomalies in the North Atlantic Ocean, the eastern Arctic Ocean and the South Atlantic Ocean (Fig. 16a). The multi-decadal precipitation anomalies in Oceania are associated with a negative ENSO-like SST pattern in the Pacific Ocean, a negative IOD-like pattern in the tropical Indian Ocean, positive SST anomalies to the south of Australia, negative SST anomalies in the tropical Atlantic Ocean, and positive SST anomalies in the North Atlantic Ocean and the South Atlantic Ocean (Fig. 16b). The multi-decadal precipitation anomalies in Africa are associated with a positive ENSO-like SST pattern in the Pacific Ocean, a basin-wide warming in the tropical Indian Ocean, positive SST anomalies in the whole Arctic Ocean, the North Atlantic Ocean and the tropical Atlantic Ocean, and negative SST anomalies in the South Atlantic Ocean (Fig. 16c). For the multi-decadal precipitation anomalies in North America, the associated SST patterns (Fig. 16d) are similar to those for Africa (Fig. 16c). For the multi-decadal precipitation anomalies in Europe, the associated SST patterns are also similar to those for Africa but with weaker correlations (Fig. 16e). The multi-decadal precipitation anomalies in South America are also associated with SST patterns similar to those for Africa (Fig. 16f), suggesting that they have the same SST forcing.

Fig. 16 Correlation of the 40–100-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the 40–100-year filtered area-averaged precipitation anomalies for (a) Asia, (b) Oceania, (c) Africa, (d) North America, (e) Europe, and (f) South America. Black stars denote the grids with correlation coefficient above the 95% confidence level. Magenta lines denote the region for precipitation averaging.

These patterns can be compared with the SST patterns associated with the interannual precipitation anomalies over each continent (Figs 5 and 9a). The interannual precipitation anomalies are mainly associated with SST anomalies in the Pacific Ocean, tropical Indian Ocean and sometimes tropical Atlantic Ocean, while the multi-decadal precipitation anomalies are also associated with significant SST anomalies in the entire Atlantic Ocean and the Arctic Ocean. In general, the interannual precipitation anomalies over the different continents are often affected by the ENSO, which is the dominant climate mode of interannual variability of the global climate system, while the multi-decadal precipitation anomalies are often affected by the AMO.

The dominant mode of multi-decadal variability of the global climate system is the AMO (Enfield et al., 2001; Schlesinger & Ramankutty, 1994). Analysis of modern instrument data, which covers two and a half AMO cycles, suggested that the AMO affects droughts in Africa (Folland et al., 1986; Zhang & Delworth, 2006), Asia (Lu et al., 2006; Wang et al., 2009; Zhang & Delworth, 2006), Europe (Sutton & Hodson, 2005), North America (Curtis, 2008; Hu & Feng, 2008; Nigam et al., 2011; Sutton & Hodson, 2005), and South America (Chiessi et al., 2009; Kayano & Capistrano, 2014). The paleoclimate studies reviewed before generally used paleoclimate records from a few stations, from which it was difficult to derive the spatial and global structures. A recent progress was Baek et al. (2017) who used the reconstructed North American Drought Atlas, Monsoon Asia Drought Atlas and Old World Drought Atlas for the past millennium, and found robust teleconnection patterns of the AMO over those three continents.

The first generation of paleoclimate reanalysis makes it possible to examine with statistical confidence the global structures of the AMO lifecycle and the underlying physical mechanisms. Figure 17 shows the lag-correlations of 40–100-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the AMO index for lags (a) −30 years, (b) −15 years, (c) 0 year, and (d) +15 years. The global SST pattern at lag 0 year (Fig. 17c), which is the peak phase of the AMO, is similar to that derived from modern instrument data (e.g. Enfield et al., 2001, their Fig. 1), which is characterized by warm SST anomalies in the North Atlantic Ocean, the tropical Atlantic Ocean, the North Pacific Ocean, the western Pacific Ocean, the southeast Pacific Ocean and the tropical Indian Ocean. Note that the AMO SST pattern is characterized by a region without significant SST anomalies in central/eastern Pacific, which was also found in modern instrument data (e.g. Enfield et al., 2001). At lag +15 years, which is a quarter cycle after the warmest SST anomalies in the North Atlantic Ocean, there is a negative ENSO-like SST pattern in the Pacific Ocean and a negative IOD-like pattern in the tropical Indian Ocean (Fig. 17d). The SST pattern at lag −30 years (Fig. 17a), which is the negative phase of the AMO, mirror that at lag 0 year (Fig. 17c). The SST pattern at lag −15 years (Fig. 17b), which is a quarter cycle before the peak of the AMO, mirror that at lag +15 year (Fig. 17d). Therefore, the ENSO-like pattern is in quadrature with the AMO pattern, and the two patterns correspond to two stages of the AMO lifecycle. This explains why both the North Atlantic SST anomalies (the AMO index; Fig. 15b) and the Nino3.4 SST anomalies (Fig. 15c) show significant multi-decadal oscillations. In fact, the global oceans are linked together in the AMO lifecycle.

Fig. 17 Lag-correlations of the 40–100-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the AMO index for lags (a) −30 years, (b) −15 years, (c) 0 year, and (d) +15 years. Lag-correlations of the 40–100-year filtered global land surface precipitation anomalies in the past 2000 years from the LMR reanalysis with the AMO index for lags (e) −30 years, (f) −15 years, (g) 0 year, and (h) +15 years. Black stars denote the grids with correlation coefficient above the 95% confidence level. For land surface precipitation, the number in the upper-left corner is the fraction of global continent area with correlation coefficient above the 95% confidence level.

The AMO strongly affects the multi-decadal mega-droughts around the world. As shown in Fig. 16, the multi-decadal precipitation anomalies in different continents are affected by the ENSO-like SST pattern and/or the AMO pattern, as well as the IOD-like pattern. Figure 17 also shows the lag-correlations of 40–100-year filtered land surface precipitation anomalies from the LMR reanalysis with the AMO index for lags (e) −30 years, (f) −15 years, (g) 0 year, and (h) +15 years. Black stars denote the grids with correlation coefficient above the 95% confidence level. The number in the upper-left corner is the fraction of global continent area with correlation coefficient above the 95% confidence level. Note the reversed red-blue colour bar with blue colour for wet anomalies and red colour for dry anomalies. At lag 0 year (Fig. 17g), there are wet anomalies over tropical Africa, entire Asia, Alaska, north Canada, Greenland and eastern United States, but dry anomalies over central Europe. 52% of the global continent area are significantly affected by the AMO at this stage, which are associated with the SST patterns in Fig. 17c. At lag +15 years (Fig. 17h), there are dry anomalies over eastern Europe, Middle East, central Asia, United States and Antarctica, but wet anomalies over Australia, eastern Brazil and Canada. 41% of the global continent area are significantly affected by the AMO at this stage when the SST pattern is the negative ENSO-like pattern (Fig. 17d). The precipitation patterns at lag −30 years (Fig. 17e) mirror those at lag 0 year (Fig. 17g). The precipitation patterns at lag −15 years (Fig. 17f) mirror those at lag +15 years (Fig. 17h) but with weaker correlations, which is associated with the weaker SST anomalies at that stage (Fig. 19b).

Figure 18 summarizes of the global impacts of the AMO on (a) the land surface precipitation, and (b) the land surface air temperature in the LMR reanalysis. Shadings show the maximum lag-correlation with the AMO index at each grid. Arrows denote the grids with the maximum lag-correlation above the 95% confidence level, and arrow directions represent the time lag of the maximum correlation with respect to the AMO index. Phase clock is shown on the upper-right corner. The number in the upper-left corner is the fraction of global continent area with the maximum lag-correlation above the 95% confidence level. The AMO strongly affects the land surface precipitation over all the continents, with 88% of the global continent area having significant lag-correlations with the AMO index. The AMO has stronger impacts on global land surface air temperature, with 97% of the global continent area including 100% of the continents outside the Antarctica having significant lag-correlations with the AMO index. Comparison with Fig. 9b shows that the global impacts of the AMO cover a wider range of latitudes than those of the ENSO. The ENSO impacts are mainly over the tropics and subtropics, while the AMO impacts are over all the global continents, which is consistent with the wider latitude range of AMO SST anomalies (Fig. 16). The AMO impacts in mid-latitudes and high-latitudes may not be generated by Rossby wave propagations. Other physical mechanisms such as the response of zonal mean circulations, stratospheric pathways and polar sea ice changes may contribute to the observed AMO impacts in mid-latitudes and high-latitudes.

Fig. 18 Summary of the global impacts of the AMO on (a) the land surface precipitation, and (b) the land surface air temperature in the LMR reanalysis. Shadings show the maximum lag-correlation with the AMO index at each grid. Arrows denote the grids with the maximum lag-correlation above the 95% confidence level, and arrow directions represent the time lag of the maximum correlation with respect to the AMO index. Phase clock is shown on the upper-right corner. The number in the upper-left corner is the fraction of global continent area with the maximum lag-correlation above the 95% confidence level.

What drives the AMO? As reviewed by Lin and Qian (2022), the existing AMO theories can be categorized into seven groups: (1) the stochastic forcing theory (Delworth et al., 1993; Griffies & Tziperman, 1995; Delworth & Greatbatch, 2000), (2) the North Atlantic Oscillation (NAO) coupled oscillator (Timmermann et al., 1998; Latif et al., 2004, 2007; Sun et al., 2015; Delworth et al., 2017; Sutton et al., 2018), (3) the zonal-meridional oscillator (Te Raa & Dijkstra, 2002; Dijkstra et al., 2006), (4) the delayed advective oscillator (Lee & Wang, 2010; Zhang & Wang, 2013), (5) the ice coupled oscillator (Dima & Lohmann, 2007; Frankcombe et al., 2010; Böning et al., 2016), (6) the cloud-radiation feedback theory (Clement et al., 2015, 2016; Bellomo et al., 2016; Brown et al., 2016; Yuan et al., 2016; Cane et al., 2017), and (7) the aerosol forcing theory (Booth et al., 2012; Knudsen et al., 2014; Yuan et al., 2016; J. Wang et al., 2017; Birkel et al., 2018). Theories (1)–(5) emphasize internal dynamics associated with the AMOC, theory (6) emphasizes internal heat fluxes associated with ocean mixed layer, while theory (7) emphasizes external forcing. Using the paleoclimate reanalysis and other observational datasets, Lin and Qian (2022) found that the AMO is associated with significant oscillation of the AMOC flow, which is likely driven by tidal gravitational forcing and enhanced by volcano forcing and cloud-radiation feedback. To examine the relationship between the AMOC and the AMO, Fig. 19 shows the lag-correlations of 40–100-year filtered global SST anomalies with the AMOC index for lags (a) −20 years, and (b) −5 years. The global SST anomalies are significantly correlated with the AMOC index, and the patterns at lag −20 years (Fig. 19a) and lag −5 years (Fig. 19b) are quite similar to those for AMO at lags −30 years (Fig. 17a) and lag −15 years (Fig. 17b), respectively. Therefore, the AMO is robustly linked to the AMOC, with the AMOC leading the AMO by about 10 years.

Fig. 19 Lag-correlations of the 40–100-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the 40–100-year filtered AMOC index for lags (A) −20 years, and (B) −5 years; with the solar forcing for lags (C) 0 year, and (D) +15 years; with the volcano forcing for lags (E) 0 year, and (F) +15 years; and with the tidal gravitational force for lags (G) −15 years, and (H) 0 year. Black stars denote the grids with correlation coefficient above the 95% confidence level.

For the external forcings, Fig. 19 also shows the lag-correlations of 40–100-year filtered global SST anomalies with the solar forcing for lags (c) 0 year, and (d) +15 years; with the volcano forcing for lags (e) 0 year, and (f) +15 years; and with the tidal gravitational force for lags (g) −15 years, and (h) 0 year. The solar forcing does not have any significant correlation with the global SST anomalies (Fig. 19c, d), which is consistent with the lack of multi-decadal variability in solar forcing (Fig. 15e). The volcano forcing does have significant multi-decadal variability (Fig. 15f), but does not show strong correlations with the global SST anomalies (Fig. 19e, f), suggesting that the volcano forcing, which is mainly a radiative forcing, is inefficient in exciting the multi-decadal variability of global ocean circulation. Interestingly, the tidal gravitational force is strongly correlated with the global SST anomalies (Fig. 19g, h). This is amazing since the tidal gravitational force is derived from the NASA JPL Planetary and Lunar Ephemerides, which are totally independent from the paleoclimate records entering the paleoclimate reanalysis, while the CCSM4 climate model used by the LMR reanalysis does not have any tidal forcing. The SST patterns in Fig. 19g, h are quite similar to those in Fig. 19a, b, respectively, and are also similar to those in Fig. 17a, b, respectively, indicating that the tidal gravitational force is in phase with the ENSO-like pattern, leads the AMOC by 5 years, and leads the AMO pattern by 15 years. These results suggest that the tidal gravitational force plays a key role in driving the AMOC and the AMO.

Since the volcano and solar forcings are inefficient in exciting the multi-decadal variability in global climate models, the tidal gravitational force likely provides a new physical mechanism for exciting the multi-decadal to centennial variability. Physically, modulation of Earth’s climate system by the tidal gravitational force is direct since (1) tidal mixing is the leading physical mechanism for driving the global thermohaline circulation (Bjerknes, 1916; Huang, 1999; Jeffreys, 1925; Munk, 1966; Munk & Wunsch, 1998; Oort et al., 1994; Sandstrom, 1908, 1916; see review by Kuhlbrodt et al., 2007), and also drives the upper ocean circulation in the Pacific Ocean (Hasumi et al., 2008; Tanaka et al., 2012; Osafune & Yasuda, 2013; Osafune et al., 2014, 2020; Tatebe et al., 2018); (2) the AMOC has intrinsic multi-decadal to centennial variability as shown by the unforced climate model experiments (Delworth & Zeng, 2012; Mysak et al., 1993; Park & Latif, 2008); and (3) tidal forcing has strong multi-decadal to centennial variability (Fig. 15g). Therefore, oscillations in tidal gravitational force likely cause resonant oscillations of the AMOC and global thermohaline circulation at the multi-decadal to centennial timescales.

In summary, the multi-decadal mega-droughts over each continent are associated with significant global SST anomalies. The dominant mode of multi-decadal variability of the global climate system is the AMO, which has strong impacts on the multi-decadal variability of precipitation over all continents. The AMO is associated with significant oscillation of the AMOC flow, which is likely driven by tidal gravitational forcing and enhanced by volcano forcing and cloud-radiation feedback. There are other possible mechanisms such as the land-atmosphere feedbacks. Currently, the first generation of paleoclimate reanalysis provide a limited set of variables, and it is not possible to analyse the land-atmosphere feedbacks.

c Centennial-Scale Mega-Droughts

The centennial-scale mega-droughts, such as those happening around the Migration Period in Europe (Fig. 1a) and before the continuous wars of the Five Dynasties and Ten Kingdoms period in China (Fig. 1b, denoted as 5D&10K), have strong impacts on human history due to their long duration. Following the method of Schubert et al. (2016), we examine the global SST patterns associated with the centennial-scale mega-droughts over each continent in the past 2000 years. Figure 20 illustrates the correlations of the 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the 100–400-year filtered area-averaged precipitation anomalies for (a) Asia, (b) Oceania, (c) Africa, (d) North America, (e) Europe and (f) South America. Magenta lines denote the regions for precipitation averaging. The centennial-scale precipitation anomalies over all the continents are strongly correlated with global SST anomalies. The centennial-scale precipitation anomalies over Asia are associated with a negative ENSO-like pattern in the Pacific Ocean, and warmings in the North Atlantic Ocean, the eastern Arctic Ocean, the South Atlantic Ocean, the eastern Indian Ocean and the southern Indian Ocean (Fig. 20a). The centennial-scale precipitation anomalies over Oceania are associated with very similar global SST patterns as those for Asia (Fig. 20b). The centennial-scale precipitation anomalies over Africa are associated with a positive AMO-like global SST pattern (Fig. 20c). The centennial-scale precipitation anomalies over North America are associated with a positive ENSO-like pattern in the Pacific Ocean, and warmings in the tropical Atlantic Ocean, the western Arctic Ocean, the tropical Indian Ocean and the southern Indian Ocean (Fig. 20d). The centennial-scale precipitation anomalies over Europe are associated with a global warming pattern (Fig. 20e), bearing some resemblance to the positive AMO-like pattern (Fig. 20c). The centennial-scale precipitation anomalies over South America are associated with a similar global SST pattern (Fig. 20f), suggesting that they have similar SST forcing.

Fig. 20 Correlation of the 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the 100–400-year filtered area-averaged precipitation anomalies for (a) Asia, (b) Oceania, (c) Africa, (d) North America, (e) Europe, and (f) South America. Black stars denote the grids with correlation coefficient above the 95% confidence level. Magenta lines denote the region for PDSI averaging.

Comparing Fig. 20 with Fig. 5 shows that the centennial-scale precipitation anomalies are associated with SST anomalies over larger regions than the interannual precipitation anomalies. Comparison of Fig. 20 and Fig. 16 indicates that the centennial-scale precipitation anomalies and the multi-decadal precipitation anomalies are associated with similar SST patterns for Asia, Oceania and North America, but highly different SST patterns for Africa, Europe and South America, suggesting the existence of different physical mechanisms at different timescales.

The modern instrument records cover two AMO cycles and provide the AMO structures for two events. However, for the centennial-scale variability, the modern instrument records are too short to cover one cycle. Therefore, there has not been any observational evidence concerning the dominant modes of the centennial-scale variability. The first generation of paleoclimate reanalysis makes it possible, for the first time, to examine the dominant modes of the centennial-scale variability. The dominant modes of global SST anomalies at the centennial timescale are shown in Fig. 21, which illustrates the first two EOFs of the 100–400-year filtered global SST anomalies from the LMR reanalysis. EOF1 is an ENSO-like pattern which explains 47.8% of the total variance, while EOF2 is an AMO-like pattern which explains 29.9% of the total variance. Lag-correlation analysis shows that PC1 leads PC2 by about a quarter cycle, suggesting that the two leading modes are associated with a coherent inter-centennial oscillation. We hereby define this oscillation the Global Inter-Centennial Oscillation (GICO), and use PC1 as the GICO index.

Fig. 21 The first two EOFs of 100–400-year filtered global SST anomalies from the LMR reanalysis. EOF1 explains 47.8% of variance, while EOF2 explains 29.9% of variance.

Figure 22 demonstrates the lag-correlation of the 100–400-year filtered SST anomalies in the past 2000 years from the LMR reanalysis with the GICO index for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. At lag −75 years (Fig. 22a), there are significant cold SST anomalies in tropical Pacific, northeast Pacific, southeast Pacific, southwest Pacific, south Indian Ocean, tropical Atlantic and southwest Atlantic. At lag −60 years (Fig. 22b), cold SST anomalies start to fade away in eastern Pacific, but grow in northwest Pacific, which amplify at lag −45 years (Fig. 22c). An El Nino-like SST pattern appears in Pacific Ocean at lag −30 years (Fig. 22d), with significant warm SST anomalies in tropical central and eastern Pacific, northeast Pacific and south Pacific, but is associated with significant cold SST anomalies in north Atlantic and south Atlantic. The pattern amplifies at lag −15 years (Fig. 22e), with warm SST anomalies emerging in tropical Atlantic Ocean and south Indian Ocean. The pattern reaches the peak at lag 0 year (Fig. 22f), with warm SST anomalies covering whole Indian Ocean and tropical Atlantic.

Fig. 22 Lag-correlations of the 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the GICO index for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. Black stars denote the grids with correlation coefficient above the 95% confidence level.

The GICO SST anomalies significantly affect the global land surface precipitation. Figure 23 shows the lag-correlation of the land surface precipitation anomalies in the past 2000 years with the GICO index for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. With the 100–400-year oscillations, the dry/wet periods last for 50–200 years. At lag −75 years (Fig. 23a), there are dry periods in Siberia, but wet anomalies in east China, west Asia and central America, which persist to lag −60 years (Fig. 23b). At lag −45 years (Fig. 23c), the dry anomalies over Siberia still persist, while new dry periods develop in Australia, India and east Africa. At lag −30 years (Fig. 23d), dry anomalies amplify in Australia, India, southeast Asia and Africa, while wet anomalies appear in North America and Argentina. At lag −15 years (Fig. 23e), dry anomalies persist in Australia, India, southeast Asia and Africa, and also appear in Brazil and Alaska, but are replaced by wet anomalies in Russia. The patterns persist to lag 0 year (Fig. 23f).

Fig. 23 Lag-correlations of the 100–400-year filtered land precipitation anomalies in the past 2000 years from the LMR reanalysis with the GICO index for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. Black stars denote the grids with correlation coefficient above the 95% confidence level. The number in the upper-left corner is the fraction of global continent area with correlation coefficient above the 95% confidence level.

Figure 24 summarizes the global impacts of the GICO on (a) the land surface precipitation, and (b) the land surface air temperature in the LMR reanalysis. Shadings show the maximum lag-correlation with the GICO index at each grid. Arrows denote the grids with maximum lag-correlation above the 95% confidence level, and arrow directions represent the time lag of maximum correlation with respect to the GICO index. Phase clock is shown on the upper-right corner. The number in the upper-left corner is the fraction of global continent area with the maximum lag-correlation above the 95% confidence level. The GICO strongly affects the land surface precipitation and surface air temperature around the world. For the land surface precipitation, the GICO is associated with significant dry/wet periods over all the continents except Middle East, central Asia, South Africa, southwestern United States and Mexico (Fig. 24a), which is quite different from the significant ENSO and AMO impacts over these regions (Figs. 9b and 18a). Precipitation anomalies over 67% of the global continent area have significant lag-correlations with the GICO. For the surface air temperature, there are significant lag-correlations over all the continents except the United States and most of Canada, which cover 80% of the global continent area (Fig. 24b).

Fig. 24 Summary of the global impacts of the GICO on (a) the land surface precipitation, and (b) the land surface air temperature in the LMR reanalysis. Shadings show the maximum lag-correlation with the GICO index at each grid. Arrows denote the grids with the maximum lag-correlation above the 95% confidence level, and arrow directions represent the time lag of the maximum correlation with respect to the GICO index. Phase clock is shown on the upper-right corner. The number in the upper-left corner is the fraction of global continent area with the maximum lag-correlation above the 95% confidence level.

What drives the centennial-scale variability is still unclear. The proposed physical mechanisms include ocean internal dynamics, solar forcing, and volcano forcing. Ocean internal dynamics refers to unforced oscillations of the AMOC (Delworth et al., 2007; Delworth & Zeng, 2012; Latif et al., 2013; Martin et al., 2013; Mysak et al., 1993; Park & Latif, 2008; Vellinga & Wu, 2004). Figure 25a shows the correlation of 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the AMOC index at lag +45 years. The AMOC is significantly correlated with a negative ENSO-like pattern in the Pacific Ocean, and warming in the North Atlantic Ocean, the South Atlantic Ocean, the eastern Indian Ocean, and the southern Indian Ocean, suggesting that the GICO is connected to the AMOC. The AMOC is also significantly correlated with the centennial-scale precipitation anomalies over South America, eastern United States, Europe, Africa, east Asia and Australia (Fig. 25d), which are likely through SST-driven teleconnections.

Fig. 25 Correlations of the 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with (a) the AMOC index, (b) the solar forcing, and (c) the volcano forcing. Correlations of the 100–400-year filtered land precipitation anomalies with the above indices are shown in (d) (e) and (f), respectively. Black stars denote the grids with correlation coefficient above the 95% confidence level.

The solar forcing was proposed by many studies to drive the observed centennial-scale variability of climate (Asmerom et al., 2013; Eichler et al., 2009; Novello et al., 2016; Poore et al., 2004; Sagawa et al., 2014; Wang et al., 2005; Waple et al., 2002; Wiles et al., 2004; Yu & Ito, 1999). The solar forcing has strong centennial-scale variability with a sharp spectral peak around 200 years and a secondary spectral peak around 120 years (Fig. 15e). However, Fig. 25b shows the correlations of 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the solar forcing. The correlations are very weak, suggesting that the solar forcing is inefficient in exciting the centennial-scale variability of the global oceans. The correlations of the solar forcing with land surface precipitation anomalies are also weak (Fig. 25e).

The volcano forcing has been proposed by several modelling studies to drive the centennial-scale variability (Crowley & Kim, 1999; Fan et al., 2009; Hartl-Meier et al., 2017; Man et al., 2012). The volcano forcing has strong centennial-scale variability with a sharp spectral peak around 180 years and a secondary spectral peak around 300 years (Fig. 15f). However, Fig. 25c shows the correlations of 100–400-year filtered global SST anomalies with the volcano forcing. Again, the correlations are very weak, suggesting that the volcano forcing is inefficient in exciting the centennial-scale variability of the global oceans. As reviewed before, climate models driven by solar and volcano forcings generally underestimate the variance of the observed centennial-scale variability of the ocean temperature and cannot reproduce the timing of the observed centennial-scale variability. However, the volcano forcing is significantly correlated with the centennial-scale precipitation anomalies over western Europe, tropical Africa, north America and southern South America (Fig. 25f), which are possibly caused by the radiative and/or microphysics effects of volcano aerosols since there isn’t significantly correlated SST anomalies (Fig. 25c). Note that the first generation of paleoclimate reanalysis employed climate model simulations forced by the solar forcing and the volcano forcing, which may contribute to the significant correlations in Fig. 25f. In the future generations of paleoclimate reanalysis, it may be better to use unforced climate model simulations, or provide both forced and unforced versions for comparison.

Contrary to the solar forcing and the volcano forcing, the tidal gravitational force demonstrates strong correlations with the centennial-scale variability of the Earth’s climate system. The tidal gravitational force also has strong centennial-scale variability with four spectral peaks around 105, 130, 190 and 310 years (Fig. 15g). Figure 26 shows the lag-correlations of the 100–400-year filtered SST anomalies in the past 2000 years from the LMR reanalysis with the tidal gravitational force for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. The tidal gravitational force shows significant correlations with the global SST anomalies, and the evolutions of global SST pattern in Fig. 26 are quite similar to the evolutions of global SST pattern associated with the GICO lifecycle in Fig. 22, with the tidal gravitational force leading the GICO by about 15 years (e.g. Fig. 26h for lag +15 years of the tidal gravitational force has similar patterns as Fig. 22h for lag 0 year of the GICO). These results suggest that the tidal gravitational force plays a key role in driving the GICO. Once again, note that the NASA JPL Planetary and Lunar Ephemerides used for calculating the tidal gravitational force are totally independent of the paleoclimate records and climate models used for the paleoclimate reanalysis.

Fig. 26 Lag-correlations of the 100–400-year filtered global SST anomalies in the past 2000 years from the LMR reanalysis with the tidal gravitational forcing for lags (a) −60 years, (b) −45 years, (c) −30 years, (d) −15 years, (e) 0 year, and (f) +15 years. Black stars denote the grids with correlation coefficient above the 95% confidence level.

The tidal gravitational force is also strongly correlated with the centennial-scale precipitation anomalies (Fig. 27). The correlation patterns in Fig. 27 are quite similar to those in Fig. 23, but with even stronger correlations. At different time lags, centennial-scale precipitation anomalies over 40–60% of the global continent area are significantly correlated with the tidal gravitational force, which is likely through SST-driven teleconnections. Figure 28 summarizes the global impacts of tidal gravitational force on the centennial-scale (a) land surface precipitation, and (b) land surface air temperature. The tidal gravitational force strongly affects the centennial-scale precipitation anomalies over all the continents, with 89% of the global continent area having significant lag-correlations with the tidal gravitational force, which is likely through SST-driven teleconnections (Fig. 26). The tidal gravitational force has stronger impacts on centennial-scale land surface air temperature, with 100% of the global continent area having significant lag-correlations with the tidal gravitational force.

Fig. 27 Lag-correlations of the 100–400-year filtered land surface precipitation anomalies in the past 2000 years from the LMR reanalysis with the tidal gravitational force for lags (a) −75 years, (b) −60 years, (c) −45 years, (d) −30 years, (e) −15 years, and (f) 0 year. Black stars denote the grids with correlation coefficient above the 95% confidence level. The number in the upper-left corner is the fraction of global continent area with correlation coefficient above the 95% confidence level.

Fig. 28 Summary of the global impacts of the tidal gravitational force on the centennial-scale 100–400-year filtered (a) land surface precipitation, and (b) land surface air temperature in the LMR reanalysis. Shadings show the maximum lag-correlation with the tidal gravitational force at each grid. Arrows denote the grids with the maximum lag-correlation above the 95% confidence level, and arrow directions represent the time lag of the maximum correlation with respect to the tidal gravitational force. Phase clock is shown on the upper-right corner. The number in the upper-left corner is the fraction of global continent area with the maximum lag-correlation above the 95% confidence level.

In summary, the centennial-scale mega-droughts over each continent are associated with significant global SST anomalies. The dominant mode of centennial-scale variability of the global climate system is the GICO, which has strong impacts on the centennial-scale variability of precipitation over all continents. The solar forcing and volcano forcing are inefficient in exciting the centennial-scale variability of SST, while the tidal gravitational force plays a key role in exciting the GICO. The tidal gravitational force is also strongly correlated with the centennial-scale precipitation anomalies, which is likely through SST-driven teleconnections. There are other possible mechanisms such as the land-atmosphere feedback and the sea ice teleconnections, which need further studies in the future.

4 Summary

This paper reviews the history of research on droughts and mega-droughts. The interannual droughts are often driven by the ENSO, the multi-decadal mega-droughts are often driven by the AMO, while the centennial mega-droughts are often driven by the GICO. Droughts over land are affected by multiple factors such as global sea surface temperature anomalies, local land-atmosphere feedbacks, internal atmospheric variability and external forcings from outside the climate system.

The multi-decadal mega-droughts over each continent are associated with significant global SST anomalies. The dominant mode of multi-decadal variability of the global climate system is the AMO, which has strong impacts on the multi-decadal variability of precipitation over all continents. The AMO is associated with significant oscillation of the AMOC flow, which is likely driven by the tidal gravitational force and enhanced by the volcano forcing and the cloud-radiation feedback.

The centennial-scale mega-droughts over each continent are associated with significant global SST anomalies. The dominant mode of centennial-scale variability of the global climate system is the GICO, which has strong impacts on the centennial-scale variability of precipitation over all continents. The solar forcing and volcano forcing are inefficient in exciting the centennial-scale variability of SST, while the tidal gravitational force plays a key role in exciting the GICO. The tidal gravitational force is also strongly correlated with the centennial-scale precipitation anomalies over all continents, which is likely through SST-driven teleconnections.

Currently, global climate models driven by solar and volcano forcings have significant difficulties in simulating the observed amplitude, timing and teleconnections of multi-decadal- to centennial-scale mega-droughts, as well as the observed multi-decadal- to centennial variability of ocean temperature. Unforced climate model experiments showed that the AMOC, which is part of the global thermohaline circulation, has intrinsic oscillations at multi-decadal to centennial timescales. However, the solar forcing and volcano forcing are inefficient in exciting the multi-decadal- to centennial-scale oscillations of the AMOC and the global thermohaline circulation.

Since the volcano forcing and solar forcing are mainly radiative forcings, increase of SST response would require an increase of climate sensitivity. Right now the state-of-the-art global climate models already have an unrealistically high climate sensitivity, which jumped from 2.0 to 4.7 K in CMIP5 models to 1.8–5.6 K in CMIP6 models (Flynn & Mauritsen, 2020; Meehl et al., 2020). Such a high climate sensitivity is not supported by well-established paleoclimate constraints (Sherwood et al., 2020; Zhu et al., 2020, 2021), and is likely caused by overly large low cloud feedback and aerosol cooling (Flynn & Mauritsen, 2020; Zelinka et al., 2020). Moreover, the solar forcing and volcano forcing are not correlated with the observed multi-decadal and centennial-scale SST anomalies (Figs. 19c–f and 25b–c), which means that even with higher climate sensitivity, they would not be able to reproduce the timing of the observed SST oscillations. Therefore, it is not possible to use the volcano forcing and solar forcing to drive the observed multi-decadal- to centennial-scale variability.

Physically, modulation of Earth’s climate system by the tidal gravitational force is direct since tidal mixing is the leading physical mechanism for driving the global thermohaline circulation, and also affects the upper ocean circulation in the Pacific Ocean. Therefore oscillations in tidal gravitational force likely cause resonant oscillations of the global thermohaline circulation and upper ocean circulation in the Pacific Ocean at the multi-decadal to centennial timescales.

There are other possible mechanisms for multi-decadal and centennial-scale mega-droughts such as the land-atmosphere feedback and the sea ice teleconnections. However, due to the lack of paleoclimate records, the first generation of paleoclimate reanalysis provide a limited set of variables, and it is not possible to analyse the land-atmosphere feedback and sea ice teleconnections. Currently, we have very limited paleoclimate records on soil moisture, dust and paleo-vegetation. Future efforts are needed to significantly expand the archive of paleoclimate soil moisture, dust and vegetation records and incorporate them into the paleoclimate reanalysis.

Regarding the internal processes of the climate system, physics-based evaluations of climate models are important for understanding the relative contributions of SST-driven teleconnections, atmospheric internal dynamics and land-atmosphere feedbacks. Simple comparisons with the observed precipitation could not guarantee the correct partitioning of contributions from the different physical mechanisms. It is important to evaluate step-by-step the physical processes and feedbacks using the corresponding observational data. For example, the SST-driven teleconnections can be divided into three steps: (1) generation of convective heating over the SST anomalies; (2) forcing of atmospheric circulation anomalies and Rossby waves by the convective heating; and (3) generation of precipitation anomalies over the continents by the atmospheric circulation anomalies and Rossby waves. The state-of-the-art climate models can simulate the second step well. However, the first step and the third step both involve the convective parameterization, which is one of the key bottlenecks of global climate modelling. For the first step, model intercomparisons showed that the models have difficulty in simulating the position of precipitation anomalies associated with ENSO SST anomalies, and the model precipitation anomalies over equatorial Pacific shift westward comparing to observations (Bellenger et al., 2014; Wittenberg et al., 2006). When the heating anomalies shift westward, the Rossby wavetrains will be affected. For the third step, currently, the models have difficulty in simulating the diurnal cycle of continental precipitation with the precipitation occurring 6 h earlier than in observations (Covey et al., 2016). Therefore, improving the convective parameterization over the ocean and over the continent is important for improving the SST-driven teleconnection in climate models.

We propose the following future directions for research on droughts and mega-droughts:

1. Use the next generation of Earth system models with tidal forcing in ocean models to study and predict the droughts and mega-droughts. Currently the global ocean modelling community is installing tidal mixing into the global ocean models.

2. Conduct step-by-step evaluations of the physical processes and feedbacks in global climate models for understanding the relative contributions of the SST-driven teleconnections, atmospheric internal dynamics and land-atmosphere feedbacks.

3. Understand the role of internal atmospheric variability in the development of drought, and implications for predictability. What are the physical mechanisms for the Rossby waves associated with extratropical droughts? Are the Rossby waves stationary waves or travelling waves? What force these waves? How do the Rossby waves affect the droughts?

4. Understand the link between heat waves and drought, and how a warming world impacts drought on all time scales.

5. Evaluate the strength of land-atmosphere feedback in climate models using observational data.

6. For paleoclimate reanalysis, use unforced climate model simulations to avoid any possible artificial effects from the prescribed solar forcing and volcano forcing.

7. Expand the paleoclimate reanalysis from the surface to the upper air and the ocean subsurface, and extend the time coverage from the past 2000 years to the past 10,000 years, which will provide the foundation for the research on the structure and physical mechanisms of mega-droughts.

8. Conduct paleoclimate model intercomparison experiments to examine the SST-driven teleconnections for multi-decadal and centennial-scale mega-droughts. Force the AGCMs using the reconstructed SST datasets for the past 2000 years to understand the global teleconnections and impacts of the AMO and GICO. Examine the possible roles of Rossby wave propagations, response of zonal mean circulations, stratospheric pathways and polar sea ice forcing.

9. Examine the land-atmosphere feedback mechanisms for multi-decadal and centennial-scale megadroughts, such as the vegetation feedback (Zeng et al. 1999) and dust aerosol feedback (Forman et al., 2001; Hanson et al., 2010; Miao et al., 2007). For this purpose, we suggest to significantly expand the global archive of paleo-vegetation reconstructions (Buso et al., 2019; Leavitt et al., 2007), and paleo-dust reconstructions (Albani et al., 2015; Albani & Mahowald, 2019; Kohfeld & Harrison, 2001; Maher et al., 2010).

Interannual droughts and multi-decadal- to centennial-scale mega-droughts have profound impacts on human society, and require long-lead predictions in order for the governments to have time for preparation. Currently, the leading global climate forecasting centres are issuing 3-month official forecasts, which provides too short lead time for disaster preparation and mitigation. The tidal gravitational force is highly predictable and the NASA JPL Planetary and Lunar Ephemerides provide predictions for the next 500 years. Therefore, incorporating the tidal gravitational force might significantly extend the lead time of drought predictions.

Drought is the deadliest natural disaster on Earth due to its long duration, wide spatial coverage and direct connection with food supply for human beings. Currently, the different forecasting centres are using different models which sometimes yield quite different predictions, but an international drought prediction bulletin is still lacking. We suggest that an international drought prediction bulletin, which is similar to the current ENSO prediction bulletin, is set up to summarize and inter-compare all the available predictions, and promote international efforts in drought prediction. This will especially benefit the developing countries which do not have the experts for making good drought prediction, and may save many people’s lives.

Acknowledgements

Part of this work was done by TQ when working in the Department of Geography, The Ohio State University. The LMR reanalysis datasets are from https://atmos.washington.edu/~hakim/lmr/. The PHYDA reanalysis datasets are from https://zenodo.org/record/1198817#.YZ5i2S-cZaM. The JPL Planetary and Lunar Ephemerides are from https://ssd.jpl.nasa.gov/planets/eph_export.html. JL designed Fig. 8a and advised Andrea S. Taschetto to make an earlier version of it.

Disclosure statement

No potential conflict of interest was reported by the author(s).