Model and proxy evidence for coordinated changes in the hydroclimate of distant regions over the Last Millennium

The Medieval Climate Anomaly (MCA; ca. 950–1250 CE) and the Little Ice Age (LIA; ca. 1450–1850 CE) were periods generally characterized by respectively higher and lower temperatures in many regions. However, they have also been associated with drier and wetter conditions in areas around the Intertropical Convergence Zone (ITCZ) and the Asian Monsoon region and in areas impacted by large-scale climatic modes like the Northern Annular Mode and Southern Annular Mode (NAM and SAM respectively). To analyze coordinated changes in large-scale hydroclimate patterns and whether similar changes also extend to other periods of the Last Millennium (LM) outside the MCA and the LIA, reconstruction-based products have been analyzed. This includes the collection of tree-ring-based drought atlases (DAs), the Paleo Hydrodynamics Data Assimilation product (PHYDA) and the Last Millennium Reanalysis (LMR). These analyses have shown coherent changes in the hydroclimate of tropical and extratropical regions, such as northern and central South America, East Africa, western North America, western Europe, the Middle East, Southeast Asia, and the Indo-Pacific, during the MCA, the LIA and other periods of the LM. Comparisons with model simulations from the Community Earth System Model – Last Millennium Ensemble (CESM-LME) and phases 5 and 6 of the Coupled Model Intercomparison Project (CMIP5 and CMIP6) show that both external forcing and internal variability contributed to these changes, with the contribution of internal variability being particularly important in the Indo-Pacific basin and that of external forcing in the Atlantic basin. These results may help to identify not only those areas showing coordinated changes, but also those regions more impacted by the internal variability, where forced model simulations would not be expected to successfully reproduce the evolution of past actual hydroclimate changes.

1 Introduction

The climate of the Last Millennium (LM) was characterized in many regions by periods of warmer and cooler conditions with respect to the pre-industrial mean climate (Diaz et al., 2011; Graham et al., 2010), two well-known examples being the Medieval Climate Anomaly (MCA; ca. 950–1250 CE) and the Little Ice Age (LIA; ca. 1450–1850 CE; Graham et al., 2007; Laird et al., 2012; Ledru et al., 2013). Such temperature changes are often associated with changes in external forcing factors like solar variability and volcanic activity (Mann et al., 2009; Schurer et al., 2013; Fernández-Donado et al., 2013). The impact of these changes on hydroclimate is not straightforward, but reconstructions show that episodes of severe drought also took place during the LM in many regions (Cook et al., 2022). Periods of wetter and drier conditions existed in large areas of North America (Cook et al., 2010b), Europe (Luterbacher et al., 2012), South America (Vuille et al., 2012) and Monsoon Asia (Hu et al., 2008), with occurrences that coincide with the MCA and the LIA.

These similarities in the evolution of the hydroclimate of distant regions suggest coordinated changes at a global scale (Graham et al., 2010; Atwood et al., 2021). For example, reconstructions from southwestern North America (Graham et al., 2007; Meko et al., 2001; Cook et al., 2004, 2010b; Hughes and Funkhouser, 1998; Anderson, 2011) and the Mediterranean basin (Luterbacher et al., 2012; Martín-Puertas et al., 2010; Morellón et al., 2009) show drier conditions during the MCA and wetter conditions during the LIA. In contrast, reconstructions from northwestern North America (Steinman et al., 2013; Stevens and Dean, 2008), central Europe (Büntgen et al., 2010, 2011) and the British Isles (Wilson et al., 2012; Proctor et al., 2000) tend to show wetter conditions during the MCA and drier conditions during the LIA.

An antiphased relationship between the MCA and LIA is also observed in tropical areas of South America. Wetter conditions are found during the MCA in areas of northern South America, like the Cariaco Basin in Venezuela (Haug et al., 2001), and drier conditions in areas of central South America, like Peru (Bird et al., 2011; Reuter et al., 2009; Vuille et al., 2012; Thompson et al., 1986; Apaéstegui et al., 2014) and eastern Brazil (Vuille et al., 2012; Taylor, 2010; Novello et al., 2012). In East Africa, Anchukaitis and Tierney (2013) showed coordinated changes for Lake Challa (Verschuren et al., 2009; Wolff et al., 2011), Naivasha (Verschuren et al., 2000; Verschuren, 2001), Tanganyika (Tierney et al., 2010a), Masoko (Gilbert et al., 2002; Garcin et al., 2006, 2007) and Malawi (Johnson et al., 2001; Brown and Johnson, 2005; Johnson and McCave, 2008), with an opposite behavior between the lakes in the south and the north of the Intertropical Convergence Zone (ITCZ). Changes in the ITCZ have also been associated with coordinated changes in the hydroclimate of the Indian Monsoon region. Wetter conditions during the MCA and drier conditions during the LIA are found for wide areas of Pakistan (von Rad et al., 1999; Treydte et al., 2006), India (Sinha et al., 2007; Berkelhammer et al., 2010; Bryson and Swain, 1981; Phadtare and Pant, 2006) and the Arabian Sea (Burns et al., 2002; Fleitmann et al., 2003; Anderson et al., 2002; Gupta et al., 2003). Alterations of the ITCZ have also been associated with changes in the hydroclimate of the Indo-Pacific basin (Atwood et al., 2021). A marked transition between the MCA and LIA is found in areas of China (Wang et al., 2005; Zhang et al., 2008; Hu et al., 2008), Southeast Asia (Buckley et al., 2010, 2007), Indonesia (Griffiths et al., 2016; Newton et al., 2006; Tierney et al., 2010b), northern Australia (Denniston et al., 2016) and the eastern Pacific islands (Sachs et al., 2009; Higley et al., 2018).

