2.1 Dynamical systems metrics

In this study, we use a dynamical systems approach to compute two metrics: θ⁻¹ and α. The metric θ⁻¹, which we term persistence, is very intuitively a measure of the average residence time of the system around a state of interest. Hence, the higher the value of θ⁻¹, the more likely it is that the preceding and future states of the system will resemble the current state over relatively long timescales (Faranda et al., 2017b; Messori et al., 2017; Hochman et al., 2019). The metric α, which we term co-recurrence ratio, is a measure of the dynamical coupling between two variables, independently of their values (e.g. wet or dry) or in other terms their dependence structure.

The calculation of the dynamical systems metrics stems from the combination of Poincaré recurrences with extreme value theory (Lucarini et al., 2012; Freitas et al., 2010; Faranda et al., 2020). By recurrences we refer to the system being analysed returning arbitrarily close to a previously visited state in the phase space. Given an atmospheric variable x, we consider a state of interest ζ_x. In our case, this would be an instantaneous configuration of that variable, such as a latitude–longitude temperature map on a given day over the MED. We then consider recurrences to be those states that are close to ζ_x, namely other time steps at which the selected variable takes a very similar configuration. In order to quantify how close two configurations are to one another, we use the Euclidean distance (dist) between latitude–longitude maps. To compute the recurrences we first define an observable via logarithmic returns as follows:

where x(t) represents the time series of x. We then define a threshold s(q, ζ_x) as a function of high qth quantile of the time series g(x(t), ζ_x). Next, ∀g(x(t), ζ_x)>s(q, ζ_x) we define an exceedance u(ζ_x)=g(x(t), ζ_x)−s(q, ζ_x). The cumulative probability distribution F(u, ζ_x) then converges to the exponential member of the generalized Pareto distribution (Freitas et al., 2010; Lucarini et al., 2012):

where ϑ is the extremal index (Moloney et al., 2019), and we estimate it here following Süveges (2007). The dynamical systems persistence is computed as ${\mathit{\theta }}^{-\mathrm{1}}\left({\mathit{\zeta }}_x\right)=\mathrm{\Delta }t/\mathit{\vartheta }\left({\mathit{\zeta }}_x\right)$. In our case, Δt=1 d and ${\mathit{\theta }}_x^{-\mathrm{1}}$ has the units of the time step of the data being analysed (i.e. days). For conciseness, we hereafter adopt the notation ${\mathit{\theta }}_x^{-\mathrm{1}}$ to refer to the persistence of state ζ_x.

To extend the analysis to two variables, x and y, we can compute joint logarithmic returns around a state of interest ζ={ζ_x, ζ_y} as follows:

where $\Vert .\Vert$ represents the average root mean square norm of a vector's coordinates. Once joint logarithmic returns are defined, we compute the co-persistence ${\mathit{\theta }}_{x,y}^{-\mathrm{1}}$ based on the recurrences around ζ. This effectively amounts to a weighted average of ${\mathit{\theta }}_x^{-\mathrm{1}}$ and ${\mathit{\theta }}_y^{-\mathrm{1}}$ (Faranda et al., 2020; Abadi et al., 2018). In our analysis, the joint state ζ={ζ_x, ζ_y} would correspond to two instantaneous latitude–longitude maps: one for precipitation and one for temperature.

We further define the co-recurrence ratio (Faranda et al., 2020) α between x and y as

where s_x(q) and s_y(q) are high qth quantiles (or thresholds) of the univariate logarithmic returns g(x(t)) and g(y(t)) and ν[−] represents the number of events that satisfy condition [–]. Given a state ζ={ζ_x, ζ_y}, the co-recurrence ratio $\mathrm{0}\le \mathit{\alpha }\le \mathrm{1}$ measures the number of events where x resembles ζ_x given that y resembles ζ_y versus the number of cases when only x resembles the relevant reference state. When α=0, there are no co-recurrences of ζ={ζ_x, ζ_y} when we observe a recurrence of ζ_x. When α=1, recurrences of ζ_x are always also co-recurrences of ζ={ζ_x, ζ_y}. Hence, α may be interpreted as a measure of the dynamical coupling between x and y. However, α does not indicate causality: indeed, the order of x and y may be exchanged without affecting the value of α.

