{'fraction_non_alphanumeric': 0.039873961615582924, 'fraction_numerical': 0.017158407333142365, 'mean_word_length': 4.725930785632278, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical Systems Theory Sheds New Light on Compound Climate Extremes in Europe and Eastern North America

1Geography and Environment, Loughborough University, Loughborough, UK 2Department of Earth Sciences, Uppsala University, Uppsala, Sweden 3Centre of Natural Hazards andDisaster Science (CNDS), Uppsala, Sweden 4Department ofWater and Climate Risk, Vrije Universiteit Amsterdam, Amsterdam, The Netherlands 5Department ofMeteorology, Stockholm University and Bolin Centre for Climate Research, Stockholm, Sweden 6Laboratoire des Sciences du Climat et de l’Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France 7LondonMathematical Laboratory, London, UK

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 (10m wind gusts) extremes in Europe and concurrent cold (2m 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 north-western 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 extra-tropical cyclones impacting north-western Europe. For the Eastern North America-Europe case, α extremes once again reflect concurrent climate extremesin this case cold extremes over North America and wet extremes over Europe. Our analysis provides detailed spatial information on regional hotspots for these compound ex-1 treme 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 et al., 2018;Zscheischler and Seneviratne, 2017;Leonard et al., 2014). 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 (1) the selection of multiple major hazards that the country faces, and (2) the specific contexts where hazardous events may occur simultaneously, cascadingly or cumulatively over time, and taking into account the potential interrelated effects (UNDRR, 2017).
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. Bevacqua et al., 2017;Brunner et al., 2018;Lee and Joe, 2018;Shiau, 2006), to max-stable models (e.g. Oesting and Stein, 2018;Wang et al., 2014), conditional exceedance models (e.g. Keef et al., 2009Keef et al., , 2013Neal et al., 2013;Speight et al., 2017;Zheng et al., 2014), Bayesian models (Kuczera, 1999;Kwon et al., 2016;Madadgar and Moradkhani, 2014;Yan and Moradkhani, 2014) and the multivariate skew-t distribution (Ghizzoni et al., 2010(Ghizzoni et al., , 2012. 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 (Hao et al., 2018;Genton et al., 2015), 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;Faranda et al., 2017b;Messori et al., 2017). 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 occurrences (i.e. same-day) of extreme values in climate variables.
We begin by providing a description of our data (Section 2) and methodology (Sections 3, 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;van Oldenborgh et al., 2015;Huntingford et al., 2014;Wild et al., 2015) and cold spells in North America (e.g. Trenary et al., 2015;Lee et al., 2015;Messori et al., 2016;Harnik 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).
ERA-Interim outperforms other reanalysis products, such as MERRA, NCEP-NCAR, ERA-40, CFSR and GLDAS, 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 12 and 00 UTC (12h step) of total precipitation (mm) and 10m wind gusts (ms −1 ). To obtain a unique daily value, these were then summed and averaged, respectively. We also used 2m temperature (K ) at 00, 06, 12 and 18 UTC, which were averaged as daily observations. From now on we will refer to these variables simply as precipitation, wind and temperature. We consider two domains: Europe (18 • W-51 • E, Oscillation Index (NAOI) (Barnston and Livezey, 1987), the dominant mode of climate variability within the North Atlantic region, were downloaded from the NOAA -Climate Prediction Center website (available here). The analyses were also repeated by making use of daily maximum wind gust observations (ms −1 ), i.e. the daily maximum of the 12 and 00 UTC values for the 1979-2018 period. The results are qualitatively very similar compared 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 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 (Faranda et al., 2017a(Faranda et al., ,b, 2019a(Faranda et al., , 2020Messori et al., 2017;Rodrigues et al., 2018;Hochman et al., 2019a,b).
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 timescales (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 viceversa. 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 system near a state of interest ζ. To compute our dynamical systems metrics, we first define our observable, via logarithmic returns as: Where d i st 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 q -th 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 Generalized 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 codimension 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 viceversa), d x,y = min(d x , d y ).
Following the definition of joint logarithmic returns given in Eq. 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 (Faranda et al., 2020;Abadi et al., 2018).
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 q -th 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 the 19 t h February is the mean of all 19 t h February temperatures at that location over the 40 years considered here.
Precipitation and wind extremes are defined as daily values > 99 t h quantile, whereas cold temperature extremes are daily observations < 1 st 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 (> 99 t h 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 nonparametric 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 remaining 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, i.e. 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: i) the two samples are merged and ranked; ii) the sum of the ranks is computed separately for data belonging to sample x and y , and denoted R x and R y respectively; iii) the Mann-Whitney U test statistics is obtained as U x = R x − n x (n x + 1)/2. In 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 mediansignificantly 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 cold 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 as per 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 g r . 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 N A 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 g r − n N A ).
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, while 34% display a negative anomaly. Clearly, the higher the % of members sharing same-sign anomalies, the more robust the results are.

| 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 (Dec-Feb) 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 extra-tropical 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 (De Luca et al., 2017;Fink et al., 2009;Lavers et al., 2011;Matthews et al., 2016). 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 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  Nobre et al., 2017;Huntingford et al., 2014). This agrees with the strong negative and positive anomalies respectively in d and θ −1 , discussed above. Indeed, 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-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), proving 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 monovriate precipitation or wind extremes, are highest in Northern/Western Europe, although for wind significant α anomalies span virtually the whole continent (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 systems 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 3c-f were respectively repeated by using daily wind maxima observations ( Figure S4) and a different definition of extremes (see 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 north-eastern part of the domain

| 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 respectively in Eastern North America and Europe. 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  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 (Faranda et al., 2019b;Messori et al., 2017;Scher and Messori, 2018). 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 hot-spots of compound extremes identified. There is also the hope that this methodology may help investigating outstanding scientific questions; for example the role of Arctic Amplification in driving northern hemisphere mid-latitude climate extremes (e.g. Cohen et al., 2014;Barnes, 2013;Coumou et al., 2018;Screen and Simmonds, 2014;Mori et al., 2014;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, 10m wind gusts, 2m 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 north-eastern 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) Fig. 4d). The North American region is, however, shifted to the north-east 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. Messori et al., 2016;Harnik et al., 2016;Overland and Wang, 2017;Trenary et al., 2015). 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 Sect. 5.2.
The spatial information on the extreme events intrinsic to our approach can thus provide new insights into multihazards (or compound events) research. Our dynamical systems perspective is flexible, links directly to concepts of predictability and coupling and may be useful in both practical (e.g. numerical forecasting of extremes) and more theoretical (e.g. understanding the atmospheric drivers of extremes) contexts.

AU T H O R C O N T R I B U T I O N
GM and PDL conceived and designed the study. PDL performed the analyses, created the figures and wrote the first manuscript draft. GM and DF provided the MATLAB code to compute the dynamical systems metrics and contributed to the methods. FP designed and performed the statistical tests. All the authors contributed to preparing the manuscript.

C O N FL I C T O F I N T E R E S T
The authors declare that they have no conflict of interest.

F I G U R E S
F I G U R E 1 (a) Daily climatological means of the precipitation-wind co-recurrence ratio (α). (b) Cumulative distribution functions (CDFs) of local dimension (d ) anomalies associated with α daily extremes (> 99 t h quantile); (c) 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). Blue shaded areas in (a) represent one standard deviation from the mean. The data cover 1979-2018 over Europe.  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 < 1 st quantile.