Dynamical systems theory sheds new light on compound climate extremes in Europe and Eastern North America

We propose a novel approach to the study of compound extremes, grounded in dynamical systems theory. Specifically, we present the co‐recurrence ratio (α), which elucidates the dependence structure between variables by quantifying their joint recurrences. This approach is applied to daily climate extremes, derived from the ERA‐Interim reanalysis over the 1979–2018 period. The analysis focuses on concurrent (i.e., same‐day) wet (total precipitation) and windy (10 m wind gusts) extremes in Europe and concurrent cold (2 m temperature) extremes in Eastern North America and wet extremes in Europe. Results for wet and windy extremes in Europe, which we use as a test‐bed for our methodology, show that α peaks during boreal winter. High α values correspond to wet and windy extremes in northwestern Europe, and to large‐scale conditions resembling the positive phase of the North Atlantic Oscillation (NAO). This confirms earlier findings which link the positive NAO to a heightened frequency of extratropical cyclones impacting northwestern Europe. For the Eastern North America–Europe case, α extremes once again reflect concurrent climate extremes – in this case cold extremes over North America and wet extremes over Europe. Our analysis provides detailed spatial information on regional hotspots for these compound extreme occurrences, and encapsulates information on their spatial footprint which is typically not included in a conventional co‐occurrence analysis. We conclude that α successfully characterises compound extremes by reflecting the evolution of the associated meteorological maps. This approach is entirely general, and may be applied to different types of compound extremes and geographical regions.