In order to compute the dynamical metrics we use a quantile q=0.98 to determine s. In previous studies (e.g. Faranda et al., 2011; Lucarini et al., 2012; Faranda et al., 2017b, 2019), this value has provided good estimates of the dynamical indicators, as it is high enough to select only genuine recurrences of ζ, while also ensuring a sufficiently large sample of recurrences for analysis. Tests further showed little sensitivity of the results to q in the range $\mathrm{0.95}<q<\mathrm{0.99}$ (Faranda et al., 2017b).

Finally, the dynamical systems approach rests on a number of theoretical assumptions, not all of which are strictly fulfilled by climate data. Specifically, the framework assumes the existence of an underlying chaotic attractor for the dynamics and was derived for ergodic systems (Freitas et al., 2010). However, recent applications have shown that weak nonstationarities do not preclude the validity of the results (e.g. Faranda et al., 2019), provided that they do not lead to bifurcations of the system. Unlike common statistical techniques (e.g. copulas), which rely on extrapolation of extreme values from statistical distributions, the metrics we use here are grounded in the underlying dynamics of the system being analysed.

In our analysis, we consider each daily time step in our datasets in turn as the state of interest ζ. The final result of our analysis is therefore a value for each metric and time step for the MED domain. This allows us to relate specific values of the metrics to the corresponding geographical anomaly patterns. We use the term compound dynamical extremes (CDEs) for the days characterized by α>90th quantile of the full-year distribution over the 1979–2018 period. We selected the 90th quantile as a good balance between an extreme value threshold and obtaining a sufficiently large sample of events. As a sensitivity test we repeated the analysis in Sect. 4.2 for a 95th quantile threshold, obtaining similar results (not shown). The two dynamical metrics successfully reflect large-scale features of atmospheric motions and have recently been applied to a range of different climate variables over different geographical domains (Faranda et al., 2017a, b, 2019, 2020; Messori et al., 2017; Rodrigues et al., 2018; Hochman et al., 2019, 2020; De Luca et al., 2020a; Scher and Messori, 2018).

2.2 Data

We use the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis over 1979–2018, with a spatial horizontal resolution of 0.25^∘ and a 6-hourly temporal resolution (C3S, 2017). Our MED domain follows the “Full Mediterranean” region described in Giorgi and Lionello (2008). For ERA5, this corresponds to 27.75–48.00^∘ N, 9.75^∘ W–39.00^∘ E. To improve the robustness of our results, we have repeated the bulk of the analysis on ERA-Interim (Dee et al., 2011) and the ERA5 10-member ensemble (C3S, 2017) (see Supplement). We use the instantaneous 6-hourly data to compute daily maximum and minimum 2 m temperature (K) and forecasted 1-hourly data for daily total precipitation (mm), from now on termed T_max, T_min and P respectively. Warm–dry days are days displaying positive T_max and negative P anomalies relative to JJA means. Similarly, cold–wet days are DJF days displaying negative T_min and positive P anomalies relative to DJF means. These are collectively referred to as “compound events” and the corresponding anomaly means are computed individually at a grid-point level. Therefore, if for example a grid point on a given day is warm, it does not necessarily imply that it is also dry. We also analyse daily-mean SLP (hPa) anomalies relative to JJA (DJF) means, computed from instantaneous 6-hourly steps.

2.3 Statistical tests

The statistical significance of the Sen's slopes (Sen, 1968) of the α and θ⁻¹ time series is verified using the Mann–Kendall test (Mann, 1945) from the R package https://cran.r-project.org/web/packages/modifiedmk/modifiedmk.pdf (last access: 28 August 2020). The Sen's slopes provide information about the steepness of the trends. If the Sen's slope is positive (negative), the corresponding trend is increasing (decreasing).

