The role of atmospheric rivers in the distribution of heavy precipitation events over North America

. Atmospheric rivers (ARs) are ﬁlaments of extensive water vapor transport in the lower troposphere that play a crucial role in the distribution of freshwater but can also cause natural and economic damage by facilitating heavy precipitation. Here, we investigate the large-scale spatiotemporal synchronization patterns of heavy precipitation events (HPEs) over the western coast and the continental regions of North America (NA), during the period from 1979 to 2018. In particular, we use event synchronization and a complex network approach incorporating varying delays to examine the temporal evolution of spatial patterns of HPEs in the aftermath of land-falling ARs. For that, we employ the SIO-R1 catalog of ARs that landfall on the western coast of NA, ranked in terms of intensity and persistence on an AR-strength scale which varies from level AR1 to AR5, along with daily precipitation estimates from ERA5 with a 0 . 25 ◦ spatial resolution. Our analysis reveals a cascade of synchronized HPEs, triggered by ARs of level AR3 or higher. On the ﬁrst 3 d after an AR makes landfall, HPEs mostly occur and synchronize along the western coast of NA. In the subsequent days, moisture can be transported to central and eastern Canada and cause synchronized but delayed HPEs there.

With this increased interest in understanding the dynamics, impacts, and future evolution of ARs, a plethora of methodological approaches to identify and track these atmospheric features have been proposed, and multiple AR catalogs have been produced and made available to the public (Gershunov et al., 2017;Guan and Waliser, 2015;Prabhat et al., 2021;Pan and Lu, 2019;Traxl, 2022). This abundance of information has brought big challenges to the climate research community as the climatological statistics of ARs have proven to be highly dependent on the identification method chosen (Huning et al., 2017), affecting, in particular, the resultant AR climatologies and the attribution of high-impact weather and climate events to ARs (Shields et al., 2018). As a collective effort to address these issues, the Atmospheric River Tracking Method Intercomparison Project (ARTMIP) has quantified and analyzed the uncertainties in AR science based on the choice of detection/tracking methodology (Shields et al., 2018;Rutz et al., 2019;O'Brien et al., 2022;Lora et al., 2020). The results achieved by this scientific community provide guidelines to identify the most appropriate algorithm for a given scientific question or region of interest.
As a consequence, novel and relevant topics of AR science have been studied, such as the initiation and evolution of ARs, their moisture sources Waliser and Guan, 2017;Rutz et al., 2014), and the foreseen response of ARs to a warmer or different climate (Gao et al., 2015;Hagos et al., 2016;Payne et al., 2020). Due to these contributions, the following key findings have been recently revealed to the climate scientific community: 1. In the Northern Hemisphere, ARs usually originate in the midlatitude ocean basins  paired with extratropical cyclones (Zhang et al., 2019).
2. As they transit to the east, ARs accumulate and transport moisture, primarily to the western coasts of NA and Europe, where they facilitate precipitation and play a key role both in the freshwater supply and in the occurrence of HPEs.
3. In the context of ongoing climate change and as a response to the higher water vapor content in a warmer atmosphere, a poleward shift in the land-falling location, together with an increase in the frequency and intensity of ARs, can be expected in the coming decades Hagos et al., 2016;Payne et al., 2020).
4. Compared with the present, ARs affecting middle and high elevations are expected to result in more liquid than solid precipitation, exacerbating the potential risk and severe impacts of natural hazards such as floods and landslides .
In light of the scientific knowledge that has been gained in recent decades about ARs, there has been an increasing effort to characterize and predict the landfall of ARs along the North American west coast by presenting comprehensive analyses of their drivers and properties. However, the spatiotemporal synchronization patterns of HPEs induced by ARs have not yet been studied. Here, we understand spatiotemporal synchronization as a relation between pairs of precipitation time series measured at different locations, for which events in one of the time series are significantly followed or preceded by events in the other one. Such an assessment has led, among other findings (see, e.g., Boers et al., 2013;Stolbova et al., 2014;Agarwal et al., 2019;Wolf et al., 2020b), to forecasting precipitation events in the eastern Central Andes (Boers et al., 2014a) and identifying Rossby waves as one controlling mechanism of HPEs worldwide (Boers et al., 2019). In this light, it has not been examined to what extent ARs are accompanied by characteristic synchronization patterns of HPEs. Additionally, the lag-dependent spatial impact of ARs making landfall on the western coast of NA remains unrevealed. To elaborate on the issues exposed by the ARTMIP project (Shields et al., 2018;Rutz et al., 2019;O'Brien et al., 2022), in this study we address these research gaps using two AR catalogs based on different tracking schemes. Both catalogs considered here, the SIO-R1 product which was recently published by Gershunov et al. (2017) and a self-constructed one using the Image-Processing-based Atmospheric River Tracking (IPART) algorithm (Xu et al., 2020;Traxl, 2022), cover the period between 1979 and 2018. Based on the occurrence of ARs, we perform time series and complex network analyses evaluating the nonlinear spatiotemporal correlation of HPEs and their relation to ARs. To interpret our results, we furthermore study the synoptic conditions of the lower atmosphere during AR events leading to the identified synchronization patterns.
The paper is structured as follows. First, we introduce the employed data sets and methods, in particular the characteristics of the two AR catalogs, ERA5, and the event synchronization (ES) and complex network techniques. Second, we conduct an ES-based assessment of the temporal correlation between land-falling ARs and HPEs for different lags. Having revealed different timescales at which AR-related HPEs occur, we set up two climate networks based on HPEs taking place at different lags. Finally, we study composite anomalies of integrated water vapor transport, geopotential height, wind, and precipitation for the times during which we iden-tified features of synchronized HPEs to discuss our findings in the context of the guiding climatology.
2 Data and methods

Data sets
We use data from ERA5 (Hersbach et al., 2020;Hersbach et al., 2023a, b). All ERA5 data sets are available on a longitude-latitude grid with a spatial resolution of 0.25 • × 0.25 • . We construct daily estimates for integrated water vapor transport (IVT), geopotential height at 500 hPa, wind at 650 hPa, and precipitation by considering the daily mean of the hourly data sets for the period 1979-2018. To examine the synchronization of HPEs, we especially consider the 95th percentile thresholds of the daily precipitation estimates (see Supplement, Fig. S1). Only days exceeding 1 mm of total precipitation, which we refer to as wet days, are used for computing the percentiles.
Recent studies have revealed the biases present in ERA5, especially for precipitation estimations. Disagreements on the number of wet days, the co-occurrence of precipitation events, and the precipitation intensity were identified, along with a consistent pattern of decreasing agreement with increasing intensity of events, independently of the season (Rivoire et al., 2021). Larger differences were found over western NA, where ERA5 has additional difficulties in detecting and estimating orographic precipitation events (Adhikari and Behrangi, 2022). However, despite these biases, we use ERA5 because it provides a globally gridded hourly product, from which we extracted regional daily estimates of all the climatic variables considered, maintaining consistency among the data and with the large-scale circulation patterns.
In addition to the ERA5 data set, we use the SIO-R1 catalog of ARs by Gershunov et al. (2017). It includes ARs land-falling on the western coast of NA and was constructed using Lagrangian backtracking of high values of two variables, namely the vertically integrated horizontal vapor transport and the vertically integrated water vapor (IWV), on a longitude-latitude grid with a resolution of 2.5 • × 2.5 • . The catalog features a 6-hourly time series indicating whether an AR has been active as well as the grid cells covered by the AR and the IVT over the grid cell along the coast where the AR made landfall. We transform the 6-hourly catalog into a daily one by considering each day with at least one of the four 6-hourly time steps with an active AR as an AR day. Approximately one-third of all days of the analysis period are AR days (at least one AR active somewhere in the spatial domain covered by the catalog). These days are distributed relatively equally over the different years but are strongly seasonally clustered in the boreal autumn and winter. Furthermore, we create an additional catalog of ARs with features similar to the SIO-R1 catalog but using the IPART algorithm (Xu et al., 2020). As opposed to conventional detection methods that rely on thresholding of IVT or IWV fields (for instance the detection algorithm of the SIO-R1 catalog), IPART implements the detection task from a spatiotemporal-scale perspective and is, therefore, free from magnitude thresholds. The advantage of IPART's approach is that it negates the implicit assumption of thresholding approaches that the atmospheric moisture level stays unchanged throughout the analysis period. As input to the IPART algorithm, we use IVT fields of the ERA5 data set re-gridded to a spatial resolution of 0.75 • × 0.75 • and a temporal resolution of 6 h. The parameters passed to the different steps of the IPART algorithm are summarized in Table S1 in the Supplement (Traxl, 2022). We transform the 6-hourly product into a daily one in the same manner as described for the SIO-R1 catalog.
To separate the impact of rather weak ARs from strong ARs, we rank AR events of both catalogs in terms of intensity and persistence on the AR-strength scale proposed by Ralph et al. (2019) and using the notation in Eiras- Barca et al. (2021). As a result, each AR event is assigned a level that increases from AR1 to AR5, and we run our analysis repeatedly, excluding ARs from lower levels.

Directed event synchronization (ES) between land-falling ARs and HPEs
ES, as initially introduced by Quiroga et al. (2002), is a nonlinear temporal correlation measure that quantifies the covariability of events in a pair of time series. While it was originally proposed for analyzing spike trains in electroencephalographic time series, ES is nowadays an established tool employed to construct climate networks (Malik et al., 2012;Boers et al., 2013Boers et al., , 2019Ozturk et al., 2019;Wolf et al., 2020). To define directed ES, let the sequence {t µ l } µ=1,...,n l denote the time series of ARs of level l or higher, making landfall on the western coast of NA, and let the sequence {t ν i } ν=1,...,n i denote the time series of HPEs observed at grid cell i. The total number of land-falling ARs of level l or higher and the total number of HPEs at grid cell i are denoted by n l and n i , respectively. We say that the HPE ν, observed at location i, at time t ν i , is synchronized with the preceding AR event µ, of level l or higher, which made landfall at time t µ l , if and only if their temporal delay, t ν i − t µ l ≥ 0, does not exceed the dynamical delay: This adaptive lag allows us to consider events in more densely and more sparsely occupied parts of the time series in an automated manner, in contrast to the classical lead-lag approach that only allows one lead or lag for the entire time series. Since the first and last events of each time series do not have a preceding or following event, we exclude them from our calculations and only consider µ = 2, 3, . . ., n l − 1 and ν = 2, 3, . . ., n i − 1. See Fig. 1 for an illustration of the computation of the dynamical delay τ µ,ν l,i . To avoid a collapse of the dynamical delay to τ µ,ν l,i = 1 2 time step due to sequences of consecutive events, we only consider the first event of all event sequences (cluster-corrected ES) (Boers et al., 2019;Wolf and Donner, 2021). Additionally, we can limit the dynamical delay by a minimal value τ min to consider a minimum lag between synchronized events and also by a maximal value τ max to prevent an unrealistically large temporal delay between synchronized events.
Then, the synchronization condition reads and for grid cell i, we define as the total number of HPEs that can be uniquely associated with a preceding land-falling AR in a time-resolved manner, within an interval of at least τ min days and no more than τ max days.
Finally, we analyze the statistical significance of each empirical value ES l,i by means of a null model. To account for the time order of directed ES, we use the unchanged time series of land-falling ARs and incorporate 1000 surrogate time series of HPEs, preserving the original number of events n i but destroying a potential correlation structure. We calculate the value of ES l,i from the time series of land-falling ARs and each surrogate to estimate an empirical probability distribution, which we then use to infer the significance level of our ES l,i value (Boers et al., 2019). We say that ES l,i is significant at the level where denotes the Heaviside function and ET ρ l,i is the ρth percentile of the surrogate test distribution for ES l,i .

Identification of ARs highly synchronized with a specific region
One particular advantage of ES is that we can use it to identify the temporal ordering and the time delay of synchronized events. From the synchronization condition (see Eq. 2), note that S µ,ν l,i = 1 if and only if the land-falling AR event µ precedes and is synchronized with the HPE event ν observed at location i. In that case, the time delay between this pair of uniquely associated events is d µ,ν l,i = t ν i −t µ l . If we use Eq. (4) to only consider the grid cells where the synchronization between land-falling ARs and HPEs is significant, then, for a region of interest R, we can define as the total number of HPEs within R that were preceded and uniquely associated with the land-falling AR µ at the significance level α = 1 − ρ 100 during the time window [τ min , τ max ] (Boers et al., 2019). Based on this definition, we can retrieve the time delays between the AR event µ and the significantly synchronized HPEs, (6) to define and calculate the typical synchronization delay between the AR event µ and the region of interest R as the mode of D µ→R (if there are multiple modes, we take the smallest one).
We use this framework to select the ARs of level l or higher with the largest number of significantly synchronized HPEs within the region of interest R and to identify the time points that are then used to compute the composite anomalies of integrated water vapor transport, geopotential height, upper-level meridional wind, and precipitation to be shown in our results.
Even more, we can use Eq. (5) to identify ARs that were not synchronized with HPEs in the region of interest (ES ρ µ→R = 0). Based on a precedence analysis, we can also identify the dates when HPEs occurred in the region of interest without any AR making landfall on the coast during the previous 12 d (the selection of this preceding time window comes after one of our results). These time points are used to compute composite anomalies of the aforementioned climatological variables in order to reveal the particular synoptic conditions that differentiate ARs distributing synchronized HPEs into the region of interest.

Directed event synchronization (ES) between HPEs at different locations
We are also interested in investigating if there is a directed synchronization pattern between HPEs at different locations in the aftermath of land-falling ARs. Adapting the definition of directed ES in Sect. 2.2, we consider two HPEs ν and ϕ in time series describing observations made at grid cells i and j at times t ν i and t ϕ j as synchronized if and only if their temporal delay, t ϕ j −t ν i , does not exceed the dynamical delay: Again, sequences of consecutive events are counted as single events (cluster-corrected ES), and the first and last events of each time series are discarded; i.e., ν = 2, 3, . . ., n i −1 and ϕ = 2, 3, . . ., n j −1, where n i and n j denote the total number of HPEs at grid cells i and j , respectively.
In this case, when the event at i happens before the event at j , the synchronization condition reads Note the subtle but important difference with Eq. (2): we do not include events that occur simultaneously at different locations, since we cannot determine their temporal ordering.
We define the directed event synchronization from i to j as which is the total number of synchronized events where an event at i precedes an event at j by at least τ min days and no more than τ max days. The reverse time direction is given by resulting in the asymmetric matrix ES, which we use for setting up climate networks as described in the following section.

Climate networks
Functional networks are defined as graphs, where nodes represent the elements of a complex system and edges represent the interaction between them. In functional networks, edges are placed between nodes following some statistical similarity, regardless of whether the nodes are physically connected or not. A climate network, as introduced in previous work by Tsonis and Roebber (2004), Tsonis et al. (2006), and Donges et al. (2009a), is a functional network whose nodes are identified with geographical locations and whose edges account for a significant and strong correlation between the climatological time series measured at the respective locations. Recently, the climate network approach has attracted much attention after being successfully applied to reveal novel insights into the dynamics of the Earth's climate system over different spatiotemporal scales (Tsonis and Swanson, 2008;Yamasaki et al., 2008;Donges et al., 2009bDonges et al., , 2011Malik et al., 2012;Steinhaeuser et al., 2012;Boers et al., 2014bBoers et al., , 2015Agarwal et al., 2019;Boers et al., 2019;Messier et al., 2019;Wolf et al., 2020b). A network consists of N nodes connected by e edges. The topology of such a network is commonly encoded in the adjacency matrix A, with elements a ij indicating whether nodes i and j are connected. In this study, we construct climate networks based on directed ES to assess the spatiotemporal correlation structure of HPEs in NA and to unravel possibly non-linear and long-range climatic linkages associated with the land-falling of ARs on the western coast of NA. We identify the nodes of the network with the gridded time series of the ERA5 data and connect the nodes based on their statistical association evaluated by directed ES. To transform the daily ERA5 precipitation estimates to an event time series, we threshold the time series of wet days of each grid cell at the 95th percentile. Subsequently, we apply the clustercorrected ES (Boers et al., 2014a;Wolf et al., 2020) with a lower and upper threshold for the dynamical delay which is specified in the respective sections. We use the resulting synchronization matrix ES to place a directed link pointing from grid cell i to grid cell j if HPEs at i precede and are synchronized with HPEs at j .
Note that, by construction, we typically have different numbers of events at different grid cells, since we calculate the percentile of the wet days and not of the entire time series. The number of events affects the measured values of ES; therefore, to ensure that the edges placed between nodes are statistically significant and to account for the time order of directed ES, we apply a locally tailored significance testing scheme. For a pair of grid cells (i, j ), we set up a null model using the unchanged time series of HPEs for grid cell i and 1000 surrogate time series of HPEs for j that preserve the respective number of events n j . We assume that the events at j occur independently according to a uniform random distribution and compute the value ES i,j for each surrogate time series, obtaining an empirical distribution function. We then connect nodes in the network (by setting a ij = 1 in the adjacency matrix) if ES i,j exceeds the 99.5th percentile of the respective surrogate test distribution (Boers et al., 2019): where denotes the Heaviside function, ET 99.5 i,j is the 99.5th percentile of the surrogate test distribution for ES i,j , and Kronecker's delta δ ij is used to exclude self-loops.
As a remark, we want to emphasize that considering the temporal ordering of the events to calculate ES (see Eqs. 9 and 10) results in a directed climate network, for which the adjacency matrix is not symmetric (in contrast to undirected networks). Moreover, since we also do not consider edge weights, the adjacency matrix is binary: a ij = 1 if an edge points from node i to node j ; a ij = 0 otherwise. This methodology of combining ES and complex networks is based on the idea that ARs influence the way HPEs synchronize at different locations and that these effects are contained in the internal structure of the climate network, which can be accessed by appropriate complex network measures. In many cases, climate networks are a powerful alternative to more traditional approaches based on eigenvalue techniques (e.g., principal component analysis, PCA) (Boers et al., 2013). However, the methodology employed here has been specifically developed to analyze time series of climatological extreme events, for which PCA-like methods are not applicable due to the binary-like structure and the non-Gaussian distributions of the data (Malik et al., 2012;Boers et al., 2013Boers et al., , 2014aBoers et al., , c, b, 2015Boers et al., , 2019Stolbova et al., 2014).
To calculate the number of nodes to which node i is connected, we compute the in-degree k in i (out-degree k out i ) as the total number of links pointing to (from) grid cell i (Boccaletti et al., 2006). To aggregate both measures and to highlight regions of predominately outgoing (or incoming) connections, we define the network divergence as with positive (negative) values of d i indicating sources (sinks) of the network: HPEs in these locations are followed (preceded) by HPEs in other locations.
3 Results and discussion

HPEs synchronized with strong ARs
As a first step, we investigate where HPEs occur synchronized with land-falling ARs and at which lags. For that, we employ ES and evaluate the synchronization between the AR time series and the time series of HPEs at each grid point in the study area. To obtain the latter, we threshold the precipitation time series during wet days at the respective 95th percentile. For the former, we consider ARs from the SIO-R1 catalog making landfall on the North American west coast at a latitude ≥ 47.5 • N. Initially, we included all ARs, but additional analyses showed that our results are predominantly caused by ARs that landfall north of 47.5 • N (see Supplement, Fig. S2). Moreover, we want to emphasize that our results can be reproduced using an alternative AR catalog based on the IPART algorithm (also featured in the Supplement; see Fig. S3). To separate the impact of rather weak ARs from strong ARs, we differentiate between AR levels (classification based on Ralph et al., 2019, using the notation in Eiras- Barca et al., 2021) and run the analysis repeatedly, excluding ARs from lower levels. Figure 2 shows the grid points whose time series of HPEs are significantly synchronized with the AR time series, given a particular parameter setting of ES and ARs of level AR3 or higher. First, note that when τ min = 0 (left column), large areas close to the western coast of NA show significant correlations. This pattern is caused by HPEs on the coast that are directly triggered by ARs and is not affected by increasing τ max . If events occur in close succession, then a higher possible maximal delay will often not be considered and, therefore, the pattern does not change visibly. This strong synchronization between land-falling ARs and HPEs on the western coast of NA was expected as it has already been implied by findings of previous studies (Neiman et al., 2008;Gershunov et al., 2017;Waliser and Guan, 2017;Ralph et al., 2019), and serves as a proof of concept for our methodology. Excluding ARs of the lower levels AR1 and AR2 does not change this result (see Supplement, Fig. S4). Therefore, we have selected ARs of level AR3 or higher for our analysis.
When τ min = 3 (second column), the synchronization close to the coast decreases as most HPEs occur on the first days after an AR makes landfall. Additionally, most ARs do not persist longer than 3 d (Gershunov et al., 2017). The remaining synchronized events are likely associated with ARs  . Event synchronization (ES) between ARs making landfall on the western coast of NA and HPEs. We use the SIO-R1 catalog of land-falling ARs but only consider ARs making landfall north of 47.5 • N. ES is calculated with τ min = 3 and τ max = 12. From (a) to (e) the lower limit of the considered AR level increases: (a) ARs of level AR1 and higher, i.e., all ARs, (b) ARs of level AR2 and higher, etc. Color bar as in Fig. 2. of the higher levels, which have a longer persistence. Additionally, we observe a patch of synchronized events in central and eastern Canada. This pattern is strongest when τ min = 3 and τ max = 6 and stands out to a smaller extent up to τ max = 12. For elevated values of τ min (third and fourth columns), the synchronization pattern completely vanishes.
This result implies that in central and eastern Canada, HPEs occur synchronized with (but lagging) ARs making landfall on the western coast of NA. We, therefore, suspect that moisture that has been transported to the North American west coast by an AR can be channeled to central and eastern Canada and also cause HPEs there, as suggested in previous studies by Rutz et al. (2014), Waliser and Guan (2017), and Guan and Waliser (2019).

Synchronization across AR strength
Using ES, we have identified a region of synchronized HPEs in central and eastern Canada, as explained in the previous section. To further evaluate how the results depend on the selected AR-level criteria, we exclude ARs of the lower levels from the analysis stepwise, as shown in Fig. 3. Doing this, we reduce the number of events in the AR time series, which heavily affects the outcome of ES. As a result, we find that when we consider all ARs, the signal in central and eastern Canada is present but is accompanied by a more prominent synchronization pattern right next to the western coast of NA (Fig. 2). When we discard low-level ARs, the pattern next to the coast is filtered out, and when we consider ARs of level AR3 or higher, just the signal in central and eastern Canada is left. Only examining ARs of levels AR4 and AR5 leads to a vanishing of the synchronization in central and eastern Canada, although precipitation anomalies show HPEs in that region for composites based on the days after such strong ARs (see Supplement, Fig. S5). This is a result of reducing the number of events in the AR time series. Very intense ARs are rare, and when we only consider them, the AR time series contains very few events. In other words, the AR time series gets too sparse and the ES scores are not significant anymore, although HPEs might be always caused by such ARs but likely not just by them.

Network analysis of HPEs during AR days
In the previous sections, we examined the synchronization between HPEs and the singular AR time series. To elaborate further on the concept of synchronized HPEs, we assess how precipitation at different locations is organized during AR days. For that, we select the days with active ARs of level AR3 or higher, making landfall north of 47.5 • N and the respective subsequent 3 d. We acknowledge that with this approach we can only relate HPEs close to the coastline to ARs. Based on these selected days, we run a network analysis using ES with the parameters τ min = 0 and τ max = 3. There- Figure 4. Network divergence based on event synchronization (ES) between HPEs at different locations during AR days. We calculate ES with τ min = 0 and τ max = 3, and we only consider HPEs that occurred from 0 to 3 d after an AR of level AR3 or higher makes landfall north of 47.5 • N according to the SIO-R1 catalog. In the network, only nodes with significantly directed event synchronization are connected (see Sect. 2.5 for more details). Purple (orange) colors indicate regions with positive (negative) divergence, i.e., nodes with more outgoing (incoming) connections. fore, we investigate the immediate synchronization pattern of HPEs occurring simultaneously with ARs.
In Fig. 4 we show the resulting network divergence, which is characterized by a large area of negative values on the coast overland and positive values over the eastern part of the Pacific. Note that network divergence is computed by subtracting k out i −k in i (out-degree minus in-degree; see Eq. (12)); therefore, areas with positive (negative) values have more outgoing (incoming) edges and can be regarded as sources (sinks) of the network. We identify a clear source in the Pacific Ocean and a sink close to the western coast of NA. We find that HPEs over the Pacific occur first and are followed by HPEs over the western parts of NA. Additionally, there is a wave-like pattern hidden in Fig. 4 (moving further east, we observe another large area of positive network divergence, followed by negative values over the Atlantic and the east coast of NA). Since our filtering only allows for interpreting the dynamics near the western coast of NA, we can only speculate about the causes of this pattern. We suggest that this wave pattern resembles a large-scale midlatitude wave train that may indeed contain alternating moist and dry advection which may include a sequence of synchronized landfalling ARs and HPEs (Mo and Lin, 2019).

Network analysis of HPEs in the aftermath of AR events
After the proof of concept in the previous section, we now extend the spatial and temporal domains. Our initial analysis showed that HPEs that occur in central and eastern Canada are synchronized with ARs land-falling on the western coast of NA considering delays in a window between τ min = 3 and τ max = 12 (See Fig. 2). Now, to link HPEs that are related to ARs, we choose the following setup: we consider events that occurred from 0 to 12 d after the landfall of an AR of at least level AR3 and employ ES with τ min = 3 and τ max = 12.
With that, we assure that we keep HPEs on the coastline and in central and eastern Canada but only allow synchronization for temporal delays larger than 3 d. Consequently, we avoid obtaining a strong signal of synchronized events along the western coast of NA (where the main synchronization of events happens during the first 3 d after the first AR day). In other words: we examine the delayed synchronization pattern of HPEs (at least 3 d between events) occurring at any time after an intense AR makes landfall. The resulting network divergence is displayed in Fig. 5a. We identify a region of positive network divergence along the northern part of the western coast of NA and especially a region of reduced network divergence over central and eastern Canada, where the synchronization between the AR time series and the HPEs was initially discovered. To finally verify that there is a strong connection between the North American west coast, where we find a large number of outgoing edges, and central and eastern Canada, where many edges terminate, we analyze where edges that connect to central and eastern Canada originated (see out-degree and red box in Fig. 5b). Figure 5b highlights the grid cells where edges originate that terminate in central and eastern Canada (red box, 58 to 68 • N and 95 to 122 • W). This box has been chosen based on where we have found the synchronization between the AR time series and the HPEs time series (see Figs. 2 and 3). One main source of edges to this region is the north and northwest of the study area. This confirms that edges emerge from the region that is marked by positive values in the network divergence (where HPEs synchronized within the first 3 d of the landfall of an AR) and terminate in the red box.
In summary, we have identified a cascade of HPEs: on the first 3 d after an AR of level AR3 or higher makes landfall, HPEs mostly occur close to or on the coast and synchronize in this area. In the subsequent days, moisture can be transported to central and eastern Canada and cause HPEs there. This takes place between 3 and 12 d after the first ARinduced precipitation on the coastline.

Synoptic conditions facilitating AR-induced HPEs in central and eastern Canada
Measuring the synchronization between time series revealed the spatial extent as well as the temporal dimensions of heavy precipitation related to land-falling ARs over NA. A delayed synchronization pattern between ARs making landfall on the western coast of NA and HPEs in central and eastern Canada was identified. To examine the climatic drivers leading to this long-range correlation, we study the synoptic conditions of different climatological variables for three types of events, defined as follows: 1. Type I event. An AR makes landfall on the western coast of NA and synchronizes with HPEs in central and eastern Canada in the subsequent 3 to 12 d.
2. Type II event. An AR makes landfall on the western coast of NA but does not synchronize with HPEs in central and eastern Canada in the subsequent 3 to 12 d.
3. Type III event. HPEs occur in central and eastern Canada but no AR made landfall on the western coast of NA during the previous 12 d.
To carefully choose the time points corresponding to each type of event, we first identify specific times with high event synchronization between land-falling ARs and HPEs in central and eastern Canada. We do so by using Eq. (5) with l = 3, ρ = 0.9, τ min = 3, and τ max = 12 (to match the spatial pattern found in Figs. 2 and 3) and the region of interest as the box in Fig. 5b, which we denote as region B (see Fig. 6a). Note that the resulting {ES 0.9 µ→B } µ=1,...,n l is a sequence that gives the total number of HPEs in region B that were preceded and uniquely associated with the AR event µ at a significance level α ≤ 0.1. We identify ARs whose total number of associated HPEs is above the 80th percentile of the nonzero values of this sequence, and we get 35 AR events highly synchronized with HPEs in central and eastern Canada. The landfalling times of these ARs are the time points of type I events. We also use Eq. (6) to determine the typical synchronization delay between region B and each of these highly synchronized ARs, with the most common value being 5 d. We then select the 35 most intense AR events for which ES 0.9 µ→R = 0 -i.e., those that did not synchronize with HPEs in central and eastern Canada -and we select their land-falling times as the time points of type II events. Finally, we identify the time points of type III events as the 35 d with the highest number of HPEs in central and eastern Canada that occurred in the absence of land-falling ARs during a time window of 12 preceding days.
We use these three types of events to analyze the antecedent IVT over the western coast of NA and the subsequent precipitation over central and eastern Canada and thus define the two regions of interest shown in Fig. 6a. Region A, where we analyze IVT anomalies, covers the area where the ARs make landfall north of 47.5 • N (from 46 to 59 • N latitude and from 223 to 238 • E longitude). Region B, where we study precipitation anomalies, delimits the area where the delayed synchronization pattern between ARs and HPEs was identified (from 58 to 68 • N latitude and from 238 to 265 • E longitude). In Fig. 6b we show the annual cycle at daily resolution of the mean IVT in region A (solid dark line). The shading encloses 1 standard deviation from the daily mean and the dashed line indicates the multi-annual mean IVT of the region. The dots show the IVT values of the three types We calculate ES with τ min = 3 and τ max = 12 and consider all HPEs that occurred from 0 to 12 d after an AR of level AR3 or higher makes landfall north of 47.5 • N according to the SIO-R1 catalog. In the networks, only nodes with significantly directed event synchronization are connected (see Sect. 2.5 for more details). (a) Network divergence. Color bar as in Fig. 4. (b) Out-degree of a directed network. Highlighted are the grid cells where edges terminating in central and eastern Canada originated (red box, 58 to 68 • N and 95 to 122 • W). This box was chosen based on the region with significant synchronization between land-falling ARs and HPEs found in Figs. 2 and 3. First, note the clear imprint of ARs on IVT values for type I and II events (Fig 6b, blue and red dots). Whether or not the land-falling AR synchronizes with HPEs in central and eastern Canada, the IVT values in region A are above the daily and the multi-annual mean with very few exceptions. Highly synchronized ARs, corresponding to type I events, have particularly anomalous values of IVT, in most cases exceeding the mean climatology by 1 standard deviation or more. By contrast, the IVT values of most type III events (orange dots) are around or below the daily and multi-annual mean IVT of the region, confirming the absence of a preceding AR making landfall on the coast for this type of event. The delayed precipitation over central and eastern Canada associated with each type of event is displayed in Fig. 6. For type I events, as expected, the precipitation values are above the daily and the multi-annual mean, with anomalies that frequently exceed the mean climatology by more than 1 standard deviation. Conversely, type II events have average or lower precipitation values, which is consistent with the absence of a synchronization pattern between land-falling ARs and HPEs. Not surprisingly, type III events referring to HPEs with no preceding ARs also have high precipitation anomalies.
Besides this proof of concept of the effectiveness of our methodology to identify types of events and synchronization delays, Fig. 6 also reveals a key characteristic of the ARs that are highly synchronized with HPEs in central and eastern Canada: their seasonality. These ARs, which define the type I events, occur most likely during July, August, and September and can be characterized as intense, long-lasting, late-summer ARs. Type II events that do not contribute to the observed synchronization pattern in central and eastern Canada are caused by ARs with similar intensity and persistence but making landfall during the early winter, i.e., in October, November, and December. Lastly, type III events describing HPEs that are not preceded by ARs are more common during the early summer months of June, July, and August.
To reveal the synoptic conditions facilitating the delayed effect of ARs in the precipitation over central and eastern Canada, we compute composite anomalies of IVT, geopotential height at 500 hPa, wind at 500 hPa, and precipitation on the days of type I events (when the highly synchronized ARs made landfall) and for the following 3, 5, and 12 d. The results are shown in Fig. 7. Similar figures for type II and III events can be found in the Supplement (see Figs. S6 and S7).
In Fig. 7, we first present the temporal evolution of the IVT anomalies (top row) after the landfall of highly synchronized ARs on the western coast of NA. The high positive IVT anomaly on the Pacific when lag = 0 d is a clear characteristic of land-falling ARs. It is noteworthy that, in the following days, this anomalous water vapor influx is able to penetrate the continent through the topographic gap previously identified by Rutz et al. (2014) and that it traverses the mainland reaching central and eastern Canada, where the synchronization pattern was found. In the second row of Fig. 7, the geopotential height at 500 hPa is assessed for varying delays as for the IVT. The trough together with the strong negative anomalies in the northwest of the study region at the moment of landfall indicates the position of the cold front driving the ARs. In the following days, the cold front digs into Canada as a ridge builds in the central and eastern United States. This configuration of the geopotential height field during highly synchronized ARs implies a midlevel pressure dipole that traverses the continent accompanied by a southwesterly steering wind (third row of Fig. 7), bringing warm moist air from the west coast into the northern regions of NA. Under these circumstances, high precipitation anomalies occur just downwind of the trough: first at the coast and then over region B, where we identified the delayed synchronization pattern between ARs and HPEs (fourth row of Fig. 7). These synoptic conditions, which are exclusive to highly synchronized ARs (see Figs. S6 and S7) and have already been identified as conducive to the occurrence of seasonal precipitation extremes over Canada (Tan et al., 2019), explain the physical mechanism by which land-falling ARs serve as moisture sources of HPEs in the northern regions of NA.

Conclusions
In this study, we have investigated the influence of ARs on the large-scale spatiotemporal synchronization patterns of HPEs over NA. For this purpose, we have first analyzed whether there is a significant association between ARs making landfall on the western coast of NA and HPEs over the coastal and continental regions. Employing event synchronization (ES), we have revealed timescale-dependent spatial patterns of HPEs that are significantly correlated with ARs making landfall north of 47.5 • N: (i) immediately after an AR makes landfall on the coast, HPEs synchronize over the coastal areas; (ii) then, from 3 d after the landfall, the synchronization close to the coast decreases significantly, and only HPEs associated with more persistent ARs remain; and (iii) from 3 to 12 d after an AR makes landfall, a synchronization pattern between land-falling ARs and HPEs in central and eastern Canada emerges. These results have been reproduced using an alternative AR catalog, establishing the robustness of our findings.
After examining the synchronization of HPEs with the time series of land-falling ARs, we have analyzed the organization of HPEs on the day of landfall and the subsequent days. For this, we have evaluated directed ES between time series of HPEs at different locations and for different temporal lags after the landfall of ARs. Based on that, we have first constructed a complex network considering a time window from 0 to 3 d after landfall. This result confirmed the common knowledge of the association between ARs and HPEs in western NA. Initially, HPEs occurring simultaneously with ARs are synchronized over the eastern Pacific Ocean. They are then followed by synchronized HPEs on the western coast of NA. By examining a second complex network based on HPEs occurring at any time after the landfall of an intense AR but only allowing synchronization with a delay from 3 to 12 d, we have uncovered a strong connection between HPEs Figure 7. IVT, geopotential height at 500 hPa, wind at 500 hPa, and precipitation anomalies (from top to bottom) from 0, 3, 5, and 12 d (from left to right) after the landfall of ARs leading to a delayed synchronization pattern with HPEs in central and eastern Canada. Only ARs of level AR3 or higher with land-falling latitude north of 47.5 • N are considered. In the second row, the shading indicates the anomaly of the geopotential height at 500 hPa and the contours show the mean geopotential height. In the third row, the shading indicates the anomaly of the meridional wind at 500 hPa and the arrows show the mean wind field. on the North American west coast and in central and eastern Canada: moisture from ARs land-falling along the coast can be transported to central and eastern Canada and cause HPEs there.
To further investigate this result, we identified specific days with high event synchronization between land-falling ARs of level AR3 or higher and HPEs in central and eastern Canada. Then, we used these time points to analyze the composite anomalies of vertically integrated water vapor transport (IVT), geopotential height at 500 hPa, wind at 500 hPa, and precipitation on the day of landfall and during the subsequent 3, 5, and 12 d. Our approach yielded two key findings regarding the climatic conditions that facilitate AR-induced HPEs in central and eastern Canada: (i) intense, long-lasting, late-summer ARs making landfall north of 47.5 • N on the western coast of NA are the ones leading to the occurrence of delayed HPEs in central and eastern Canada, and (ii) such ARs are driven by a cold front digging from the northeast Pacific Ocean into Canada as a high-pressure region builds in central and eastern United States. This mid-level pressure dipole traverses the continent accompanied by a southwesterly steering wind, bringing the warm moist air deposited on the coast by the land-falling ARs into central and eastern Canada and facilitating synchronized but delayed HPEs there. These particular synoptic conditions, which are consistent with the seasonality of the identified ARs, explain the physical mechanism by which late-summer ARs serve as moisture sources of HPEs in the northern regions of NA. However, whether or not these ARs remain as identifiable objects following landfall remains an open question that requires a different catalog than SIO-R1 to be addressed, more specifically, a catalog that tracks the ARs not only until they make landfall but also as they penetrate the continent.
In summary, we have studied the spatiotemporal synchronization pattern of HPEs induced by ARs, revealing its extent and its temporal evolution. We have shown that the impact of ARs making landfall on the western coast of NA is not limited to these areas, since they can be accompanied by delayed but significantly synchronized HPEs in the continental regions. In particular, we have identified a cascade of synchronized HPEs: on the first 3 d after an AR makes landfall, HPEs occur and synchronize along the coast; on the subsequent days, this moisture can be transported to central and eastern Canada and cause synchronized HPEs there. Our results illustrate the role of ARs in the distribution of HPEs over NA, not only on the west coast but also over the continental regions through inland penetration of IVT. The findings presented in this work should be considered to better anticipate the evolution of the climate dynamics of the region and the associated impacts on the precipitation patterns in the context of a warming atmosphere, for which we expect an increased frequency and strength of the ARs as well as a northward shift in the locations where the ARs make landfall.
Author contributions. SMVB and FW conducted the analysis and wrote the paper. SMVB prepared the figures with the help of FW. DT prepared the section and table for the IPART catalog. All authors reviewed and improved the paper. SMVB took the lead in assessing the reviews, making the necessary changes, and coordinating the responses of all authors to the reviewers' comments.
Competing interests. The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. This research has been funded by the project ClimXtreme and the training group NatRiskChange. Sara M. Vallejo-Bernal acknowledges and appreciates the stimulating and fruitful discussions with Tobias Braun. Frederik Wolf wants to give special thanks to Lukas Nageln for his exceptional advice. Niklas Boers acknowledges funding from the Volkswagen Foundation. Finally, the co-authors thank the reviewers for their valuable comments and suggestions, which contributed to improving the quality of this study.
Financial support. This research has been supported by the Bundesministerium für Bildung und Forschung (grant no. 01LP1902J) and the Deutsche Forschungsgemeinschaft (grant no. GRK 2043/1).
The publication of this article was funded by the Open Access Fund of the Leibniz Association.
Review statement. This paper was edited by Nadav Peleg and reviewed by two anonymous referees.