INTRODUCTION
Individual natural hazards can interact (Gill and Malamud, 2014), often resulting in high-impact events leading to heavy socio-economic losses Zscheischler and Seneviratne, 2017;Zscheischler et al., 2018). These events are known as compound extremes, or multi-hazards, and they are explicitly defined by the United Nations Office for Disaster Risk Reduction (UNDRR) as (a) the selection of multiple major hazards that the country faces, and (b) the specific contexts where hazardous events may occur simultaneously, cascadingly or cumulatively over time, and taking into account the potential interrelated effects (UNDRR, 2017). In recent years, the concept of compound extremes, also extended to include multi-risks, has attracted the attention of the scientific community (e.g., Kappes et al., 2012;Liu et al., 2016;Vahedifard et al., 2016;Gill and Malamud, 2016;2017;AghaKouchak et al., 2018;Collet et al., 2018;Terzi et al., 2019). These extremes have potentially major implications for a wide range of private and public stakeholders including policy makers, (re)insurance companies, governments and local communities around the world (e.g., Fuchs et al., 2015;Franzke, 2017). Identifying compound extremes and quantifying their observed (De Luca et al., 2017;Ward et al., 2018;Couasnon et al., 2019) and future projected (Zscheischler and Seneviratne, 2017;Ben-Ari et al., 2018) spatio-temporal characteristics is thus a highly scientifically and socio-economically relevant goal, which supports disaster risk reduction through increased understanding of risks and enhancement of resilience (UNDRR, 2017;AghaKouchak et al., 2018;Zscheischler et al., 2018).
Statistical models are required to adequately assess the risk of compound extremes and to complement numerical models. However, compound extremes are inherently complex: they are rare, multivariate, and live in a sparsely sampled region of a multidimensional space with limited observational coverage. Statistical models thus rely heavily on extrapolation. A wide range of techniques have been proposed to study compound extremes, ranging from multivariate extreme value statistical models based on copula assumptions (e.g., Shiau, 2006;Bevacqua et al., 2017;Lee and Joe, 2018;Brunner et al., 2019), to max-stable models (e.g., Wang et al., 2014;Oesting and Stein, 2018), conditional exceedance models (e.g., Keef et al., 2009;2013;Neal et al., 2013;Zheng et al., 2014;Speight et al., 2017), Bayesian models (Kuczera, 1999;Madadgar and Moradkhani, 2014;Yan and Moradkhani, 2014;Kwon et al., 2016) and the multivariate skew-t distribution (Ghizzoni et al., 2010;. All these approaches can provide only partial information about the chaotic and spatial behaviour of the atmosphere. Indeed, incorporating the latter requires complex extensions of these methods (Genton et al., 2015;Hao et al., 2018), often with a risk of overfitting the statistical model.
Here, we propose a novel approach grounded in dynamical systems theory which overcomes some of these limitations, and aims to provide a complementary view to more traditional analyses issuing from statistics and climate dynamics (e.g., Martius et al., 2016;Waliser and Guan, 2017). We specifically propose an objective measure of the co-recurrence of extremes in different climate variables, which can provide both temporal and spatial information and can be linked to the underlying dynamical properties of the climate system, as well as explicitly accounting for its underlying chaotic nature. This is complemented by two other metrics characterising the evolution of large-scale climate fields (Lucarini et al., 2016;Messori et al., 2017;Faranda et al., 2017b). The approach is very flexible, and can in principle be applied to any number of variables, geographical regions and types of dataset. Within the text, we adapt the vocabulary used in the literature to fit the novel approach we propose. We will be referring to compound extremes when considering bivariate metrics calculated from dynamical systems theory and to concurrent extremes when discussing the joint (i.e., same-day) occurrencesof extreme values in climate variables.
We begin by providing a description of our data (Section 2) and methodology (Sections 3 and 4). Next, we illustrate the physical interpretation of the different dynamical systems metrics and their seasonality (Section 5.1.1), and present an application to concurrent wet and windy extremes in Europe (Section 5.1.2), along with concurrent cold extremes in Eastern North America and wet extremes in Europe (Section 5.2). We thus address concurrent extremes in both single and separate geographical regions. Our choice is motivated by the significant scientific and media attention that both wet and windy weather in Europe (e.g., Pinto et al., 2013;Huntingford et al., 2014;van Oldenborgh et al., 2015;Wild et al., 2015) and cold spells in North America (e.g., Lee et al., 2015;Trenary et al., 2015;Harnik et al., 2016;Messori et al., 2016) have attracted in recent years. The analysis in Section 5.1 primarily aims to relate the information provided by our metrics to known features of the atmospheric variability over the Euro-Atlantic region, and can be viewed as a test of the robustness of the metrics we propose. Section 5.2 instead illustrates how our approach can be used to gain novel insights into concurrent extremes in two different regions. We conclude by discussing our results in the context of the literature on these classes of extremes (Section 6).

DATA
We use the ERA-Interim reanalysis (Dee et al., 2011) from the European Centre for Medium-range Weather Forecasts (ECMWF) , at a horizontal resolution of 0.75 • , from 1979 to 2018. ERA-Interim outperforms other reanalysis products, such as MERRA, NCEP-NCAR, ERA-40, CFSR and GLDAS 1 , in reproducing temperature and wind speed observations. Precipitation is notoriously problematic to capture on a relatively coarse horizontal grid, but ERA-Interim nonetheless performs comparatively well in this respect (Decker et al., 2012). The analysis was conducted on daily data at 1200 and 0000 UTC (12 hr step) of total precipitation (mm) and 10 m wind gusts (m⋅ s −1 ). To obtain a unique daily value, these were then summed and averaged, respectively. We also used 2 m temperature (K) at 0000, 0600, 1200 and 1800 UTC, which were averaged as daily observations. From now on we will refer to these variables simply as precipitation, wind and temperature.  Barnston and Livezey, 1987), the dominant mode of climate variability within the North Atlantic region, were downloaded from the NOAA Climate Prediction Center website (https://www.cpc.ncep.noaa. gov/products/precip/CWlink/pna/nao.shtml; accessed 5 February 2020). The analyses were also repeated by making use of daily maximum wind gust observations (m⋅ s −1 ), that is, the daily maximum of the 1200 and 0000 UTC values for the 1979-2018 period. The results are qualitatively very similar to daily wind means (see Supporting Information).

Qualitative interpretation of the dynamical systems metrics
We use three dynamical systems metrics: the local dimension (d), the local persistence ( −1 ), and the local co-recurrence ratio ( ). d and −1 are calculated for both univariate and bivariate cases, whereas always refers to two different variables. All three metrics are local in phase-space, and instantaneous in time. They can thus be computed for any state on the underlying attractor -a 1 Modern-Era Retrospective Analysis for Research and Applications; National Centers for Atmospheric Prediction -National Center for Atmospheric Research; ECMWF re-analysis from September 1957 to August 2002; NCEP's Climate Forecast System Reanalysis; Global Land Data Assimilation System state in our case being a latitude-longitude map of one or more atmospheric variables at a given time. d and −1 have recently been applied to a range of different climate variables over different geographical domains and were found to successfully reflect large-scale features of atmospheric motions (Messori et al., 2017;Faranda et al., 2017a;2017b;Rodrigues et al., 2018;Hochman et al., 2019a;. The local dimension d( ) describes the evolution of the system around , and can be interpreted as a proxy for the number of degrees of freedom active around the state of interest. The local persistence −1 ( ) measures the mean residence time of the system around such state. A high (low) persistence indicates that the system's evolution will slowly (rapidly) lead to a dynamically different state. The higher the value of −1 , the more likely it is that the preceding and future states will resemble over comparatively long time-scales (Faranda et al., 2017b;Messori et al., 2017;Hochman et al., 2019a). The two dynamical systems metrics can be computed for both individual variables (e.g., sea-level pressure, Faranda et al., 2017b) and multiple variables jointly (Faranda et al., 2020). Given two different variables, we can define a state on the Poincaré section jointly spanned by x and y as = { x , y }. We will use the notation d x or −1 x to refer to the monovariate metrics computed for the atmospheric variable x and the notation d x,y or −1 x,y to refer to the bivariate metrics jointly computed on the atmospheric variables x and y.
Along with d and −1 , we adopt a phase-space local measure of dependence between variables first introduced in Faranda et al. (2020). We term this the co-recurrence ratio ( ). For a given state of interest = { x , y }, we have 0 ≤ ( ) ≤ 1. When ( ) = 0, no recurrences of = { x , y } are observed in the phase-space when we observe a recurrence of x . Intuitively, when ( ) = 1, all recurrences of x correspond to recurrences of = { x , y }, and vice versa. In other words, quantifies the co-recurrences of two (or more) variables within a hyper-ball around .

Derivation of the dynamical systems metrics
We view the atmosphere as a dynamical system with an observed phase-space trajectory x(t), and consider the properties of such a system near a state of interest . To compute our dynamical systems metrics, we first define our observable, via logarithmic returns as: where dist is the Euclidean Norm. This observable discriminates parts of the trajectory x(t) that are close to (i.e., when the distance approaches 0) from those that are further away. We then define s(q, ) a high qth quantile of the time series g{x(t), }, and we define ∀ g{x(t), } > s(q, ) an exceedance u( ) = g{x(t), } − s(q, ). According to the Freitas-Freitas-Todd theorem (Freitas et al., 2010), later modified in Lucarini et al. (2012), the cumulative probability distribution F(u, ) then converges to the exponential member of the Generalised Pareto Distribution: The parameters and depend on the chosen point . is the so-called extremal index (Moloney et al., 2019), and we estimate it here using the approach of Süveges (2007). The local dimension d is obtained as d( ) = 1∕ ( ), while the local persistence is given by −1 =Δt∕ , where Δt is the timestep of the data used and the local persistence is in the same units as the timestep. In our case, trivially Δt = 1 day.
The bivariate d and −1 metrics are derived following an analogous procedure (Faranda et al., 2020). Given two trajectories x(t) and y(t) -in our analysis two different atmospheric variables -and a state = { x , y }, we can define d x ( ) and d y ( ). These are the dimensions of the Poincaré sections defined by x and y around , with respect to the chosen distance metric. We can further define joint logarithmic returns of as: where ||.|| represents the average root mean square norm of a vector's coordinates. One can then obtain the co-dimension d x,y analogously to the derivation of the univariate d (for ease of notation we hereafter drop the dependence on ). The co-dimension has the following property: The case d x,y = d x + d y implies no coupling bewteen x and y. If instead the variables are deterministically coupled (i.e., x is a function of y or vice versa), d x,y = min(d x , d y ).
Following the definition of joint logarithmic returns given in Equation (3), one can also define the local co-persistence −1 x,y . This is defined as a weighted average of −1 x and −1 y , where the weights depend on the size of the hyper-ball around in the Poincaré sections x and y (Abadi et al., 2018;Faranda et al., 2020).
Along with d and −1 , we adopt a phase-space local measure of dependence termed co-recurrence ratio (Faranda et al., 2020). For two variables x(t) and y(t), the co-recurrence ratio 0 ≤ ( ) ≤ 1 of a state = { x , y } is defined as a ratio: 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 [−].

Definition of anomalies and extremes
We define anomalies relative to a daily climatology. For example, the climatological temperature at a given gridpoint for 19 February is the mean of all 19 February temperatures at that location over the 40 years considered here. Precipitation and wind extremes are defined as daily values >99th quantile, whereas low temperature extremes are daily observations <1st quantile of the respective distributions at each gridpoint. The selection was performed on absolute values, rather than anomalies, in the spirit of an impacts-based perspective. Indeed, a very strong windstorm is more likely to cause damage than an out-of-season (and thus extreme in terms of anomalies) but moderate (in terms of absolute wind-speed values) windstorm, although out-of-season extremes are associated with their own specific sets of risks. The same analysis was also performed for extremes based on anomalies of the climate variables and the two approaches show good spatio-temporal agreement (see Supporting Information).
The anomalies of precipitation, wind and temperature are stratified according to daily extremes (>99th quantile), thus informing on the behaviour of the climate variables under strong dynamical coupling conditions. We also consider, for each gridpoint, the composite mean anomalies associated with univariate and concurrent extremes of precipitation, wind and temperature. A list of extreme occurrences at each gridpoint was compiled, and then the corresponding daily anomalies were composited, thus providing a different composite anomaly value at each gridpoint. As an example, consider precipitation and wind extremes for a gridpoint within the European region. First, we identified days when the two extremes (co)occurred at that location, and then we computed the mean anomalies based on those days (note that the daily values are the same for both types of extremes, since it is computed using both variables simultaneously). For the case of concurrent cold and wet extremes in Eastern North America and Europe respectively, we adopted a similar procedure.
We first identified dates on which at least one gridpoint in each domain displays an extreme observation. We next identified all gridpoints displaying extreme values on each selected date. Finally, we calculated the mean anomalies for each selected gridpoint over the two regions. Note that, again, daily values are the same for both variables, regardless of the fact that they are taken over different geographical domains.

Statistical testing
We evaluate the significance of the anomalies in climate variables associated with extremes using the Mann-Whitney test (Mann and Whitney, 1947). This consists of a non-parametric procedure that, given two samples x and y of size n x and n y , tests the null hypothesis that the samples are drawn from populations with the same median. We perform the test in each grid point, in such a way that y is the sub-sample selected using extreme values of , and x is the remainder of the time series in that grid point. We choose the Mann-Whitney test since the data considered in this analysis, especially wind and precipitation, are strongly non-Gaussian and thus do not meet one of the assumptions of the Student's t-test on the sample means. Nevertheless, an exact application of the test would require independently and identically distributed samples, while time series of atmospheric variables are expected to exhibit at least significant autocorrelation, or even nonlinear dependence in time. In general, the presence of serial dependence in the data reduces the power of the test, that is, the probability to not reject the null hypothesis when the median difference is significant is higher than in the case of samples meeting all the test assumptions (e.g., Yue and Wang, 2002). If there are realistic expectations about the dependence structure in the samples, the distribution of the test statistics for the dependent case can be approximated using Monte Carlo sampling. However, this is not the present case, and we decide to apply the regular Mann-Whitney test, aware of the possible loss of power. Indeed, we believe that this issue should not excessively affect the results, since the test is performed over a large number of grid points, and our expectation is to observe regions of significant effects. While the total number of significant points in one of such regions may be reduced due to loss of power, we expect it to be still possible to draw robust general conclusions.
The procedure consists of the following steps: a the two samples are merged and ranked; b the sum of the ranks is computed separately for data belonging to sample x and y, and denoted R x and R y respectively; c the Mann-Whitney test statistic U is obtained as U x = R x − n x (n x + 1)∕2.
In the case of ties, the rank in the midpoint between the two closest non-tied ranks is used. Under the null hypothesis, U follows a known tabulated distribution, and for samples large enough (n ≳ 20) converges to a normal distribution. Notice that U can be computed indifferently using sample x or y. We state the alternative hypothesis based on prior knowledge that dynamical extremes should correspond to large positive precipitation and wind anomalies, and to large negative temperature anomalies. Then, the alternative hypothesis is that the median of the sample selected with is larger than the median of the rest of the time series for precipitation and wind, and smaller for temperature.
A similar procedure was also applied to check whether the mean anomaly values associated with univariate extreme climate variables were significant. In particular, we perform again the Mann-Whitney test, but this time the two samples x and y consists of the percentiles of each value of in the distribution of the merged samples of mean anomalies. The alternative hypothesis is now that the values of corresponding to climate extremes are -in median -significantly larger than the ones not corresponding to climate extremes. The strategy for concurrent extremes is similar to the univariate case, merging the time series of different climate variables by dates (i.e., precipitation-wind in Europe and low temperatures-precipitation over Eastern North America and Europe), so that only extreme observations co-occurring on the same days are kept, along with the associated daily observations. We then repeat the procedure according to the univariate case above.
The test is conducted at the standard level a 0 = 0.05. However, the level requires a correction, due to the fact that, for each event, we perform as many tests as the number of gridpoints, n gr . This increases the probability of erroneously rejecting the null hypothesis, leading to false positives. On the other hand, in the case of concurrent extremes, a number n NA of tests is skipped due to missing data. To overcome this issue, we apply a Bonferroni correction (Bonferroni, 1936), which consists of using the original level divided (or corrected) by the effective number of tests, a = a 0 ∕(n gr − n NA ).
The robustness of the anomaly means observed during univariate and concurrent extremes was also assessed by performing a sign test. This quantifies the fraction of individual extreme observations at a given gridpoint which have the same sign anomaly as the overall composite. For example, a sign test value of 66% at a location with a positive composite anomaly means that 66% of the days used to create the composite display a positive daily anomaly, F I G U R E 1 (a) Daily climatological means of the precipitation-wind co-recurrence ratio ( ). The blue shaded areas represent one standard deviation from the mean. (b) Cumulative distribution functions (CDFs) of local dimension (d) anomalies associated with daily extremes (>99th quantile). (c) is as (b) but for local persistence ( −1 ). In (b, c) the metrics are calculated for precipitation (brown), wind (yellow) and jointly for precipitation-wind (blue). The data cover 1979-2018 over Europe while 34% display a negative anomaly. Clearly, the higher the %age of members sharing same-sign anomalies, the more robust the results are.

5.1
Wet and windy extremes in Europe

Dynamical climatology and compound extremes
In Figure 1 we show both univariate and bivariate dynamical metrics for precipitation and wind in Europe. peaks during boreal winter (December-February) with values in excess of ∼0.17 and reaches its lowest values during spring (∼0.075) and summer (∼0.10; Figure 1a).
The winter peak follows the intuition that during these months the passage of extratropical cyclones (ETCs) over western-continental Europe and the British Isles (BI) can lead to widespread flooding episodes and comparatively frequent concurrent wet and windy episodes, whose spatial footprints resemble one another (Fink et al., 2009;Lavers et al., 2011;Matthews et al., 2016;De Luca et al., 2017). On the other hand, the lower values observed during spring and summer may reflect precipitation events which are mainly driven by small-scale convective processes (e.g., Pieri et al., 2015) and therefore have a weaker link to both the synoptic-scale and larger-scale circulation (Bisselink and Dolman, 2008).
During extremes (i.e., strong dynamical coupling), negative d anomalies are observed for both the compound and univariate cases (Figure 1b), with the former and wind displaying the more numerous negative anomalies. An almost opposite pattern of positive values is found for −1 F I G U R E 2 (a) Daily NAOI and (b) monthly counts for daily extremes (>99th quantile). In (a) the black solid line represents NAOI = 0, the orange dashed line the mean NAOI value for extremes, and the violet dot-dashed line the NAOI mean 97.5th quantile obtained by random sampling (n = 10, 000; no replacement) from the entire NAOI daily time-series. The data cover 1979-2018 over Europe anomalies (Figure 1c), with precipitation again being the variable showing the weaker deviation from climatology. Low d and high −1 are associated with predictable configurations (Faranda et al., 2017b;Messori et al., 2017;Scher and Messori, 2018), and we find here that this matches a strong coupling between the large-scale configurations of the two variables we are analysing. We next relate to the NAO index (NAOI; Barnston and Livezey, 1987). The NAOI daily values during extremes (Figure 2a) are overwhelmingly positive (mean = 0.85, orange dashed line), and well above about what may be obtained by chance (violet dot-dashed line). When the NAO is in a positive phase, Northern/Western Europe experiences enhanced storminess from the passage of ETCs, which bring concurrent flooding and windstorms (Huntingford et al., 2014;De Luca et al., 2017;Nobre et al., 2017;. This agrees with the strong negative and positive anomalies respectively in d and −1 discussed above, since the positive NAO has been previously identified as a low-d, high-−1 configuration (Faranda et al., 2017b). Moreover, the vast majority of extremes occur during late autumn and winter, with no observations in May and summer (June-August; Figure 2b), further confirming the plausibility of the link between extremes and winter ETCs.
Very similar results are obtained when repeating the above analysis on daily wind maxima (Figures S1 and S2), instead of wind means, computed over two daily observations. For the maxima, daily climatological means are lower ( Figure S1a), likely reflecting the intrinsically noisier nature of the variable.

Concurrent extremes
Here, we consider the link between extremes in computed on the temperature and wind fields and climate extremes. During extremes, positive anomalies in both the precipitation and wind maps are observed in Northern/Western Europe, over a broad region spanning the BI, the North Sea and the surrounding coastal areas. Negative anomalies are instead observed in Southern Europe and the Mediterranean (Figure 3a,b). This pattern resembles that associated with the NAO (Figure S3), further strengthening our interpretation that extremes are able to capture the impacts of ETCs in Northern/Western Europe. The number of concurrent wet and windy extremes peaks during late boreal autumn and winter (Figure 3c). This closely reflects the timing of the extremes (Figure 2b), providing additional evidence for a close link between compound dynamical and concurrent climate extremes.
Finally, we test the inverse link, namely whether large values can be recovered when conditioning on extremes in the climate variables. The anomaly means, conditioned on monovariate precipitation or wind extremes, are highest in Northern/Western Europe, although for wind significant anomalies span virtually the whole continent F I G U R E 3 (a) Precipitation and (b) wind anomaly means associated with extremes (>99th quantile). Stippling represents statistically significant values. (c) shows monthly counts of concurrent (i.e., same-day) precipitation-wind extremes. anomaly mean values observed during (d) precipitation, (e) wind and (f) concurrent precipitation-wind extremes (>99th quantile). In (d)-(f) only statistically significant values are plotted and stippling represent gridpoints with >66% of daily positive anomalies (sign test). Statistical significance is assessed via the Mann-Whitney test (Figure 3d,e). The signal for precipitation is somewhat patchy, but there is an overall good spatial correspondence with the regions displaying large anomalies in the climate variables themselves (Figure 3a,b). We also note that the anomaly means related to wind extremes are higher than those for precipitation. This is not surprising, since precipitation has a noisy and fractal structure, which is reflected in the dynamical system metrics (Faranda et al., 2017a). anomaly means calculated for concurrent wet and windy extremes mostly reflect the patterns observed for the univariate extreme cases (Figure 3f). The largest values are again distributed over the BI, Northern/Western Europe and the Nordic Seas. The anomaly means for the concurrent extremes show even higher values than for wind alone. We interpret this as a direct result of the metric reflecting the dynamical coupling of the two variables. The sign test (stippling in Figure 3d-f) largely mirrors the magnitude and spatial pattern of the anomaly means. The same analyses as shown in Figure 3 and Figure 3c-f were respectively repeated by using daily wind maxima observations ( Figure S4) and a different definition of extremes (Section 4.1 and Figure S5). The results are consistent between the three approaches.

Cold and wet extremes in Eastern North America and Europe
In this section, we consider temperature over Eastern North America and precipitation over Europe (Figure 4). is computed on these two variables, each in its respective domain. extremes are associated with negative temperature anomalies across Eastern North America, with the largest anomalies situated in the northeastern part of the domain (Figure 4a). Precipitation anomalies are instead mostly positive, and are strongest over Iberia, Western France and Iran (Figure 4b). The vast majority of the concurrent cold and wet extremes are observed during the boreal winter months (Dec-Feb). No co-occurrences are observed from April to October (Figure 4c). Statistically significant anomaly means associated with cold extremes are observed over the northeastern part of the North American domain (Figure 4d), reflecting the larger temperature anomalies shown in Figure 4a.
anomaly means linked with wet extremes in Europe are instead significant over Scandinavia, northwestern Russia and in the Middle East (Figure 4e). The difference in geographical distribution relative to the composite precipitation anomalies shown in Figure 4b likely reflects differences in the local precipitation distributions. That is, the fact that the large anomalies shown in Figure 4b are locally less extreme than the weaker anomalies found over Scandinavia. Note that Figure 4e differs from Figure 3e because of the different variables and domains over which is computed here. Lastly, anomaly means observed during concurrent cold and wet extremes show statistical significance in the northeastern North American domain and Iberia/Western France (Figure 4f). These reflect very closely the mean anomaly patterns in Figure 4a,b. Similar to what we found for concurrent wet and windy extremes in Europe (Figure 3f), these anomaly means for the concurrent extremes display larger values than for the univariate extremes. The monovariate anomaly means mostly fail the sign test, while the anomalies for the concurrent extremes display extensive sign agreement mainly over Europe (stippling in Figure 4d-f). We interpret these differences as indicating that the computed on the North American temperature and European precipitation reflect more closely the concurrent extremes in these variables than the monovariate extremes in the individual domains. As for Figure 3c-f, similar results are retrieved for a different definition of extremes ( Figure S6).

DISCUSSION AND CONCLUSIONS
In this study we applied a new approach issued from dynamical systems theory (Lucarini et al., 2016;Faranda et al., 2017b;Messori et al., 2017;Faranda et al., 2020) to the study of compound and concurrent climate extremes. We use the terms compound extremes when considering bivariate metrics calculated from dynamical systems theory and concurrent extremes when discussing the joint occurrences (i.e., same-day) of extreme values in the climate variables. We specifically focused on wet and windy extremes in Europe, and cold and wet extremes in Eastern North America and Europe, respectively. We defined a coupling parameter (the co-recurrence ratio) that measures the joint Poincaré recurrences of specific configurations in two variables. Almost by definition, is sensitive to the chosen geographical domain(s). Indeed, the dynamical coupling between the variables of interest can change significantly when including or excluding data. It is therefore important to limit the analysis only to regions of interests, whatever the criterion for "interesting" may be. Two more dynamical systems metrics were also computed for both uni-and bivariate cases, namely the local dimension (d) and the local persistence ( −1 ). These inform on the number of degrees of freedom active around a given state of the system (here the atmosphere in the Euro-Atlantic and Eastern North American sectors, as represented by the temperature, wind gust and precipitation maps) and on the mean residence time of the system around such state, respectively.
To verify whether our approach can provide new insights into the nature of concurrent extremes, we compare our results for the wet and windy extremes in Europe against those obtained from a Canonical Correlation Analysis (CCA) (Härdle and Simar, 2007), which identifies the linear combinations of two variables with maximum correlation (Figures S7-S9). The spatio-temporal patterns issuing from the CCA largely resemble the ones found in Figures 1a, 2, and 3, lending credence to the robustness of our results. Specifically, the CCA correlation coefficient F I G U R E 4 As Figure 3, but for temperature in Eastern North America and precipitation in Europe. In (c), (d), and (f) temperature extremes are daily values <1st quantile (cc) peaks during winter, albeit with a more marked seasonal pattern ( Figure S7) than (Figure 1a). The cc extreme values (>99th quantile) match positive NAO phases, mainly occurring during winter months ( Figure  S8), in very close agreement with Figure 2. Lastly, similar results are also obtained when replicating Figure 3, with the cc as a measure of dependence ( Figure S9). The main difference is that in Figure S9a-b,e, the signal is weaker than in Figure 4a-b,e, suggesting that provides a more nuanced picture of the extremes. Additional reasons motivating the use of the co-recurrence metric we propose here over a CCA approach are that is linked to other dynamical properties of the phase space underlying the atmospheric motions, such as local dimension and persistence. It therefore indirectly provides a rich set of information on the system, including on its intrinsic predictability (Messori et al., 2017;Scher and Messori, 2018;Faranda et al., 2019b). One may, for example, further investigate the predictability implications of the fact that concurrent precipitation-wind extremes over Europe display a low d and high −1 , and compare this to the empirical results from ensemble reforecast data. Moreover, the definition of can be easily extended to a multivariate case with more than two variables, whereas the CCA framework requires more complex adaptations (such as partial CCA) when working with more than two sets of variables.
The same approach could in the future be applied to a suite of renalysis products, to evaluating large ensembles of historical Atmosphere-Ocean General Circulation Models runs (e.g., Taylor et al., 2011;Eyring et al., 2016), or to projecting future changes in compound extremes under a warming climate. Moreover, the dynamical systems approach used here can be implemented for any geographical region and variable (as long as the latter shows a chaotic behaviour). A diverse set of compound and concurrent events can therefore be tested, and global hotspots of compound extremes identified. There is also the hope that this methodology may help to investigate outstanding scientific questions; for example the role of Arctic Amplification in driving Northern Hemisphere midlatitude climate extremes (e.g., Barnes, 2013;Cohen et al., 2014;Mori et al., 2014;Screen and Simmonds, 2014;Coumou et al., 2018;Ye and Messori, 2019).
The main findings of the study can be summarised as follows: • The co-recurrence ratio captures the seasonal cycle in extreme event occurrences for the domains (Eastern North America, Europe) and variables (precipitation, 10 m wind gusts, 2 m temperature) considered here.
• extremes, here termed compound extremes, match same-day extremes in the climate variables, here termed concurrent extremes. High values indicate a strong coupling between the different climate variables, which intuitively should be the case for spatially and/or temporally co-occurring extremes.
• For the European domain, the close correspondence between extremes in and in the climate variables is mediated by the strong projection of extremes on the positive phase of the NAO.
• For the concurrent North American and European extremes, highlights regional extreme event hotspots. These are regions characterised by similar footprints of the extreme events across the different co-occurrence instances.
These findings highlight the ability of dynamical systems metrics to capture the spatio-temporal variability of the atmosphere. In Europe, our results mostly confirm the well-known link between a positive NAO and stormy weather over the northeastern part of the continent, leading to widespread flooding linked with severe winds (Fink et al., 2009;Lavers et al., 2011;Huntingford et al., 2014;Matthews et al., 2016;De Luca et al., 2017). The physical relationship between low temperatures in Eastern North America and wet extremes in Europe is not as obvious. The notably wet and stormy winter of 2013/2014 in parts of Europe and the British Isles (e.g., Huntingford et al., 2014;van Oldenborgh et al., 2015;Wild et al., 2015) co-occurred with a very cold winter in the Eastern United States (Lee et al., 2015;Trenary et al., 2015). However, the longer-term climatology of concurrent winter cold spells in North America and wet and stormy events over western Europe has received rather less attention. Messori et al. (2016) explicitly addressed the question of possible physical links between the two sets of extremes, and concluded that cold spells in North America are a systematic, but not necessary, precursor to stormy conditions in western Europe. The question of whether different instances of co-occurring extremes shared similar spatial footprints -of high relevance for investigating their predictability -was left largely unanswered. Our analysis highlights southwestern Québec, southeastern Ontario and part of the western European coastline from Portugal to France as hotspot regions for concurrent extremes with spatially similar characteristics over time. The European region roughly matches that highlighted in Messori et al. (2016) (their figure 4d). However, the North American region is shifted to the northeast of the region analysed in a number of studies dealing with the winter of 2013/2014 and, more generally, eastern North American cold spells (e.g., Trenary et al., 2015;Harnik et al., 2016;Messori et al., 2016;Overland and Wang, 2017). We reiterate that our analysis specifically points to regions characterised by concurrent extremes with similar spatial footprints across the different co-recurrence instances, which may explain the discrepancy with some of the results in the literature. Finally, we note that the interpretation of the dynamical systems analysis may pose some challenges -as is for example the case of the anomalies associated with precipitation considered in Section 5.2.
The spatial information on the extreme events intrinsic to our approach can thus provide new insights into multi-hazards (or compound events) research. Our dynamical systems perspective is flexible, links directly to concepts of predictability and coupling, and may be useful in contexts which are both practical (e.g., numerical forecasting of extremes) and more theoretical (e.g., understanding the atmospheric drivers of extremes).