Given these uncertainties, the development and analysis of new large‐scale reconstructions is imperative for improving our understanding of climate dynamics in the Mediterranean. Here we use a spatially resolved tree ring‐based field reconstruction of European and Mediterranean hydroclimate to investigate the dominant spatiotemporal patterns of drought variability across the basin over nearly the entirety of the last millennium. Specifically, our analysis addresses three primary research questions: (1) What are the dominant modes of hydroclimate variability in the Mediterranean? (2) How spatially coherent are drought events across the Mediterranean Basin? (3) And how do recent droughts compare to long‐term hydroclimate variability over the last 900 years?

However, a more complete understanding of natural Mediterranean drought variability and anthropogenically forced moisture trends requires comparisons with long‐term variability that is not available from relatively short instrumental records. To this end, the paleoclimate community has been active throughout this region, developing estimates of Common Era drought variability from a variety of proxies, including tree rings [ Brewer et al. , 2007 ; Chbouki et al. , 1995 ; Glueck and Stockton , 2001 ; Touchan et al. , 2003 ; Akkemik and Aras , 2005 ; Touchan et al. , 2005 ; Esper et al. , 2007 ; Andreu et al. , 2007 ; Nicault et al. , 2008 ; Touchan et al. , 2008a ; Büntgen et al. , 2010 ; Touchan et al. , 2011 ; Köse et al. , 2011 ; Touchan et al. , 2014a ], sediment cores [e.g., Jones et al. , 2006 ; Roberts et al. , 2012 ; Moreno et al. , 2012 ], speleothems [e.g., Jex et al. , 2011 ; Wassenburg et al. , 2013 ], and networks using multiple proxies [ Carro‐Calvo et al. , 2013 ; Luterbacher et al. , 2012 ; Pauling et al. , 2006 ]. To date, however, there is little consensus across these different records regarding the character and dominant drivers of drought variability across the basin over the last millennium. In particular, there are extant uncertainties regarding how widespread droughts are in the Mediterranean [ Roberts et al. , 2012 ], the magnitude and timing of long‐term trends and centennial‐scale variability [ Esper et al. , 2007 ; Touchan et al. , 2011 ; Wassenburg et al. , 2013 ], and how seasonal signals and large‐scale climate modes are reflected in proxy reconstructions [ Touchan et al. , 2014a , 2014b ; Seim et al. , 2015 ].

Climate change impacts on water resources are a significant concern in the regions surrounding the Mediterranean Sea [ Iglesias et al. , 2007 ; García‐Ruiz et al. , 2011 ], an area including southern Europe, northern Africa, and the Levant region of the Middle East (Cyprus, Israel, Jordan, Lebanon, Palestine, Syria, and Turkey). Projections from climate models almost uniformly point toward drying in the Mediterranean from increased greenhouse gas forcing in the coming decades [ Giorgi and Lionello , 2008 ; Collins et al. , 2013 ; Seager et al. , 2014 ], part of an overall trend toward desiccation and poleward expansion of subtropical dry zones [ Held and Soden , 2006 ; Seager et al. , 2010 ]. Indeed, analyses of recent climate trends in the region suggest that this process may have already begun [ García‐Ruiz et al. , 2011 ; Gleick , 2014 ; Hoerling et al. , 2012 ; Kelley et al. , 2012 , 2015 ].

Spectral and spectral coherency analyses are conducted in two ways. The first is a multitaper approach [ Thomson , 1982 ; Chave et al. , 1987 ; Mann and Lees , 1996 ; Czaja and Marshall , 2001 ; Huybers , 2004 ] with significance levels estimated using a nonparametric Monte Carlo procedure with red noise (AR1) conditioned on the original data. We also use wavelet coherency analyses [ Maraun and Kurths , 2004 ; Maraun et al. , 2007 ] to investigate time‐varying coherency and phasing between various drought series. Following the best practice recommendations in Grinsted et al. [ 2004 ], we confirmed normality of the time series used in the wavelet coherency analyses (using Lilliefors test).