Most of the impacted regions are located in the extratropics, in the boundary between low and high pressures within the regions of the Northern Annular Mode and Southern Annular Mode (NAM and SAM; Thompson and Wallace, 2001; Jones et al., 2009; Fogt et al., 2009), or in tropical areas around the ITCZ. Analyses based on simulated data (Roldán-Gómez et al., 2020) link transitions from the MCA to LIA in extratropical areas with alterations in the variability in modes like the NAM and SAM (Ortega et al., 2015; Jones et al., 2009; Fogt et al., 2009), mainly driven by changes in external forcing factors (Roldán-Gómez et al., 2020). In tropical areas, the transitions from the MCA to LIA can be linked to alterations of the ITCZ (Atwood et al., 2021), like a contraction over the western Pacific (Denniston et al., 2016; Yan et al., 2015) and shifts over East Africa (Anchukaitis and Tierney, 2013), the Indian Monsoon region (Fleitmann et al., 2003) and the eastern Pacific (Higley et al., 2018). The analysis of the simulated hydroclimate from the Community Earth System Model – Last Millennium Ensemble (CESM-LME; Otto-Bliesner et al., 2015) performed by Roldán-Gómez et al. (2022) shows that these alterations of the ITCZ may be impacted by both external forcing factors and internal variability.

Even if individual reconstructions from different regions suggest coordinated changes in the hydroclimate during the MCA and LIA, global analyses based on proxy data are scarce and limited to certain reconstructions and products. The analysis of simulated data provides an insight into the mechanisms explaining coordinated changes in the hydroclimate of the LM, but analyses of LM simulations must account for both forced and internal variability (Roldán-Gómez et al., 2022). This work addresses these challenges by performing a comparative assessment of the LM hydroclimate at global and large continental scales from different sources, including reconstructions, model simulations, and hybrid products based on data assimilation and reanalysis. The use of several reconstruction-based products allows the identification with a robustness assessment of those regions more impacted by hydroclimatic changes during the MCA and LIA and more generally during the LM. The comparison between reconstruction-based products and model simulations from different ensembles allows an assessment of the contributions of external forcing and internal variability to the hydroclimate of each region, showing those regions where model simulations can provide representative information about changes in the hydroclimate.

2 Evidence from the literature

To survey the available evidence, we assessed 96 reconstructions reporting changes from wetter to drier or from drier to wetter conditions during the transition from the MCA to LIA (Fig. 1, with detailed references in Table 2). The reconstructions used in this compilation include proxy data from tree rings, marine and lake sediments, speleothems, ice cores, and documentary information. These sources provide information about precipitation, moisture, level of lakes and river flows, which can be linked to drier and wetter conditions. Some reconstructions represent the annual hydroclimate, while some others are sensitive to a particular season, as described in Table 2.

Figure 1 Hydroclimate reconstructions showing (dark brown) drier conditions during MCA and wetter conditions during LIA, (light brown) drier conditions during MCA and wetter conditions during early LIA (until the late 1500s CE), (dark green) wetter conditions during MCA and drier conditions during LIA, (light green) wetter conditions during MCA and drier conditions during early LIA (until the late 1500s CE), and (gray) drier or wetter conditions in MCA and LIA depending on the season. Shown within the map, according to the references included in Table 1, are the current position of the ITCZ in July and January; the ITCZ changes in the transition from MCA to LIA in the Indian Monsoon region (C1.ITCZ), eastern Pacific (C2.ITCZ), East Africa (C3.ITCZ) and western Pacific (C4.ITCZ); the boundary between low and high pressures in the pattern of NAM and SAM; and the changes in this boundary in the transition from MCA to LIA (C.NAM and C.SAM). Arrows only indicate the direction of the changes and not their magnitude. Regional maps over North America and Europe are shown on the right side of the panel. Detailed references are included in Table 2.

Newton et al. (2006)Newton et al. (2006)Fleitmann et al. (2003)Higley et al. (2018)Anchukaitis and Tierney (2013)Denniston et al. (2016)Yan et al. (2015)Li and Wang (2003)Gong and Wang (1999)Ortega et al. (2015)Jones et al. (2009); Fogt et al. (2009)

Table 1 References for the position of ITCZ, the boundary between low and high pressures in the pattern of NAM and SAM, and the changes in these two elements during MCA and LIA according to the codes included in Fig. 1.

Verschuren et al. (2009)Wolff et al. (2011)Verschuren et al. (2000)Verschuren (2001)Tierney et al. (2010a)Gilbert et al. (2002)Garcin et al. (2006, 2007)Johnson et al. (2001)Brown and Johnson (2005)Johnson and McCave (2008)Laird et al. (1996)Graham et al. (2007)Gray et al. (2004)Grissino-Mayer (1995)Meko et al. (2001)Michels et al. (2007)Case and MacDonald (2003)St.George and Nielsen (2002)Benito et al. (2010)Sheffer et al. (2007)Tyson and Lindesay (1992)Benito et al. (2011)Sinha et al. (2007)Berkelhammer et al. (2010)Burns et al. (2002)Bryson and Swain (1981)Phadtare and Pant (2006)Anderson et al. (2002)Gupta et al. (2003)von Rad et al. (1999)Bird et al. (2011)Boucher et al. (2011)Boucher et al. (2011)Buckley et al. (2010)Buckley et al. (2007)Büntgen et al. (2011)Büntgen et al. (2010, 2011)Büntgen et al. (2010, 2011)Cook et al. (2004, 2010b)Proctor et al. (2000)Ely et al. (1993)Esper et al. (2007)Treydte et al. (2006)Griffiths et al. (2016)Helama et al. (2009)Kremenetski et al. (2004)Laird et al. (2012)Ledru et al. (2013)Chepstow-Lusty et al. (2009)Reuter et al. (2009)Luoto et al. (2013)Touchan et al. (2007)Antonioli et al. (2003)Göktürk (2011)Göktürk (2011)Göktürk et al. (2011)Greenbaum et al. (2000)Jones et al. (2006)Woodbridge and Roberts (2011)Kuzucuoglu et al. (2011)Bookman et al. (2004); Neumann et al. (2007)Martín-Puertas et al. (2008)Martín-Puertas et al. (2010)Stine (1994)Hughes and Funkhouser (1998)Llasat et al. (2003)Oswald and Foster (2011)Pederson et al. (2005)Booth et al. (2006)Booth et al. (2006)Magny et al. (2007)Routson et al. (2011)Sachs et al. (2009)Sachs et al. (2009)Higley et al. (2018)Sachs et al. (2009)Sachs et al. (2009)Newton et al. (2006)Tierney et al. (2010b)Langton et al. (2008)Haug et al. (2001)Wang et al. (2005)Stager et al. (2005)Steinman et al. (2013)Steinman et al. (2013)Stevens and Dean (2008)Anderson (2011)Benson et al. (2002)Vuille et al. (2012)Thompson et al. (1986)Vuille et al. (2012)Taylor (2010)Vuille et al. (2012)Novello et al. (2012)Neukom et al. (2010)Apaéstegui et al. (2014)Wilhelm et al. (2012)Moreno et al. (2008)Morellón et al. (2009, 2011)Wilson et al. (2012)Zhang et al. (2008)Denniston et al. (2016)Hu et al. (2008)Moy et al. (2002)Rein et al. (2004)Yan et al. (2011)Kirby et al. (2014)Kirby et al. (2019)Hiner et al. (2016)Kirby et al. (2023)