The statistical significance of SLP, T_max, T_min and P composite anomalies occurring during CDEs is computed using a one-tailed Mann–Whitney test at the 5 % confidence level (Mann and Whitney, 1947). The null hypothesis is that a randomly selected median anomaly value during a CDE is equally likely to be less than or greater than a randomly selected median value from the days that are not CDEs. The alternative hypothesis is that during JJA (DJF), the SLP and T_max (T_min) median anomalies observed during CDEs are higher (lower) than those observed during other days. For P in JJA (DJF), the alternative hypothesis is that anomalies observed during CDEs are lower (higher) than those during other days. To avoid incurring Type I errors (or false positives), we apply the Bonferroni correction to all p values when considering single-grid-point data (Bonferroni, 1936). The one-tailed Mann–Whitney test is also applied to the cumulative distribution functions (CDFs) of the anomaly means occurring during CDEs versus all other days.

Lastly, we checked the statistical significance of the percentage (%) agreement between JJA (DJF) CDEs and compound events. Here, the null hypothesis is that the JJA (DJF) observed percentage agreement is due to chance, and to compute the significance the following steps have been followed: (i) create n=1000 datasets of random dates, with the same number of elements in each dataset as we have for the CDEs; (ii) compute the percentage of agreement between CDEs and compound events' days for each dataset and grid point; (iii) pool together all the random percentage values and compute their 1st and 99th quantiles for each grid point; (iv) check whether the observed percentage values fall outside the random quantile values, and if this is the case, consider the percentage values statistically significant at the 1 % level (p value < 0.01).

4.1 Seasonality of CDEs

We next investigate the temporal distribution of CDEs. For α computed on T_max and P, all three reanalysis products display most of the CDEs clustering in July and August, with a secondary peak in DJF (Figs. 3a and S7a). For α computed on T_min and P, most CDEs occur during DJF, July and August (Figs. 3b and S7b). This holds for all three reanalysis products. The large number of CDEs during July and August (Figs. 3b and S7b) reflect compound warm–dry events (not shown). We further note that, notwithstanding the previously mentioned correlation between co-persistence and α, the seasonality of ${\mathit{\theta }}_{T_{\max },\mathrm{P}}^{-\mathrm{1}}$ extremes – defined analogously to the CDEs – does not reflect that of the CDEs (not shown). For both variable combinations (i.e. T_max–P and T_min–P), the two shoulder seasons (i.e. spring and autumn) display very few CDEs. In Faranda et al. (2017a), the authors hypothesized that during autumn and spring the atmospheric flow sits on a saddle-like point of the dynamics, while winter and summer represent more stable basins of attraction. Assuming that distinct attractors does indeed exist for winter and summer, we thus interpret these low CDE counts as the result of the atmospheric flow exploring both summer and winter configurations, resulting in rarer co-recurrences.

4.2 Pressure, temperature and precipitation anomalies during CDEs

During JJA, CDEs correspond to statistically significant positive SLP anomalies over the Western MED (north-western Africa) and the Anatolia–Black Sea region. These are separated by negative SLP anomalies spanning the Aegean Sea, the Levant and northern Egypt (Figs. 4a and S8a and b). These SLP anomalies are in turn associated with significant warm T_max anomalies over most of the MED, with a particularly warm Balkan Peninsula, and a negative anomaly over central northern Africa, to the east of the positive SLP anomaly (Figs. 4c and S8c and d). Lastly, we observe weak dry P anomalies over the Black Sea (Figs. 4e and S8e and f) and stronger wet P anomalies over the Alps. The latter correspond to statistically significant convective available potential energy (CAPE, J kg⁻¹) positive anomalies (Fig. S9) and may therefore be linked to localized convective P events. We conclude that JJA CDEs are closely linked to widespread warm T_max anomalies but have a weaker footprint on P anomalies, except over the Alps.

Figure 5Histograms and cumulative distribution functions (CDFs) of anomaly means of (a) T_max and (b) P during JJA CDEs and (c) T_min and (d) P during DJF CDEs. The data are the same as in Fig. 4c–f. The distributions are statistically different from those of all other JJA and DJF days respectively (p value ≪ 0.01, Mann–Whitney one-tailed test).

Download