We use Spearman's rank correlations to assess the statistical relationship between reconstructions, observations, and indices. Spearman's is a nonparametric alternative to the Pearson product‐moment correlation that is less sensitive to outliers. In all correlation plots, regions with insignificant correlations ( p > 0.05) are masked by grey asterisks. For regional time series, uncertainties are estimated as the 95th percentile bias‐corrected and accelerated (BCa) [ Efron , 1987 ] bootstrap confidence intervals. On the relevant figures, these confidence intervals are shown by the grey‐blue shading around the time series. In some figures, the time series were smoothed using a 10 year loess fit [ Cleveland and Devlin , 1988 ] to emphasize low‐frequency variability, although all statistical analyses were conducted on the original (unsmoothed) data.

To account for the loss of grid cells with declining proxy availability (and regression model degradation) in the past, we use a frozen grid. Any missing grid cells in year 1100 C.E. (the first year of our analysis) were treated as missing in all subsequent years. This ensures stationarity in the number of observations available in all years of our analysis.

Monthly precipitation data are from the high‐resolution CRU gridded climate data (TS3.21) [ Harris et al. , 2014 ]. We also use monthly average climate indices from the Climate Prediction Center for modes that have been previously shown to have an influence on Mediterranean climate [ Sousa et al. , 2011 ]. The North Atlantic Oscillation (NAO) [ Hurrell , 1995 ] consists of a north‐south dipole in atmospheric pressure between Greenland and the subtropical North Atlantic, with positive phases typically associated with below‐average precipitation in the Mediterranean and southern Europe. The Scandinavian Pattern (SCA) [ Bueh and Nakamura , 2007 ] is centered over Scandinavia, with positive height anomalies in this region causing above‐average precipitation across central and southern Europe. The East Atlantic Pattern (EA) [ Barnston and Livezey , 1987 ] is similar to the NAO, consisting of a north‐south anomaly dipole in the Atlantic. Unlike the NAO, however, the EA has stronger connections to the subtropical ridge. Positive phases of the EA are linked to below‐average precipitation across southern Europe. For the CRU precipitation data and the climate indices, we restrict our analysis for the most recent period when data quality and availability is highest (1950–2012).

The scPDSI itself is a modification of the original PDSI formulation of Palmer [ 1965 ], a locally normalized drought index incorporating changes in supply (precipitation), demand (evapotranspiration), and storage (soil moisture). PDSI has an inherent memory timescale of 12–18 months [ Guttman , 1998 ; Vicente‐Serrano et al. , 2010 ], allowing the JJA target of the OWDA to incorporate climate information from the previous winter and spring, the main seasons of moisture supply to the Mediterranean [ Touchan et al. , 2011 ]. The scPDSI used as the target for the OWDA reconstruction [ van der Schrier et al. , 2013 ] is calculated from version TS3.21 of the high‐resolution Climatic Research Unit (CRU) climate grids [ Harris et al. , 2014 ], incorporates a snow module, and uses the Penman‐Monteith formulation for calculating potential evapotranspiration [ Xu and Singh , 2002 ].

The longest available chronologies (start year before 1100 or 1200 C.E.) are well distributed across the Mediterranean region, with greatest densities in Anatolia, the western Mediterranean, and Northern Italy. For regions where these longer chronologies are not available locally (e.g., the Levant and Middle East) the 1000 km search radius (roughly equivalent to the correlation decay e ‐folding distance of the PDSI being reconstructed; see the Supplemental Materials in Cook et al. [ 2015 ]) allows for reconstructions further back in time than the shorter local chronologies would allow. For example, the longest pre‐1100 C.E. chronologies from Anatolia are within 1000 km of the Levant and Middle East region and used in the OWDA to reconstruct scPDSI in these regions prior to the start date of the local chronologies, which mostly begin in the 1300s and 1400s. In such cases, the grid point reconstruction is still required to pass minimum calibration and validation thresholds to be included in the drought atlas. Up to 1978, scPDSI values in the OWDA are from the tree ring reconstruction, which are merged with the instrumental data from 1979 to 2012. Complete details on the OWDA calibration and validation can be found in the Supplemental Materials of Cook et al. [ 2015 ]. Given these caveats, we therefore have confidence that the temporal coverage (1100–2012 C.E.) and spatial domain of our analyses are well supported.