Table 2 Location, references, season and type of reconstruction (sediments, tree rings, speleothems, documentary, ice cores or multi-proxy) for the reconstructions represented in Fig. 1. The first column indicates the code used in Fig. 1 for each reconstruction.

All the reconstructions included in the compilation cover the MCA and LIA. Some of them show consistent drier or wetter conditions during the whole LIA (dark green and dark brown), while some others show drier or wetter conditions until the late 1500s CE, during the early LIA (light green and light brown), as reported for areas of Pakistan, western India and southern China by Graham et al. (2010).

The reconstructions from Fig. 1 are used as reference for the analysis of reconstruction-based and model-based products included in the following sections.

3 Data and methods

The following reconstruction-based products are employed: the 100-member ensemble of the Paleo Hydrodynamics Data Assimilation product (PHYDA; Steiger et al., 2018b); the 20-member ensemble of the Last Millennium Reanalysis (LMR; Tardif et al., 2019); and the tree-ring-based drought atlases (DAs), including the Old World Drought Atlas (OWDA; Cook et al., 2015), the North American Drought Atlas (NADA; Cook et al., 2010b) and the Monsoon Asia Drought Atlas (MADA; Cook et al., 2010a). Due to the fact that they do not extend back into the MCA and LIA, other DAs have not been included in our analyses, and since the MADA does not extend back into the MCA, it has only been considered for the analysis of time series of Sect. 3.3. Even if the PHYDA and LMR provide data for the whole millennium, some of the proxies used for the assimilation do not extend so far in the past. Figure 2 shows the distribution of proxies from Steiger et al. (2018a), including those that cover the MCA and LIA. There is a large number of proxies covering these periods in regions like North America and Europe, so for these regions the products would be expected to follow the behavior of the proxies. However, for other regions like central Africa and eastern South America, the lack of proxies could increase the contribution of the model used for the assimilation and thus the uncertainty in the products.

Figure 2 Coverage of the MCA and LIA by the proxies from Steiger et al. (2018a), considered for the generation of PHYDA and LMR. Positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1.

Regarding the model-based products, we use the 13 all-forcing simulations of the CESM-LME, together with 10 LM simulations from the Coupled Model Intercomparison Project Phase 5–Paleoclimate Modelling Intercomparison Project Phase 3 (CMIP5–PMIP3; Taylor et al., 2007; Schmidt et al., 2011, 2012; Stocker et al., 2013) and 4 LM simulations from the Coupled Model Intercomparison Project Phase 6–Paleoclimate Modelling Intercomparison Project Phase 4 (CMIP6–PMIP4; Eyring et al., 2016; Jungclaus et al., 2017), all of which are described in Table 3. All these simulations have been interpolated to a common grid resolution of $\mathrm{6}{}^{\circ }\times \mathrm{6}{}^{\circ }$, the coarsest among the analyzed simulations, to allow for the calculation of ensemble averages across the CESM-LME, CMIP5 and CMIP6 products. The CESM-LME and CMIP5 simulations are forced according to the CMIP5–PMIP3 protocol (Schmidt et al., 2011, 2012), while the CMIP6 simulations consider the CMIP6–PMIP4 protocol (Jungclaus et al., 2017). Agreement between reconstruction-based and model-based products would only be expected for those areas dominated by the external forcing, as may be the case for extratropical areas like North America and Europe (Roldán-Gómez et al., 2020) and tropical areas in the Atlantic basin like South America and Central Africa (Roldán-Gómez et al., 2022), but not in those areas dominated by the internal variability, like the tropical Indo-Pacific (Roldán-Gómez et al., 2022).

Cook et al. (2015)Cook et al. (2010b)Cook et al. (2010a)Tardif et al. (2019)Steiger et al. (2018b)Otto-Bliesner et al. (2015)Schmidt et al. (2006, 2014)Watanabe et al. (2011)Schurer et al. (2014)Landrum et al. (2013)Dufresne et al. (2013)Phipps et al. (2012)Yukimoto et al. (2011)Adachi et al. (2013)Giorgetta et al. (2013)Hajima et al. (2020)Yukimoto et al. (2019)Mauritsen et al. (2019)

Table 3 Reconstruction-based and model products considered in this work. When referring to an ensemble of experiments, the abbreviation used for the ensemble and the number of members (N) are provided. The starting year, the available seasons and references describing each product are provided in the last columns.