In DJF we observe an east–west dipole in SLP over the MED, that favours cold-air advection from northern Europe to the Balkans, parts of the Italian Peninsula and the southern and Eastern MED (Figs.4b and S10a and b). Indeed, negative and significant T_min anomalies are observed over most of the MED region (Figs. 4d and S10c and d). The Eastern MED also displays significant positive P anomalies (Figs. 4f and S10e and f). The statistically significant (p value < 0.05) negative SLP anomalies over the Eastern MED are reminiscent of the footprint of Cyprus Lows, which are the main rain-bearing systems over the region (Alpert et al., 2004; Saaroni et al., 2010) (Figs. 4b and S10a and b). Cyprus Lows are also associated with the majority of wintertime cold spells over the Eastern MED (Hochman et al., 2020), and we do indeed find that some of the P anomalies over the Eastern MED are snowfall events, particularly over the Balkans, Turkey and Lebanon (Fig. S11). We thus conclude that CDEs are associated with wintertime cold–wet compound events over the Eastern MED.

As a proxy for the variability within our composites in Fig. 4, we compute the standard deviations (SDs) of the anomalies (not shown). We observe that SLP SDs are larger over the northern and central MED, while temperature SDs are larger over land compared to the sea – the latter a natural consequence of the sea's large thermal inertia. Finally, precipitation SDs are larger where the higher anomaly mean values are reported (i.e. the Alps in JJA and south-eastern MED in DJF), which may be linked to the prevailingly dry summertime conditions in the MED which yield low SDs where little or no rain falls. Similar results are obtained when computing Fig. 4 using only extreme anomalies (anom > 90th and anom < 10th quantiles) matching CDEs, although the JJA positive SLP anomalies are less geographically extensive (Fig. S12).

4.3 Distributions of temperature and precipitation anomaly means

We next test empirically whether the CDEs highlighted above have a systematic link to compound JJA warm–dry and DJF cold–wet events. During JJA, T_max and P daily anomaly means, computed for each grid point during CDEs, are predominantly warm (85 %) and dry (79 %) respectively (Fig.5a and b). Similar results are also obtained for ERA-Interim and the ERA5 ensemble (Fig. S13). P anomalies tend to cluster around zero, owing to the overall dry summertime climate of the region, although as noted above they do show a preference for negative (dry) values (Figs. 5b and S13b and d). A Mann–Whitney one-tailed test between the anomaly means during CDEs versus all other days in JJA results in statistically significant differences (p value ≪ 0.01) for all reanalysis products for both T_max and P. This implies that CDEs are significantly warmer and drier than other JJA days.

In DJF, most of the T_min and P anomaly means are cold (78 %) and wet (58 %) respectively for ERA5 (Fig. 5c and d) and the other reanalysis products (Fig. S14). Again, a Mann–Whitney one-tailed test between anomaly means during CDEs and all other DJF days highlights statistically significant (p value ≪ 0.01) differences for all reanalysis products, except ERA-Interim's P (p value < 0.05). This implies that CDEs are significantly colder and wetter than all other DJF days. The CDEs therefore present a somewhat mirror image of the preferred anomalies seen in the geographical anomaly composites for both JJA and DJF.

4.4 Spatial patterns of compound warm–dry and cold–wet events

We next complement the statistical information provided by the histograms and CDFs with spatial distributions of percentage (%) match between CDEs and compound events. Put simply, for each grid point in Fig. 6 we identify the days reporting compound events and CDEs and then divide the total number of these days by the total number of CDEs and multiply the resulting number by 100 to obtain the percentage agreement value. Across the MED, a high fraction of CDEs coincide with compound warm–dry events during JJA. Values locally exceed 70 %, meaning that >70 % of all JJA CDEs occur during compound warm–dry events (Figs. 6a and S15). The highest percentages occur in southern Spain, the Balearic Islands, Italy and the Balkans. During DJF, the percentage match between CDEs and compound cold–wet events is lower than that seen for warm–dry JJA events (<50 %) (Figs. 6b and S16). The highest percentage occurs over the Eastern MED sea, between the coastlines of Libya, Egypt, Greece and Turkey. In both JJA and DJF the vast majority of observations (%) are statistically significant at the 1 % level (p value < 0.01; Figs. 6 and S15 and S16).

Figure 6Percentage (%) of CDEs occurring during compound (a) JJA warm–dry and (b) DJF cold–wet events. The data are from the ERA5 reanalysis during 1979–2018. Stippling represent values not statistically significant at the 1 % level (p value ≥ 0.01).