The Old World Drought Atlas (OWDA) [ Cook et al. , 2015 ] is a new, tree ring‐based reconstruction of summer season (June‐July‐August, JJA) self‐calibrating Palmer Drought Severity Index (scPDSI) [ van der Schrier et al. , 2013 ]. The OWDA uses 106 tree ring chronologies to reconstruct scPDSI at 5414 half‐degree grid points for the entire Common Era over the European‐Mediterranean domain (27°N–71°N, 12°W–45°E). The reconstruction uses the point‐by‐point method [ Cook et al. , 1999 , 2013 ] with a proxy search radius of 1000 km around each target grid point. Grid cell scPDSI is reconstructed from the proxy predictor series within the search radius, but the tree ring proxies are effectively weighted by the principal components regression such that those that covary most strongly with observations have the greatest influence on the reconstructions. Proxy site locations in the Mediterranean portion of the OWDA (30°N–47°N, 10°W–45°E) are shown in Figure 1 , along with the approximate start dates of the various records.

3 Results and Discussion

3.1 Climate Signals in the OWDA Correlations between winter (January–March; JFM) and spring (April–June; AMJ) precipitation and the tree ring reconstructed JJA scPDSI are uniformly positive across the basin (Figure 2). The strongest correlations with JFM precipitation are localized in Spain and Morocco in the western Mediterranean and the Levant region in the eastern basin. The AMJ precipitation correlations are more uniform and strongly positive across nearly the entire Mediterranean, suggesting that the summer (JJA) season soil moisture variability, reflected in the OWDA and the underlying tree growth, is driven primarily by spring precipitation [cf. Touchan et al., 2011, 2014a]. Figure 2 Open in figure viewer PowerPoint Point‐by‐point Spearman's rank correlation coefficients between CRU 3.21 precipitation totals (JFM and AMJ) and OWDA summer season (JJA) scPDSI. Correlations are calculated over the period 1950–2012 C.E. Regions of insignificant correlation (p > 0.05) are masked by grey asterisks. Climate mode correlations with the CRU precipitation data or the OWDA scPDSI (Figure 3) are consistent in sign but generally stronger in the precipitation data. The weaker OWDA scPDSI correlations are expected for at least two reasons. First, these climate modes reflect shifts in atmospheric circulation that have a direct impact on precipitation by modulating, for example, storm track positions and moisture convergence. Circulation impacts on the scPDSI will be by definition more indirect, as the scPDSI is computed based on the anomalies of a variety of variables that influence soil moisture. And, as mentioned previously, through year 1978 the scPDSI is reconstructed as a scaled linear function of the underlying tree ring proxies, which imparts additional uncertainties on the estimates of scPDSI. Figure 3 Open in figure viewer PowerPoint Spearman's rank correlation coefficients between the teleconnection indices (NAO, SCA, and EA) and (left column) simultaneous season CRU 3.21 precipitation totals and (right column) OWDA summer season scPDSI. Correlations are calculated over the period 1950–2012 C.E. Regions of insignificant correlation (p > 0.05) are masked by grey asterisks. The influence of the NAO is strongest in winter and early spring (January–April). Positive phases of the NAO cause drying across the northern reaches of the Mediterranean Basin, from Spain and Morocco across to the Balkans and western Turkey, while favoring wetter conditions in coastal regions of Libya, Egypt, and the Levant. Consistent with both the instrumental observations [e.g., Lamb et al., 1997; Knippertz et al., 2003] and the underlying controls on tree growth [Touchan et al., 2008a, 2008b, 2011; Panayotov et al., 2010], the influence of the NAO in the OWDA is strongest in the far western part of the domain and the Balkans. The largest influence of the SCA pattern is during the winter (JFM), with positive phases increasing moisture across most of the basin. The SCA correlation with the OWDA is substantially weaker, consistent with previous analyses (Figure 2) demonstrating the much stronger connection between the OWDA and spring, rather than winter, precipitation. Unlike the previous two modes, the influence of the EA is highest during the spring (AMJ), causing widespread drying across the basin with, notably, similar magnitude correlations with both precipitation and scPDSI. The results from our teleconnection analyses are similar to previous studies quantifying connections between these climate modes and drought variability in the region. For example, Sousa et al. [2011] found comparable patterns of significant correlation between PDSI and the NAO. Their correlations with the SCA and EA are similar in sign, though with larger magnitudes and greater significance than in our analyses, possibly because they focused on different seasonal composites. Other comparisons between the NAO and precipitation over the Mediterranean are also similar [Cullen and deMenocal, 2000; Roberts et al., 2012; Xoplaki et al., 2004], showing a broad pattern of significant negative correlations across Southern Europe and out of phase anomalies over coastal Africa in the eastern part of the basin.