Considering the wide bibliography of single-point reconstructions showing wetter or drier conditions during MCA and LIA (Fig. 1), most analyses have been based on these two periods. Following the approach of Masson-Delmotte et al. (2013), the MCA has been considered from 950 to 1250 CE and the LIA from 1450 to 1850 CE, even if these temporal intervals might not be suitable for all the regions and for all products (Neukom et al., 2014, 2019). To assess the robustness of our results to the specifics of this definition, we have performed additional analyses based on the trend during the whole period from 950 to 1750 CE and report our findings in Appendix A.

Our analyses focus on reconstructions of the Palmer drought severity index (PDSI; Palmer, 1965). Even if other variables like the soil moisture, the precipitation, or the difference between precipitation and evaporation may also provide information about the hydroclimate conditions, the use of a standardized index like the PDSI is generally preferred when combining different products (Cook et al., 2014). The PDSI takes into account the soil moisture conditions from the previous months, estimated based on temperature and precipitation and considering the water storage capacity of the soil. The cumulative behavior of the PDSI makes its variations smoother than those of atmospheric variables, making it a suitable variable for analyses on decennial to centennial scales, while the direct impact of temperature on the PDSI makes it more sensitive to changes in external forcing than variables exclusively based on precipitation (Roldán-Gómez et al., 2020).

The PDSI is directly provided by the DAs, the PHYDA and the LMR. For the CMIP5, CMIP6 and CESM-LME simulations, PDSI has been computed based on atmospheric variables and soil parameters, following the approach of a self-calibrating PDSI (scPDSI; Wells et al., 2004). The potential evapotranspiration is calculated using the Thornthwaite method (Thornthwaite, 1948). While the Thorntwaite formulation has been shown to have limitations for projections of 21st-century soil moisture conditions because of its strong dependence on temperature, assessments of the Thornthwaite formulation in model simulations of the LM have shown that such limitations are not a concern for the LM period (Smerdon et al., 2015). For the 20th century, the PDSI from the different sources have been compared to the instrumental observation-based PDSI from the Climatic Research Unit (CRU; Barichivich et al., 2021).

The analysis of hydroclimate based on proxy data is sensitive to the seasonality of the proxies. As shown in Table 2, most proxies provide annual values, but certain proxies for the Northern Hemisphere are linked to the hydroclimate of the boreal summer (June, July and August, JJA), while certain proxies in the Southern Hemisphere are associated with the boreal winter (December, January, and February, DJF). The selection of the season considered in the different products is then a critical point for the analyses. As shown in Table 3, the PHYDA provides water-year annual averages and seasonal averages for JJA and DJF. The LMR only provides annual averages, while the Northern Hemisphere DAs exclusively target the JJA period. For all the analyses, JJA values are considered for all the products except for the LMR, for which the annual values are used. To confirm that the annual results from the LMR can be meaningfully compared to the seasonal results from the other products, the annual PDSI from the PHYDA and the ensembles of CMIP5, CMIP6 and CESM-LME are analyzed in Appendix B.

Figure 3 LIA–MCA differences in PDSI (shading) or precipitation (contours) from (a) DAs, (b) PHYDA and (c) LMR. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the precipitation, contours of −1 mm per month (brown) and 1 mm per month (green) are shown. The average of JJA is considered for the PHYDA and the DAs, while the annual average is considered for the LMR. For the PHYDA and the LMR, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Coordinated changes during the MCA and the LIA are characterized by computing the difference between each variable and product during the LIA and MCA periods. The agreement between products is also assessed by computing the global spatial correlation between PDSI LIA–MCA maps from each product (Coats et al., 2013a) and by counting the number of products showing positive or negative PDSI differences for each location. Finally, time series of PDSI for the whole millennium are extracted from different locations and different products to analyze the correlations between them beyond the periods of the MCA and LIA.

4 Results

4.1 Reconstructions of hydroclimate during the MCA and LIA

PDSI differences between the LIA and MCA are shown for the DAs, the PHYDA and the LMR (Fig. 3). Some wetting and drying patterns are consistent across the DAs (Fig. 3a), the PHYDA (Fig. 3b) and the LMR (Fig. 3c). For example, all products estimate drier conditions during the LIA and wetter conditions during the MCA over the Middle East. There is also good agreement across the PHYDA and LMR in South America, with drier conditions in the north and wetter conditions in the central continent during the LIA, as well as in northern Canada, Scandinavia and northern Eurasia, with wetter conditions during the LIA. The PHYDA and the LMR agree well with the reconstructions of Fig. 1, which also show drier conditions during the LIA in southeastern Europe and Anatolia (Göktürk, 2011; Jones et al., 2006; Woodbridge and Roberts, 2011; Kuzucuoglu et al., 2011), the Indian Monsoon region (Burns et al., 2002; Bryson and Swain, 1981; Phadtare and Pant, 2006; Anderson et al., 2002; Gupta et al., 2003), and South America north of the ITCZ (Haug et al., 2001). The LMR also agrees with the reconstructions of Fig. 1 in regions of East Africa (Anchukaitis and Tierney, 2013).

In western North America there is good agreement between the DA and PHYDA and between these two products and the reconstructions of Fig. 1, most of them showing wetter conditions in the southwest and drier conditions in the northwest during the LIA (Steinman et al., 2013; Cook et al., 2004, 2010b). Each of these products is nevertheless in contrast to the LMR, which shows drier conditions during the LIA for most of the American west, in sharp contrast with a broad range of literature on megadroughts during the MCA period (Cook et al., 2010b; Coats et al., 2013b). For the case of central and western Europe, the PHYDA estimates drier conditions during the LIA, in agreement with the reconstructions of Fig. 1, while the DA and LMR estimate wetter conditions. The same happens for areas of central and eastern North America, especially for those areas with a small number of proxies covering the MCA (Fig. 2). Differences between products are particularly relevant in the Indo-Pacific basin, where the PHYDA and LMR show opposite behavior in Australia, Southeast Asia and Patagonia, and at high latitudes, with opposite behavior in areas of Alaska and northeast Asia.

