Atmospheric river detection

ARs are identified using an updated version of an objective detection algorithm (as documented in ref.32). The algorithm uses gridded fields of positive anomalous vertically integrated water vapor transport (IVT), together with a series of intensity and geometric tests (e.g., mean intensity, total area, length, length-to-width ratio), to identify features that are of the appropriate spatial scale and are sufficiently plumelike in nature. This detection algorithm employs an occurrence-based approach (i.e., an AR occurrence is recorded for each period during which the criteria are satisfied), wherein each time step is scrutinized independently. As a result, the calculations regarding AR “hits” described in this study are based on the number of days during which AR-like conditions exist over a given landfall boundary. The updated detection algorithm used in this study does not contain the “multiple peak” logic that scrutinizes connected features within fields of anomalous IVT. We find that the mid- and high-latitude results are generally insensitive to the removal of this logic test. The majority of the results presented in this work are based on this detection algorithm applied to IVT calculated via the mass-weighted vertical integration of component winds and specific humidity from 1000–250 hPa from MERRA-2.35 Daily means are calculated and the data set is regridded to a 1.5° × 1.5° latitude–longitude grid before calculating IVT. With this data set, a static anomalous IVT magnitude threshold of ~173 kg m−1 s−1 is used to isolate features of interest. This value represents the 94th percentile of the all-season distribution of daily IVT anomaly values over the North Pacific Ocean. For the model skill comparison, this AR detection scheme is also applied to gridded IVT anomalies calculated from a set of retrospective forecasts from the ECMWF reforecast ensemble prediction system.

Atmospheric river activity

To construct time series of anomalous AR activity along the landfall boundaries, we first create a continuous time series of boolean AR “hits” for each boundary by recording a hit whenever the spatial extent of any AR feature overlaps at least one grid point of a given landfall boundary. Second, we remove the seasonal cycle of AR activity by subtracting the mean and first two harmonics calculated via fast Fourier transform applied to the calendar-day means of each AR time series. Third, we apply a 5-day running mean to each anomaly time series. Based on the premise that subseasonal predictions are founded on the presence and modulation of large-scale circulation anomalies, the running mean transitions the time series from representing only individual transient hits to capturing the larger-scale propensity of the anomalous flow pattern to influence the landfalling activity. The resulting time series are used to create the conditional composites for Fig. 2 and to train and verify the prediction scheme described throughout this work.

Predictors

Two potential sources of subseasonal predictability are used as initial conditions for the prediction scheme: the MJO and the QBO. We characterize the MJO according to the strength and location of the enhanced near-equatorial convection and the associated anomalous circulation, as determined by the components of the RMM index.36 This MJO index is a combination of two component indices, RMM1 and RMM2, representing the two leading principal components from a multivariate (equatorially averaged tropical outgoing longwave radiation and 200- and 850-hPa zonal winds) empirical orthogonal function analysis. When combined and considered in terms of their two-dimensional phase space, these component indices provide daily phase (1–8) and amplitude values (see ref.36). We consider the MJO as active when the amplitude meets or exceeds a value of one and also apply the basic logic test that the index must remain in that phase (e.g., location of active MJO signal) for at least two, but less than 20, days. See the Supplementary Information for alternative characterizations of the MJO (e.g., different index or more stringent phase event criteria). The QBO is characterized by the standardized monthly 50-hPa zonal wind index provided by the National Oceanic and Atmospheric Administration (NOAA) National Weather Service (NWS) Climate Prediction Center (CPC). We apply this index as a continuous time series, such that all months within the period of record are categorized as either EQBO (monthly mean standardized anomaly < 0) or WQBO (>0). We show the impact of introducing a more restrictive anomaly threshold, which has the possible effect of bolstering the prediction scheme, but also reducing the number of conditional samples, in the Supplementary Fig. S13. The QBO is not considered alone based on the presumption that the influence of the QBO on anomalous AR activity will primarily manifest via the modulation of the MJO’s convection and its ability to elicit an extratropical response. Hence, the predictors investigated in this study are simply the daily MJO phase (1–8) and the monthly QBO phase (EQBO or WQBO).

Empirical prediction