3.2 Drought Variability in the OWDA Figure 4 shows OWDA scPDSI averaged over the Mediterranean domain and the fraction of the land area in drought conditions (scPDSI≤−1) each year from 1100 to 2012. Interannual variability (standard deviation) in the scPDSI is 0.54 standardized PDSI units, and, on average, 29% of the basin experiences drought conditions in any given year. Several periods of persistent, pan‐Mediterranean drought are apparent in the record, including in the 1100s, 1200s, and 1300s. There is also a particularly strong pluvial event in the early 1100s. Noticeably absent, however, are any extended multidecadal (30 year or longer) “megadroughts,” a characteristic feature of North American drought variability during the Medieval Climate Anomaly (approximately before 1300 C.E.) [e.g., Cook et al., 2010a] and previously inferred from Moroccan tree rings [Esper et al., 2007]. Figure 4 Open in figure viewer PowerPoint (top) Area average scPDSI for the entire Mediterranean domain in the OWDA (30°N–47°N, 10°W–45°E) and (bottom) percent land area in drought (scPDSI≤−1) from 1100 to 2012 C.E. Red curves are smoothed versions of the time series using a 10 year loess smooth. In Figure 4 (top), regional average scPDSI calculated from the instrumental target data set for the OWDA reconstruction is overlain in orange. The horizontal line in Figure 4 (bottom) is the long‐term average fractional area in drought from 1100 to 2012 C.E. (29%). Grey‐blue shading indicates the 95th confidence intervals, estimated using a BCa bootstrap. The spatial patterns of five periods of widespread drought (and one pluvial) are shown in Figure 5. A nearly two decade long pluvial occurred from 1125 to 1142 C.E., with sustained wet conditions in modern day Spain, Morocco, Algeria, Tunisia, Italy, the Balkans, and Turkey. This event occurred at the same time as an extended period (1118–1179 C.E.) of low flows in the Upper Colorado River Basin in North America (the most persistent dry period in that region over the last millennium) [Meko et al., 2007], severe drought in the Sacramento River Basin [Meko et al., 2012], and widespread drought across most of North America [Cook et al., 2014a]. This synchrony is suggestive of a possible large‐scale persistent anomaly in atmospheric circulation across the Atlantic, possibly related to shifts in Atlantic sea surface temperatures [Hu and Feng, 2012]. For the moment, however, such a relationship is speculative, requiring additional data to determine whether this co‐occurrence is random or reflective of some real underlying dynamics. With Common Era drought atlases now available across all three continents of the Northern Hemisphere [Cook et al., 2010a, 2010b, 2015], there may be new opportunities to investigate such a question in more detail. Figure 5 Open in figure viewer PowerPoint Multiyear average scPDSI for different pan‐basin drought and pluvial events in the OWDA. Among the five highlighted persistent drought events (Figure 5), there is a tendency for simultaneous drought in the extreme western (Spain, Morocco, Algeria, and Tunisia) and eastern (Balkans, Greece, and Turkey) ends of the basin. This suggests a degree of spatial coherency and synchrony during major drought events in the basin. This is further confirmed by compositing the most widespread drought years in the Mediterranean, years when scPDSI values ≤−1.0 cover at least 40% or 50% of the basin (Figure 6). Severe droughts occur synchronously at the opposite ends of the basin. Additionally, Figures 5 and 6 also show evidence for antiphasing behavior in the eastern end of the basin, where there is a tendency for Libya, Egypt, and the southern Levant to be wet or near normal when Greece and Anatolia are in drought. This meridional dipole structure bears some resemblance to the NAO correlation patterns with precipitation and scPDSI discussed previously (Figure 3). Figure 6 Open in figure viewer PowerPoint Composite average scPDSI for years in the OWDA with drought area (DA; scPDSI≤−1) exceeding 40% (n = 168 years) and 50% (n = 45 years) of the total land area in the Mediterranean domain. These spatial patterns of variability have been documented in other studies of Mediterranean drought variability. For example, Xoplaki et al. [2004] previously noted north‐south antiphasing in hydroclimate variability over the eastern Mediterranean and a tendency for zonally widespread drought events across the basin. These observations of spatial coherence are, however, contrary to previous analyses of other proxy records in the region, which have suggested instead a tendency for an antiphased dipole in hydroclimate variability between these two regions [Roberts et al., 2012]. Whether these differences are due to intrinsic qualities of the proxies themselves or differing interpretations of the underlying data is still unclear. Independent reconstructions of Mediterranean climate over the Common Era are available from a variety of proxies, including lake sediments, speleothems, and historical records [Luterbacher et al., 2012]. Comparisons between these records and the OWDA, however, are difficult for several reasons, including temporal resolution and time uncertainty, spatial coverage, frequency biases, the reconstruction methodologies used, and the variables targeted for reconstruction (e.g., lake levels, precipitation, and raw isotope time series). The review by Luterbacher et al. [2012] highlights several centennial‐scale periods of enhanced aridity or wetness in the Mediterranean during the Common Era from various different proxy reconstructions. Some of these records are at least qualitatively comparable to the OWDA. For example, the carbon isotope record from Uzuntarla Cave [Fleitmann et al., 2009], in northwestern Turkey, identifies a wet period in the 1100s followed by a drought period in the 1200s, consistent with the variability identified in the OWDA (Figure 5). Cook et al. [2015] provide several comparisons with previous records relevant for the Mediterranean region. The OWDA scPDSI is positively correlated with the historical precipitation reconstruction of Pauling et al. [2006] over the western Mediterranean (southern Spain, Morocco, Algeria, and Tunisia), Greece and the Balkans, and Anatolia [Cook et al., 2015, Supplemental Figure 12]. The OWDA also identifies eight major historically documented droughts [White, 2006] that occurred during the Ottoman Empire in the late 1500s and early 1600s [Cook et al., 2015, Supplemental Figure 17].