4.2 Simulated hydroclimate during the MCA and LIA

The analysis of the model ensembles allows an estimate of those areas that are potentially more impacted by external forcing as estimated by the models because individual simulations include different realizations of internal variability but share the external forcing changes imposed as boundary conditions (Fernández-Donado et al., 2013; Jungclaus et al., 2017).

Figure 4 shows the mean PDSI differences between the LIA and MCA for the CESM-LME, CMIP5 and CMIP6 LM ensembles. For the areas of the Middle East and southern Asia the agreement between the DAs, the PHYDA, the LMR and the reconstructions of Fig. 1 also extends to the simulations from CESM-LME (Fig. 4a), CMIP5 (Fig. 4b) and CMIP6 (Fig. 4c). There is also good agreement in central South America, where the CESM-LME (Fig. 4a), CMIP5 (Fig. 4b) and CMIP6 (Fig. 4c) simulations all show a wet LIA for the areas south of the ITCZ. For northern South America, the dry LIA observed in reconstructions (Fig. 3b and c) is also found, with a limited extension, in CMIP5 (Fig. 4b) and CMIP6 (Fig. 4c) simulations. The agreement found for these areas is consistent with an impact of external forcing on the position of the Atlantic ITCZ (Roldán-Gómez et al., 2022) and the intensity of the monsoon system (Roldán-Gómez et al., 2020).

Figure 4 LIA–MCA differences in PDSI (shading) or precipitation (contours) from the ensemble average of (a) CESM-LME, (b) CMIP5 and (c) CMIP6 LM simulations. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the precipitation, contours of −1 mm per month (brown) and 1 mm per month (green) are shown. In all the cases, the average of JJA is considered. For each product, regional maps over North America and Europe are shown on the right side of each corresponding panel.

The CMIP5 simulations (Fig. 4b) show a wetter LIA in most areas of western North America and western Europe, which is also found with a more limited extension in the reconstructions of Fig. 1 and the PHYDA (Fig. 3b). However, this consistency is not so clear in the simulations from CMIP6 (Fig. 4c), for which a dry LIA is found in most regions of North America. Disagreement also exists in the CESM-LME (Fig. 4a), which simulates a dry LIA in some areas of southwestern North America and a wet LIA in western and northern Europe. As shown in Appendix C, this disagreement can be explained by a different behavior in terms of temperatures and atmospheric dynamics.

In the Indo-Pacific basin, the ensemble average of CESM-LME (Fig. 4a), CMIP5 (Fig. 4b) and CMIP6 (Fig. 4c) LM simulations shows small PDSI differences between the MCA and LIA, with no clear agreement with the PHYDA (Fig. 3b) and the LMR (Fig. 3c). This behavior is consistent with a larger impact of internal variability in these areas (Roldán-Gómez et al., 2022), with simulation-dependent changes that are filtered out when working with ensemble averages.

4.3 Agreement between products

To assess the agreement between different pairs of products, spatial pattern correlations for the global maps of Figs. 3 and 4 have been included in Table 4. Despite the regional differences found in areas of North America and the Indo-Pacific basin, the ensembles of CESM-LME, CMIP5 and CMIP6 show global spatial correlations between 0.29 and 0.39. The global correlation of reconstruction-based products is strongly impacted by low correlations in some areas of the Southern Hemisphere and at high latitudes, with a lower density of proxy data. If these regions are excluded (Table 5), correlations larger than 0.29 are also found between DAs and PHYDA and between PHYDA and LMR. The disagreement between LMR and the other reconstruction-based products in areas of western Europe and western North America explains the small global correlation between the DAs and LMR. It is important to note that the spatial correlation between reconstruction-based and model-based products is impacted by internal variability. The use of ensemble averages removes the contribution of internal variability as it was manifest in the climate trajectory from the simulated data. For the reconstruction-based products the contribution of internal variability is not removed, even if the multi-century averages for the LIA and MCA periods emphasize the contribution of the external forcing.

Table 4 Correlation between the maps of LIA–MCA in PDSI from DAs (Fig. 3a), PHYDA (Fig. 3b), LMR (Fig. 3c), CESM-LME (Fig. 4a), CMIP5 (Fig. 4b) and CMIP6 (Fig. 4c). Significant correlations (p<0.05) accounting for autocorrelation are shown in bold.

Table 5 Same as Table 4 but considering only the data between 20^∘ S and 50^∘ N. Significant correlations (p<0.05) accounting for autocorrelation are shown in bold.

The agreement between sources is particularly poor in areas of the Indo-Pacific basin. As shown in Fig. 5a, there is no clear agreement regarding the sign of the PDSI LIA–MCA differences for areas of Southeast Asia and Australia. There is also no clear agreement for these areas among the individual simulations comprising the CMIP5 and CMIP6 ensembles (Fig. 5b), suggesting that the hydroclimate of these areas is strongly impacted by simulation-dependent processes, likely those associated with internal variability. The link to the internal variability is clear when analyzing the simulations of the CESM-LME (Fig. 5c), which share the same forcing conditions and differ only in the internal variability. These results are in line with the results obtained with the CESM-LME by Roldán-Gómez et al. (2022), showing a relevant impact of internal variability on the position and extension of the Indian and western Pacific ITCZ.

Figure 5 (a) Number of products, from the six analyzed products (DAs, PHYDA and LMR as well as the ensemble averages of CESM-LME, CMIP5 and CMIP6), showing positive (green) or negative (brown) LIA–MCA differences in PDSI. (b) Number of simulations, from the 14 simulations in the ensembles of CMIP5 and CMIP6, showing positive (green) or negative (brown) LIA–MCA differences in PDSI. (c) Number of simulations, from the 13 simulations of the CESM-LME, showing positive (green) or negative (brown) LIA–MCA differences in PDSI. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. Regional maps over North America and Europe are shown on the right side of each corresponding panel.