We generate a two-class prediction scheme from the DJFM AR anomaly time series for each region; that is, we assess the probability of being above and the probability of being below an “equal chances” 50th percentile of the time series, given the initial state of the MJO and the QBO. As constructed, the probability of being above (increased AR activity) and below (decreased AR activity) are equal when all verification dates in DJFM are considered. However, the probability distribution may shift as a function of MJO phase, QBO phase, and forecast lead. For the example shown in Supplementary Fig. S5, the distribution of 5-day average anomalous AR activity near British Columbia is shifted toward higher values 18 days following MJO phase 1 dates during EQBO conditions, relative to the mean DJFM distribution. As a result, given the initial conditions of MJO phase 1 and EQBO, the empirical scheme will predict increased AR activity for the British Columbia landfall boundary around 18 days following these initial conditions. In evaluating the two-class scheme, a prediction is considered correct when the observed response, in terms of above or below the 50th percentile threshold, from the 5-day running mean of the independent verification time series matches the predicted response; the prediction is considered incorrect otherwise. As illustrated by this example, the predictors are not explicitly weighted as they may be in a scheme based on some form of regression; in contrast, the predictors (MJO phase and QBO phase) are simply used to parse the training data in order to assess the conditional shift in the likelihood of increased or decreased AR activity relative to seasonal climatology. Though the results presented herein are based on the use of a single DJFM 50th percentile threshold, we find that the overarching conclusions (time-lagged response, patterns of increased/decreased AR activity, and so on) remain even if the threshold is allowed to vary by day-of-season. We use a leave-one-out cross-validation approach to conditionally construct and evaluate this prediction scheme.42 Specifically, the verification statistics for a given season are based on distributions constructed from historical AR activity parsed by phase of the MJO and the QBO for all DJFM seasons excluding the one “left out” verification season, ensuring independence of the verification subset. As a given season is left out, we use the training data to generate forecasts for all 121 days within the left out season and for all possible forecast leads. In so doing, we perform this leave-one-out procedure 36 times, each time leaving out just one DJFM season from the available MERRA-2 record. Because the training periods differ during the leave-one-out process, the 50th percentile threshold is recalculated each time; however, the threshold value is nearly unchanged throughout the cross-validation. The output of the cross-validation procedure is the number of correct and incorrect predictions parsed by initial conditions and forecast lead times.

Skill assessment

The skill of the prediction scheme is evaluated using the HSS, a measure of the proportion of correct forecasts.13,42 The HSS is calculated as:

$${\mathrm{HSS}} = \frac{{(H - E)}}{{(T - E)}} \times 100,$$ (1)

where H is the number of correct forecasts, T is the total number of forecasts evaluated, and E is the number of correct forecasts expected by chance (T/2 in this two-class scenario). With the multiplication by 100, the two-class HSS ranges from −100 to 100. A set of perfect forecasts garners a HSS of 100, forecasts equivalent to the reference forecast (i.e., climatology) score 0, and forecasts less skillful than the reference forecast receive negative scores. The HSS values may be interpreted in terms of value added relative to a climatological reference forecast. For example, a HSS of 33 indicates twice as many correct forecasts as incorrect forecasts and a skill score of 50 indicates three times as many correct as incorrect forecasts; whereas, the reference forecast with a skill score of 0 indicates an equal number of correct and incorrect forecasts.

Significance of skill

A block bootstrapping approach is used to assess the statistical significance of the HSS values. For every conditional combination (i.e., MJO phase, QBO phase, and forecast lead), we generate a distribution of 1000 skill score values by randomly reassigning the calendar year and shifting the day-of-year indices of the “blocks” of the occurrences of the conditional data. We then perform the verification calculations on the random data in order to construct a distribution of resampled HSS values against which the actual conditional HSS may be compared. In doing so, each block bootstrap sample retains the sample size and potential autocorrelation associated with the conditional data.

Atmospheric river response assessment

For Fig. 3 and the associated discussion, we aim to not only communicate the skill within the prediction scheme, but also characterize the conditional AR response. To achieve this, we evaluate the shift in the probability of above and below normal AR activity for each conditional combination throughout the leave-one-out training and verification process. For example, if a given MJO phase, QBO phase, and forecast lead combination consistently produces a probability of above normal AR activity greater than the probability of below normal activity, we record the condition as resulting in above normal, or increased, AR activity and shade the combination green in Fig. 3.

Model skill comparison

In order to provide a numerical weather prediction skill benchmark, we calculate HSS values for a set of 46-day ECMWF retrospective forecasts from 1995 to 2016. The ECMWF reforecast ensemble prediction system data set consists of 11 members (1 control and 10 perturbed) that are created “on-the-fly.” That is to say that the database is comprised of output from different versions of the ECMWF model, as reforecasts are produced progressively as the the operational model is updated.2 We obtain instantaneous 0000 UTC variables on a 1.5° × 1.5° latitude–longitude grid, from which we calculate IVT and identify AR-like features. We use the first 7 days from every available control run reforecast to calculate calendar-day means from which we calculate the seasonal cycle via fast Fourier transform. Time series of anomalous AR activity are created for each ensemble member’s reforecasts by removing the control run’s seasonal cycle and applying a 5-day running mean. The 50th percentile threshold used to evaluate AR activity is based on anomalies from the control run, subset for the boreal winter. Each member is evaluated separately and for all verification dates within DJFM. Following this approach, 1919 reforecasts are used.

Data availability

MERRA-2 data were obtained from the NASA Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/). ECMWF reforecast ensemble system output was acquired from the World Weather Research Programme/World Climate Research Programme Subseasonal-to-Seasonal Prediction Project database2 (http://apps.ecmwf.int/datasets/data/s2s/). The QBO index was provided by the NOAA NWS CPC (http://www.cpc.ncep.noaa.gov/data/indices/qbo.u50.index). MJO indices were obtained from the Australian Bureau of Meteorology (http://www.bom.gov.au/climate/mjo/) and NOAA Earth System Research Laboratory (https://www.esrl.noaa.gov/psd/mjo/mjoindex/).