3.3 Spatial Synchrony To further investigate the tendency for zonally synchronous hydroclimate variability across the basin, and meridional antiphasing in the eastern basin, we averaged scPDSI from the OWDA over three regions: the western Mediterranean (WestMED; 32°N–42°N, 10°W–0°E), the eastern Mediterranean (EastMED; 36°N–41°N, 20°E–37°E), and an area encompassing coastal Egypt, the southern Levant, and other areas of the Middle East (MidEast; 30°N–34°N, 33°E–47°E). The first two areas (WestMED and EastMED) correspond to approximately the same regions used in the analysis of Roberts et al. [2012]. Our MidEast box corresponds generally to the area of out‐of‐phase anomalies in Figures 5 and 6. Point‐by‐point correlations (Figure 7) between these three indices and the OWDA are strongly positive at the local level, as expected, decaying in magnitude outside of the averaging regions. For WestMED the correlations remain mostly positive, or near zero, across the entire basin. EastMED, however, shows strong negative correlations over Libya, Egypt, and the Levant, while MidEast is negatively correlated with a large region surrounding the Black Sea. As with the drought composites (previous section), these results again suggest a strong tendency for meridional antiphasing in hydroclimate in the eastern Mediterranean Basin, with a pattern similar to what would be expected due to precipitation responses to variations in the NAO. Results are similar when the correlations are calculated on single‐century subsets (1101–1200 C.E., 1201–1300 C.E., etc.) of the full data range (not shown), suggesting that the correlations across regions have been relatively stationary over time. Figure 7 Open in figure viewer PowerPoint Point‐by‐point Spearman's rank correlations (1100–2012) between OWDA scPDSI and the western Mediterranean (WestMED; 32°N–42°N, 10°W–0°), eastern Mediterranean (EastMED; 36°N–41°N, 20°E–37°E), and Middle East (MidEast; 30°N–34°N, 33°E–47°E) regional average scPDSI time series. Regions of insignificant correlation (p > 0.05) are masked by grey asterisks. We further investigate the nature and strength of drought variability across these regions through various spectral coherency analyses. The spectra of the WestMED and EastMED time series show significant power at interannual and decadal to multidecadal frequencies (Figure 8). Both regions have significant peaks at about 3–4 years, with EastMED additionally peaking at decadal frequencies and WestMED significant across a broader range of multidecadal bands. The two regions also overlap in their power at around 70 years, though EastMED is only marginally significant at the 90th percentile. A coherency spectra analysis between the two regional indices demonstrates highly significant coherence between the two regions on interannual and decadal timescales (Figure 9). Notably, there is also a broad band of coherence between the two regions at multidecadal and centennial timescales (30 to 130 years), despite no significant spectral peaks at wavelengths longer than ∼75 years in either WestMED or EastMED. This may be due to overlaps in shared frequencies: both series have peaks at around 70–75 years, approximately the first harmonic of the 130 year coherence in the coherency spectra. An analysis of the relative phasing for these spectral bands (not shown) indicates that the null hypothesis (i.e., WestMED and EastMED are in phase) cannot be rejected at the p≤0.05 significance level. Figure 8 Open in figure viewer PowerPoint Power spectral density (multitaper method, three tapers) for the WestMED and EastMED regional average scPDSI series. Red and black dashed lines are the 95th and 90th percentile confidence limits, respectively, estimated from 10,000 AR(1) series generated from the original time series. Figure 9 Open in figure viewer PowerPoint (top) Smoothed versions (10 year loess smooth) of the WestMED and EastMED time series and (bottom) the coherency spectra between unsmoothed versions of the two time series. Numbers in Figure 9 (bottom) highlight regions of significant coherency. These results are further confirmed through a cross‐wavelet coherency analysis (Figure 10), which shows that EastMED and WestMED share significant (p≤0.05) in‐phase variance (i.e., black arrows point to the right) in decadal to centennial frequency bands across the record. Further, EastMED and MidEast have significant coherence at interannual to multidecadal frequencies but are largely 180° out of phase (i.e., blacks arrows are primarily pointing to 90° to the left). These results indicate that over most of the last millennium there is a reasonably strong tendency for in‐phase drought in the zonal direction across the Mediterranean and out‐of‐phase variability in the eastern Mediterranean. Interestingly, given the association between the North Atlantic and western Mediterranean drought variability (Figure 3), there is no suggestion from the OWDA (Figure 9) that the NAO was in a persistent positive phase during the Medieval era. This is in contrast with conclusions from Esper et al. [2007] and Trouet et al.[2009] but is consistent with other regional climate reconstructions [Touchan et al., 2008a, 2011] and in agreement with a recent multiproxy reconstruction of the NAO by Ortega et al.[2015]. Figure 10 Open in figure viewer PowerPoint Squared wavelet coherence for (top) EastMED versus WestMED and (bottom) EastMED versus MidEast. Black arrows indicate where the two time series have significant (p≤0.05) shared variance. Relative phasing is indicated by the direction of the arrows; right pointing arrows indicate that series are in phase, and left pointing are 180° out of phase.