In other areas, most products (Fig. 5a) and most individual simulations (Fig. 5b and c) show the same sign of PDSI differences such as in central South America, with a drier MCA and a wetter LIA, or the Sahel, extending from western Africa to the Arabian Peninsula, for which most products and simulations show a wetter MCA and a drier LIA. Consistency also extends to areas of western North America, with at least four products showing wetter conditions in the south and drier conditions in the north during the LIA, or to northern India, where most products show a wetter MCA and a drier LIA. The agreement for these areas suggests a contribution of external forcing that is more important than that of internal variability.

Figure 6 (a–e) Time series of PDSI for (a) southwestern North America, (b) northwestern North America, (c) southwestern Europe, (d) Pakistan and (e) East Africa, obtained from DAs, PHYDA and LMR as well as the ensemble average of CESM-LME, CMIP5 and CMIP6 LM simulations. Time series from the CRU PDSI data set for the 20th Century have also been included for comparison. Time series of DAs are not available for those regions not considered in OWDA, NADA and MADA. Periods of MCA and LIA are highlighted in red and blue in the middle panels. Dots in the map of the left side represent the locations of Fig. 1, with the same color code considered in that figure. (f–j) Correlations between time series of PDSI for each region obtained from DAs, PHYDA and LMR as well as the ensemble average of CESM-LME, CMIP5 and CMIP6 LM simulations for both the full millennium (LM) and the pre-industrial period (PRE). Color code for the correlations is the same as included in Fig. 8.

The MCA and the LIA are periods extensively addressed in the literature, but they are not the only periods showing relevant droughts (Cook et al., 2022). The analyses have then been extended to other periods by extracting time series of PDSI for the whole millennium from the same regions that show coordinated changes during the MCA and LIA (Figs. 6a–e and 7a–e), including western North America, southwestern Europe, northern and central South America, East Africa, Pakistan, India, Southeast Asia, and Indonesia. For each region, the correlations between the DAs, the PHYDA and the LMR as well as the ensemble average of CESM-LME, CMIP5 and CMIP6 LM simulations have been computed (Figs. 6f–j and 7f–j).

Figure 7 Same as Fig. 6 for (a, f) Indonesia, (b, g) India, (c, h) Southeast Asia, (d, i) northern South America and (e, j) central South America. Time series of DAs are not available for those regions not considered in OWDA, NADA, MADA and MXDA. Time series of CMIP6 are not available for Indonesia, since most CMIP6 simulations do not provide the necessary variables to compute the PDSI for this region.

The correlation between the DAs and PHYDA is significant for western North America, while the correlation between the PHYDA and LMR is significant for central South America, northwestern North America, southwestern Europe, Pakistan, Southeast Asia and Indonesia. Some of these correlations exceed 0.7 and are clearly visible in the time series of PDSI from these regions. This is the case for the correlation between the DA and the PHYDA in southwestern North America (Fig. 6a and f) and the correlation between PHYDA and LMR in Pakistan (Fig. 6d and i), Indonesia (Fig. 7a and f) and central South America (Fig. 7e and j). Regions showing high correlations between reconstruction-based products are in general those regions showing the largest hydroclimate changes during the LM, like those around the ITCZ (Yan et al., 2015; Anchukaitis and Tierney, 2013), those in the Monsoon region (Fleitmann et al., 2003) and those in the area of influence of the NAM (Ortega et al., 2015).

Figure 8 Correlations between time series of PDSI for different regions (solid lines), obtained from (a) PHYDA, (b) LMR and (c) CESM-LME. Only significant correlations (p<0.05, $\left|r\right|>\mathrm{0.34}$) accounting for autocorrelation are shown, following the color code on the bottom left side of the panels. The color of the points for each region indicate whether the LIA was wetter (green) or drier (brown) than the MCA in the reconstructions of Fig. 1.

Significant correlations are also found between the ensemble averages of CESM-LME, CMIP5 and CMIP6 for regions of South America, India and Southeast Asia; between CMIP5 and CMIP6 ensembles for areas of southeastern Europe, East Africa and Pakistan; and between CESM-LME and CMIP5 ensembles in areas of southwestern Europe and Indonesia. As for the case of reconstruction-based products, some of these correlations exceed 0.7 and are visible in the time series of PDSI, like those of northern South America (Fig. 7d and i) and Southeast Asia (Fig. 7c and h) for CESM-LME, CMIP5 and CMIP6 and those of East Africa (Fig. 6e and j) and India (Fig. 7b and g) for CMIP5 and CMIP6. In these cases, similarities are found during the whole millennium, including a trend in the 20th century not found in the reconstructions (Ljungqvist et al., 2016). If only the pre-industrial period is considered, the correlations are smaller but still significant for most of these areas. The correlation between simulated products is particularly high for those areas with a larger impact of external forcing, like the Atlantic ITCZ (Roldán-Gómez et al., 2020), consistent with the fact that the use of the ensemble average emphasizes the forcing signal (Fernández-Donado et al., 2013; Roldán-Gómez et al., 2020).

Considering this, the comparison between time series from reconstruction-based products and from model simulations would only be relevant for those areas where the contribution of the forcing dominates. Significant correlations are found between the ensemble average of CESM-LME, CMIP5 and CMIP6 LM simulations and the PHYDA and LMR in northern and central South America and, to a lesser extent, in regions of southwestern Europe, northern North America and Southeast Asia, suggesting a relevant contribution of the forcing to the hydroclimate of these areas. Some of these correlations between reconstruction-based products and model simulations exceed 0.5 and are visible in the time series of PDSI, like those in central South America (Fig. 7e and j) and those between the LMR and the ensemble average of CMIP5 in southwestern Europe (Fig. 6c and h). These results can then be linked to a particularly important contribution of the forcing in these regions, present in both proxy data and model simulations.

4.4 Correlation between distant regions

As described in the previous sections, distant regions showed similar hydroclimate behaviors during the MCA and LIA. These similarities suggest the presence of coordinated changes between distant regions that go beyond these two periods. To analyze these changes, correlations between the time series of different regions from Figs. 6 and 7 have been computed for the PHYDA (Fig. 8a) and the LMR (Fig. 8b). Both reconstruction-based products show significant positive and negative correlations between distant regions, some of them in agreement with the behavior of reconstructions during the MCA and LIA (Fig. 1).

In particular, both the PHYDA and the LMR show significant negative correlations between central South America and the regions of East Africa and Pakistan. These correlations can be linked to coordinated changes in the Atlantic ITCZ and the monsoon system, in line with those observed in model simulations and associated with shifts in the Atlantic ITCZ in response to changes in the external forcing (Roldán-Gómez et al., 2020, 2022). Both reconstruction-based products also show significant positive correlations across the Indian and Pacific basins, such as those between Indonesia and East Africa and between Southeast Asia and northwestern North America, consistent with the changes in the Pacific Walker circulation found in simulations (Yan et al., 2015).

The correlations across the Pacific basin are emphasized in the PHYDA, which also shows significant negative correlations between southwestern North America and the regions in India and Indonesia. In contrast, the LMR tends to emphasize the variability in the Atlantic basin, with positive correlations between southwestern North America and southwestern Europe following the pattern of the North Atlantic Oscillation (NAO).

The correlations for the PHYDA and LMR shown in Fig. 8a and b are impacted by the climatological covariance of the model used for the data assimilation. However, if the correlations are computed for the CESM-LME (Fig. 8c) without any data assimilation, only positive correlations associated with the impact of temperature in the PDSI are obtained, showing that a certain impact of the proxies may also be present in the results for the PHYDA and LMR. The fact that the PHYDA and LMR respectively emphasize the variability in the Pacific and Atlantic basins also shows that the regional variability given by the proxies is consistent with the spatial patterns from the models. This increases the confidence on the ability of the models to reproduce the behavior of the hydroclimate at a global scale, in particular for those areas and those mechanisms more impacted by the external forcing.

5 Conclusions

The analysis of different reconstruction-based and model-based products shows consistent coordinated changes during the LM in the hydroclimate of distant regions. The areas more impacted by these changes are those around the ITCZ, like northern and central South America, East Africa and the Indo-Pacific; those in the area of influence of the NAM/NAO, like western North America and western Europe; and those of the Indian Monsoon region, extending from the Middle East to Southeast Asia. Even if these changes are particularly important during the MCA and the LIA, for which a high number of reconstructions show drier or wetter conditions, significant correlations between distant regions are also found when considering the whole millennium.

The agreement between reconstruction-based and model-based products for areas of South America and East Africa is indicative of a relevant contribution of the forcing in these regions, consistent with the shifts in the Atlantic ITCZ obtained in simulations in response to the forcing (Roldán-Gómez et al., 2022). Agreement is also high in the Indian Monsoon region, consistent with an alteration of the monsoon system as a consequence of changes in the forcing (Roldán-Gómez et al., 2020). Certain agreement is also found for areas of western North America and western Europe, consistent with an alteration of the NAM/NAO in response to changes in the forcing (Roldán-Gómez et al., 2020), even if the different products show discrepancies at regional scales. These discrepancies may be linked to the contribution of internal variability in these areas, but also to uncertainties in the reconstructed data, limitations in the physics of the models and inaccuracies in the assessment of external forcing considered for the simulations. The largest discrepancies between reconstructed and simulated products are found in the Indo-Pacific basin, suggesting that this area is more impacted by internal variability processes than by changes in the forcing.

The combination of reconstructions, model simulations and hybrid products is a powerful technique to overcome the limitations of each individual source. On one hand, the sparse hydroclimate records available from proxy data can be completed with gridded information of atmospheric dynamics and hydroclimate variables provided by model simulations and hybrid products to extend the analysis to areas not fully covered by the proxies. On the other hand, the representativity of the mechanisms found in simulations, often involving temperatures, atmospheric dynamics and hydroclimate, like the shifts in the ITCZ, the alteration of variability modes like the NAM/NAO, and the changes in the Walker circulation during the MCA and LIA, could be confirmed with the hydroclimate information obtained from proxy data and hybrid products.

In this work, we considered an exhaustive and up-to-date compilation of all the available reconstruction-based and model-based sources to analyze the hydroclimate of the LM at a global scale, with a particular focus on the MCA and LIA as periods with particular forcing conditions. This exercise, perhaps more typical in the analysis of temperatures but systematically applied for the first time to the global hydroclimate, provided a novel assessment of global hydroclimate changes. The focus on the MCA and LIA periods and the comparison with model simulations allowed the isolation of those mechanisms potentially impacted by the forcing, while the novel approach based on combining reconstruction-based and model-based sources showed a coherence in the large-scale changes in the hydroclimate of the LM.

Appendix A: Impact of the definition of the MCA and LIA

For the analyses described in the previous sections, the periods of the MCA and LIA have been respectively defined from 950 to 1250 CE and from 1450 to 1850 CE. This definition may depend on the product and is not applicable at a global scale (Neukom et al., 2014, 2019). For example, reconstructions of the Indo-Pacific basin (Fig. 1, light brown and light green) tend to show stable conditions only during the early LIA until the late 1500s CE rather than during the whole LIA as found in regions of Europe and North America (Fig. 1, dark brown and dark green). To make the analyses independent from the particular selection of periods for the MCA and LIA, the LIA–MCA differences of Figs. 3 and 4 have been replaced in Figs. A1 and A2 by the slope of PDSI computed over the whole period from 950 to 1750 CE. The results obtained with this method show no major differences with respect to those shown in Figs. 3 and 4, confirming that the conclusions obtained in the paper are independent from the particular definition of the MCA and LIA.

Figure A1 Slope of PDSI (shading) or precipitation (contours) from (a) DAs, (b) PHYDA and (c) LMR, obtained with a linear regression for the period 950–1750 CE. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the slope of precipitation, contours of −0.001 mm per month per year (brown) and 0.001 mm per month per year (green) are shown. The average of JJA is considered for the PHYDA and the DAs, while the annual average is considered for the LMR. For the PHYDA and the LMR, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Figure A2 Slope of PDSI (shading) or precipitation (contours) from the ensemble average of (a) CESM-LME, (b) CMIP5 and (c) CMIP6 LM simulations, obtained with a linear regression for the period 950–1750 CE. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the slope of precipitation, contours of −0.001 mm per month per year (brown) and 0.001 mm per month per year (green) are shown. In all the cases, the average of JJA is considered. For each product, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Appendix B: Seasonal changes in the PDSI

In the results presented in the previous sections, the JJA PDSI values from the DAs and the PHYDA as well as the simulations from CESM-LME, CMIP5 and CMIP6 are compared with the annual PDSI provided by the LMR. To confirm that the use of different averaging periods does not have a significant impact on the results, the annual PDSI from the PHYDA and the different ensembles of simulations have been also analyzed. The difference between the LIA and MCA for the annual PDSI of the PHYDA and the ensembles of CESM-LME, CMIP5 and CMIP6 is respectively shown in Figs. B1a, b and B2a, b. When comparing these results with those obtained with the average of JJA (Figs. 3b and 4) only minor differences are observed, mostly limited to regional scales. This shows that despite the relevant seasonal changes found in the atmospheric variables, the cumulative behavior of the PDSI makes it less sensitive to the season.

Figure B1 LIA–MCA differences in PDSI (shading) or precipitation (contours) from (a) PHYDA and (b) the ensemble average of CESM-LME, considering annual averages. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the precipitation, contours of −1 mm per month (brown) and 1 mm per month (green) are shown. For each case, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Figure B2 LIA–MCA differences in PDSI (shading) or precipitation (contours) from the ensemble average of (a) CMIP5 and (b) CMIP6 LM simulations, considering annual averages. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the precipitation, contours of −1 mm per month (brown) and 1 mm per month (green) are shown. For each case, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Appendix C: Link to temperature and atmospheric dynamics

For many regions, the agreement in the PDSI differences between the MCA and LIA from different products could be linked to agreement in temperatures. This is the case for central South America, where all the products show a wet LIA (Figs. 3 and 4), linked to cooler temperatures during this period (Figs. C1 and C2), or the Middle East, where the PHYDA and LMR show a drier LIA and a wetter MCA (Fig. 3), consistent with a cooler MCA and a warmer LIA shown by these products (Fig. C1).

The analysis of temperatures can also explain some of the discrepancies found in the analysis of PDSI. For example, the different behavior obtained from PHYDA and LMR for areas of Australia and Patagonia (Fig. 3) can be explained by an opposite behavior in terms of temperature, with the PHYDA (Fig. C1a) and LMR (Fig. C1b) respectively showing warmer or cooler conditions in Patagonia and cooler or warmer conditions in Australia during the LIA. The drier LIA shown by CMIP6 simulations for regions of North America (Fig. 4c) could also be associated with a limited cooling during that period in the CMIP6 ensemble (Fig. C2c) compared to that shown in CMIP5 (Fig. C2b) and CESM-LME (Fig. C2a), while the wetter LIA shown by CESM-LME simulations in western and northern Europe (Fig. 4a) could be associated with the fact that CESM-LME shows cooler conditions during the LIA for these regions (Fig. C2a) compared to those simulated in CMIP5 (Fig. C2b) and CMIP6 (Fig. C2c) ensembles.

These differences can be linked to differences in atmospheric dynamics, with the CMIP5 ensemble showing important differences in sea level pressure (SLP) between the MCA and the LIA in most regions of North America and Europe (Fig. C2b), while the CMIP6 ensemble (Fig. C2c) and the CESM-LME (Fig. C2a) limit these changes to some areas of southern and eastern Europe and northern North America. Despite these differences, all the models show relevant changes in these regions during the MCA and LIA, in line with the changes in the NAO shown in Roldán-Gómez et al. (2020).

For the Indo-Pacific basin, the CESM-LME (Fig. C2a) tends to show cooler conditions during the LIA for most regions, but the temperature difference is small compared to that of North America and Europe. For the case of CMIP5 (Fig. C2b) and CMIP6 (Fig. C2c), warmer conditions are found during the LIA for some areas of India and northern Australia. Regarding the atmospheric dynamics, only the LMR shows relevant changes in SLP for the Indo-Pacific, with positive anomalies in the Indian and western Pacific basin and negative anomalies in the eastern Pacific basin during the LIA. These relatively small changes found in temperature and SLP could explain the limited impact observed in the analysis of PDSI for this area.

Figure C1 LIA–MCA differences in temperature from (a) PHYDA and (b) LMR. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the LMR, differences in SLP are also shown with contours of −5 Pa (blue) and 5 Pa (red). The average of JJA is considered for the PHYDA, while the annual average is considered for the LMR. For each product, regional maps over North America and Europe are shown on the right side of each corresponding panel.

Figure C2 LIA–MCA differences in temperature (shading) or SLP (contours) from the ensemble average of (a) CESM-LME, (b) CMIP5 and (c) CMIP6 LM simulations. Dots represent the locations of Fig. 1, with the same color code considered in that figure; positions of the modern ITCZ, NAM and SAM boundaries are also shown, as described in Fig. 1. For the SLP, contours of −5 Pa (blue) and 5 Pa (red) are shown. In all the cases, the average of JJA is considered. For each product, regional maps over North America and Europe are shown on the right side of each corresponding panel.