In the coastal mountains of western North America, most extreme precipitation is associated with atmospheric rivers (ARs), narrow bands of moisture originating in the tropics. Here we quantify how interannual variability in atmospheric rivers influences snowpack in the western United States in observations and a model. We simulate the historical climate with the Model for Prediction Across Scales (MPAS) with physics from the Community Atmosphere Model, version 5 [CAM5 (MPAS-CAM5)], using prescribed sea surface temperatures. In the global variable-resolution domain, regional refinement (at ~30 km) is applied to our region of interest and upwind over the northeast Pacific. To better characterize internal variability, we conduct simulations with three ensemble members over 30 years of the historical period. In the Cascade Range, with some exceptions, winters with more atmospheric river days are associated with less snowpack. In California’s Sierra Nevada, winters with more ARs are associated with greater snowpack. The slope of the linear regression of observed snow water equivalent (SWE) on reanalysis-based AR count has the same sign as that arrived at using the model, but is statistically significant in observations only for California. In spring, internal variance plays an important role in determining whether atmospheric river days appear to be associated with greater or less snowpack. The cumulative (winter through spring) number of atmospheric river days, on the other hand, has a relationship with spring snowpack, which is consistent across ensemble members. Thus, the impact of atmospheric rivers on winter snowpack has a greater influence on spring snowpack than spring atmospheric rivers in the model for both regions and in California consistently in observations.
Most extreme precipitation on the west coast of North America is associated with atmospheric rivers (ARs; Ralph and Dettinger 2011). When ARs make landfall on the coast of a continent, they can contribute a substantial amount of precipitation in a short period of time (Leung and Qian 2009). In some regions, like California, they can contribute one-third to one-half of total annual precipitation (Ralph and Dettinger 2011). This may also be true farther north as well (Rutz et al. 2014). Focusing on these events makes it possible to investigate how variability in large-scale atmospheric circulation relates to the conditions of extreme precipitation. Snowpack in the western United States is an important seasonal water reservoir, which can be modified by these infrequent transient events. Seasonal snowpack also integrates the influence of all weather in the water year. Both ARs and snowpack have high interannual variability. Using models with multiple ensemble members, more events can be simulated than have been observed over the historical period. This increased sampling of extreme events makes it possible to test whether observed historical relationships could be the result of chance, while also quantifying model biases.
ARs are events in which narrow bands of moisture occasionally extend from the tropics into the midlatitudes. Zhu and Newell (1998), who first quantified “tropospheric rivers” as departures from the zonal-mean meridional moisture transport, showed that most of the total meridional moisture transport occurs in narrow transient bands. Despite the implication that the moisture originates in the tropics, subsequent work has identified ARs as features associated with extratropical cyclones (Bao et al. 2006) rather than a distinct phenomenon. Dacre et al. (2015) composited 200 extratropical cyclones in a cyclone-centric frame of reference to show that ARs ahead of cold fronts are locally generated filaments of high moisture. This places wintertime ARs in the warm sector, which overlaps with the concept of the warm conveyor belt in the context of an extratropical cyclone (Ralph et al. 2004) and generally involves warmer temperatures as well as enhanced moisture (suggesting enhanced sensible as well as latent heat transport). Regardless of the synoptic processes involved in their formation, identifying ARs provides an event-based method for studying both the large-scale conditions that lead to extreme precipitation and the regional impacts.
Neiman et al. (2008), the first to consider the influence of ARs on snow, studied California and another region farther north that includes Oregon, Washington, and part of British Columbia from 1997 to 2005. Using observations and reanalysis, they noted that lower-elevation sites in the Cascades can receive less snowfall than normal in association with ARs, and those sites sometimes lose snowpack. Nonetheless, Neiman et al. (2008) concluded that ARs do not significantly modify snowpack in their northern domain. They found California, however, has enhanced snowpack with more ARs in winter and reduced snowpack with more ARs in spring.
Using the AR events identified in Neiman et al. (2008) for California, Guan et al. (2010) determined that precipitation days during ARs add 4 times as much snow water equivalent (SWE) compared with days during storms not associated with ARs. They found that the ARs with a higher impact on SWE occur when there is colder surface air temperature than is present over the Sierra Nevada during lower-impact ARs. In Guan et al. (2013), they associated the colder temperatures of higher-impact ARs with the higher frequency of ARs in years with negative indices for the Arctic Oscillation (AO) and the Pacific–North American (PNA) pattern. Finally, Guan et al. (2016) considered the contribution of ARs to rain-on-snow events in the Sierra Nevada, finding that while 50% of rain on snow is due to ARs, only 15% of ARs cause rain on snow. It occurs most frequently in association with ARs in March, when temperatures have begun to warm, but snow still remains on the ground. We would expect rain on snow to influence snowpack more in spring than winter in the northwestern United States as well. Huning et al. (2017) found that the AR-detection criteria can affect the percent of seasonal snow accumulation attributed to ARs in the Sierra Nevada because of differences in the temperature on the days that are detected by one method versus another.
Most work on the influence of ARs on snowpack for western North America has focused on California (Guan et al. 2010, 2013; Kim et al. 2013; Guan et al. 2016; Huning et al. 2017). This study aims to further investigate the relationship between North Pacific ARs and mountain snowpack, with a focus on the U.S. Northwest—defined as Washington, Oregon, Idaho, and western Montana. It is worth including the eastern part of the Northwest region because according to Rutz et al. (2014) as much as a quarter of cool-season precipitation is associated with ARs in the inland Northwest, which is also quite mountainous. Thus, we expect to see some influence of ARs throughout the Northwest, although the precipitation enhancement during an AR event will be lower in the inland parts of the region than in the coastal ranges (Neiman et al. 2008). We will separately consider California for comparison with the Northwest and to check for consistency with previous work.
While much of the work on ARs uses a mixture of observational data and reanalysis, a number of studies have also considered how ARs are simulated in models. Some have used regional models to study a particular location like California (e.g., Kim et al. 2013), the broader western United States (e.g., Leung and Qian 2009), Europe (Lavers et al. 2011; Lavers and Villarini 2013), or South America (Viale and Nuñez 2011). Comparing the historical period in reanalysis to 24 models from phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012), Gao et al. (2015) found that the number of AR days in the models and reanalysis are generally in reasonable agreement, as does Warner et al. (2015) with 10 CMIP5 models.
To address the interaction of large-scale atmospheric features like ARs with small-scale features like mountain snowpack requires a climate model with both the spatial coverage of a large domain and the resolution to capture sufficient mountain topography. Hagos et al. (2013) compare AR-like features in aquaplanet simulations produced by the Model for Prediction Across Scales (MPAS)–Atmosphere (MPAS-A; Skamarock et al. 2012; an earlier version from that used in this study) at multiple resolutions and find significant sensitivity of AR frequency to model resolution, partly a result of the sensitivity of the simulated jet stream to model resolution. The same MPAS-A model with realistic land surface and historical prescribed sea surface temperature (SST) is evaluated in detail in Sakaguchi et al. (2015) by comparing global simulations at quasi-uniform low and high resolutions with global variable-resolution simulations with refined resolution over the Americas. The simulations in the refined domain at 30-km grid spacing reproduce the climatic features of the global quasi-uniform high-resolution simulation at 30-km grid spacing, but at a fraction of the computational cost.
Thus, global variable-resolution modeling may be a computationally efficient approach for high-resolution modeling in regions such as the western United States with complex topography where model grid spacing has important effects on the simulation quality. Computational efficiency is important because to study low-frequency extreme events, it is necessary to simulate many years. Using a multiscale model with higher resolution over the area of interest, we can simulate the interaction of extreme precipitation and snowpack in mountainous areas in a computationally feasible fashion. With a global model, forced by prescribed historical SST, we simulate a sufficient number of years to study variability. A small ensemble of simulations with perturbed initial conditions is used to test the extent to which SSTs force variability in ARs. This approach allows us to explore relationships between ARs and snowpack in our region of interest, identify features that are consistent with the historical observations, and interpret relationships observed over the short historical period in the context of internal variance. We will identify whether atmospheric rivers tend to enhance or deplete snowpack depending on location and season, as well as the extent to which these relationships found in a small model ensemble can be detected in nature above the internal variance.
In this study, we use a nonhydrostatic version of MPAS-A (Skamarock et al. 2012) combined with Community Atmosphere Model, version 5 (CAM5) physics [(MPAS-CAM5) or “the model”], as described in Zhao et al. (2016). The specific version of MPAS-A used is version 2, and the version of CAM5 is 5.3.18. The nonuniform mesh has a region of ~30-km grid spacing centered at 40°N, 125°W. The grid spacing of the rest of globe is ~120 km. The higher-resolution region spans most of the mountainous regions of western North America, as well as the eastern North Pacific (see Fig. 1).
SST and sea ice are prescribed on the basis of historical observations using the merged Hadley–OISST and sea ice concentration datasets (Hurrell et al. 2008). The land surface model is the Community Land Model, version 4.5.50 (CLM4), which includes the Snow, Ice, and Aerosol Radiation (SNICAR) model (Flanner et al. 2007, 2009), and other parameterizations discussed in (Oleson et al. 2013). Three ensemble members are produced that are identical except for the perturbed initial conditions. Initial conditions were varied by starting the simulations on different days in early January 1981 or 1982. The simulations all span 1982–2012. For analyses, we drop the period prior to December 1982. All output is remapped from the MPAS mesh to a 0.25° latitude–longitude grid before analysis. Then, AR days are detected using a method to be described below.
Standalone CLM4 simulates snow adequately compared with other products (Toure et al. 2016). Snow accumulates and ages in up to five layers, and melts, sublimates, and interacts with the vegetation to influence albedo and near-surface stability. Simulation of mountain snowpack is sensitive to model resolution and physics parameterizations (e.g., Leung and Qian 2003). Using the same CAM5 physics as used in this study, Rhoades et al. (2016) found that the simulated snowpack in the western United States improves with model resolution increasing from 55 to 7 km and with the use of the Morrison–Gettleman microphysics version 2 relative to version 1 that is used in CAM5. The MPAS-CAM5 modeling framework has the advantage of offering greater resolution in our region of interest, employing the abovementioned versions of the same atmospheric physics and land surface model of a recent generation of the Community Earth System Model (CESM; Hurrell et al. 2013). While we have not tuned the parameterizations for the higher-resolution region, our focus is on mountainous coastal locations where large-scale precipitation dominates, so we expect that tuning would have a minimal influence.
We compare the model with the 0.9° × 1.3° CAM5.1 simulation forced with historical SSTs from the CESM contribution (Hurrell et al. 2013) to CMIP5, to illustrate the effect of the higher-resolution region of the MPAS-CAM5 grid. We refer to this simulation hereafter as the 1° CESM AMIP and use climatologies for the period from December 1982 through the end of 2005, when the simulation ends.
We also examine observations and reanalysis to evaluate the model. For precipitation, we use the North American Land Data Assimilation System (NLDAS-2) precipitation dataset. At 0.125° resolution, the NLDAS-2 product assimilates gauge data into an analysis of precipitation and adjusts for elevation using the Parameter-Elevation Regressions on Independent Slopes Model (PRISM; Daly et al. 1994) climatology.
The gridded Snow Data Assimilation System (SNODAS; National Operational Hydrologic Remote Sensing Center 2004) product assimilates the Snowpack Telemetry (SNOTEL) station data and satellite-based SWE estimates into a model of the land hydrology, which is important to build a self-consistent data product even where stations are not available, for example, at higher elevations. The SNOTEL stations use snow pillows that derive SWE from the weight of the snowpack (Serreze et al. 1999). Because SNOTEL stations are sparsely located in remote mountainous areas at relatively high elevation, it is useful to complement the SNOTEL data using gridded data product such as SNODAS. The model used to produce SNODAS has three snow layers and two soil layers and the energy/mass balance is not fully coupled with the atmosphere (Barrett 2003). In the SNODAS product, large areas far from stations have substantial uncertainties: satellite-derived estimates for snow-covered areas at higher spatial scale are not independently calibrated for every region (Salomonson and Appel 2006) and can be quite uncertain in the presence of vegetation (Barrett 2003). Additionally, topography, slope aspect, and vegetation can affect the snowpack on the scale of catchments in ways that are not properly captured by the model used to assimilate data and create the product. It is useful to compare the MPAS-CAM5 setup with SNODAS for realism as an observationally informed estimate of the state of the snowpack. We use both SNODAS and SNOTEL to evaluate the model snowpack climatology for SWE.
The higher-resolution observational data (NLDAS-2 precipitation, SNODAS SWE) and lower-resolution CESM AMIP 1° data are regridded with bilinear interpolation to the same 0.25° grid as the MPAS data before calculating root-mean-square error (RMSE) or pattern correlations. For a higher-resolution elevation profile to correspond with the SNODAS data, we query the 30-m U.S. Geological Survey National Elevation Dataset (Sugarbaker and Carswell 2011).
For further snowpack analysis, we use the SNOTEL station SWE for the Northwest and California. We use a period beginning in 1982 to cover a longer period than the SNODAS data. For both SNOTEL and SNODAS, we use data through 2016. Although this period does not entirely overlap with the MPAS-CAM5 simulated period, it provides the greatest possible number of seasons.
We also use the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis [ERA-Interim (ERA-I); Dee et al. 2011] for large-scale fields: 2-m air temperature, 500-hPa heights, and total column-integrated precipitable water. We construct daily means from the mean of 6-hourly instantaneous output. Neiman et al. (2008) showed that AR identification using the NCEP–NCAR reanalysis agreed well with an AR catalog for the eastern North Pacific over 1997–2005 based on satellite measurements of integrated water vapor from the Special Sensor Microwave Imager (SSM/I). Rutz et al. (2014) were the first to use ERA-I to detect ARs.
With daily data from either MPAS-CAM5 or ERA-I, we identify ARs using our algorithm described here. For each day, a two-dimensional snapshot of column-integrated precipitable water is examined in the domain from 170° to 110°W longitude and 25° to 65°N latitude. Distinct features exceeding a 2-cm threshold are identified. Features are considered to be landfalling in a region if some portion of the feature overlaps a grid cell in the appropriate latitude range with land fraction greater than 0.2. The latitude range for the Northwest is 41.5°–49.5°N and for California is 32.5°–41.5°N.
If a labeled feature intersects the specified section of coast, then it is tested for several other criteria. The feature must exceed 2000 km in length (between its northeast-most and southwest-most grid cell) and the total area divided by the length (as a proxy for width) must be less than 1000 km. These dimensions, together with the 2-cm column-integrated precipitable water threshold, were determined to be empirically appropriate for AR identification in the eastern North Pacific by Ralph et al. (2004). The landfall extent must also be less than 1000 km, where the landfall extent is the distance between the highest and lowest latitude of landfall of the feature. This method is applied to all time slices in a DataArray data structure in the Python library xarray (Hoyer and Hamman 2017), which makes the operation fast and efficient (algorithm code is available at http://github.com/theNaomiG).
In this study, we identify AR days using daily data, although the algorithm can be applied to data at other time frequencies because it looks at each time independently. Consecutive AR days are considered to be a part of the same AR event in some composites of days before and after AR events. Thus, the count of AR days for a month can exceed the total number of events, especially when stronger ARs occur. Primarily we focus on total AR days, creating the AR indices (described below) as a measure of frequency of days with AR conditions.
We use daily data here because it is sufficient to study interannual variability in AR occurrence. Another effect of choosing daily data is that we effectively exclude all short-duration events but the most intense from consideration; an AR that lasts only part of a day would have to be more intense to cause the 24-h precipitable water to exceed the threshold. The highest hourly precipitation rates are associated with the longer-duration ARs (Ralph et al. 2013), so the subdaily ARs that are excluded are less likely to impact monthly or seasonal means of hydrologic fields. Additionally, we note that an AR that spanned close to 24 h but was split between two days might be excluded, whereas a similar event contained entirely in one model output day would exceed the threshold. This issue should minimally bias the overall distributions of AR days if such an event is just as likely to occur in any month of any year.
The dates on which landfalling ARs are identified are tabulated for use in composites. The total number of landfalling AR days per month is also counted, where i indicates ensemble member (1, 2, 3) or ERA-I (4) and m is the month (November–May). Each month’s count is standardized into a unitless AR index by subtracting the mean and dividing by the standard deviation for that month:
We also construct an index for winter (December–February) and spring (March–May) seasons as above but with m denoting season instead of each month. We use this process to generate AR indices for each region separately to study the influence of ARs on snow in the mountains of the Northwest and California.
To analyze ARs across all ensemble members at once, we start with the counts for a given month from each ensemble member and concatenate them into one list with 3 times the length:
before standardizing as in Eq. (1). Similarly concatenated meteorological and surface fields can be regressed against this longer AR index.
Indices are also defined for the cumulative number of ARs for December through the month in question of the same water year. In those cases, the time series of the number of ARs in Eqs. (1) and (2) is replaced by the total across months:
We employ these indices to explore how the interannual variability in snowpack and large-scale circulation relates to the interannual variability in ARs. In that portion of the analysis, we fit a linear regression between a monthly or seasonal-mean 2D field and an AR index for one of the regions and focus on the sign and magnitude of the slope of that relationship at each point. We do this for each ensemble member. If all three are consistent with one another and the observations, this suggests that the relationship may be robust, although we cannot prove that it would always be so if we had more ensemble members. Where the three ensemble members do not agree, internal variance likely plays the dominant role, and observational relationships are interpreted with greater caution.
We first evaluate how the model performs at simulating AR counts, precipitation, and snowpack. We examine the conditions that are associated with AR days and the large-scale features of seasons with more ARs, comparing model and observations. Then, we use the small ensemble of model simulations to examine the potential for AR predictability and finally explore the relationship between AR frequency and snowpack.
3. Model simulation of key fields
We begin by evaluating how the model performs at simulating precipitation and snowpack. Then, we test the model against reanalysis using the AR algorithm and compare large-scale conditions between model and reanalysis on AR days and across years with varying AR counts.
a. Precipitation and snowpack
Given the resolved elevation shown in Fig. 1, we expect MPAS-CAM5 to better capture orographic precipitation in the Olympic Mountains, Cascades, and Rocky Mountains and some amount of the rain shadow to the east of the Cascades than the typical global model used for large-scale climate studies. Figure 2 shows the model precipitation climatologies for December–May compared with the NLDAS-2 gridded dataset. The RMSE for the whole area (32.5°–49°N, 111°–125°W over land) is 1.1 mm day−1. MPAS-CAM5 is on average biased high with respect to NLDAS-2 (see Fig. 3a for an example over the Cascades). In Canada, the NLDAS-2 dataset does not capture the full extent of precipitation; note the discontinuity at the border that is not present in the model. The discontinuity in NLDAS-2 is related to the use of two precipitation datasets, one for the United States and one for Canada, without homogenization of the datasets to ensure spatial continuity. The CESM AMIP simulation at 1° resolution does not capture the major features of the spatial distribution of precipitation in the western United States (Fig. 2c), missing the Sierra Nevada of California and the Willamette valley in Oregon, for example, and is slightly more biased than MPAS-CAM5 compared to NLDAS-2 on the whole, with RMSE of 1.2 mm day−1.
The influence of the higher spatial resolution is even more notable for snowpack (see Figs. 2d–f, 3b) than for precipitation. CESM AMIP only simulates substantial snow cover where the model’s surface geopotential height is high enough to simulate the cold temperatures required for a seasonal snowpack, with RMSE compared to SNODAS of 11.2 cm for the December–May climatology. It completely misses features such as the Oregon and California mountain snowpacks. The MPAS-CAM5 (Fig. 2e) captures the major snowpack in mountain areas of the West, with an RMSE of 8.24 cm compared to SNODAS, although, in some locations, one can see from Fig. 2e that MPAS-CAM5 does not resolve distinct peaks in snowpack where topography is not resolved. In Fig. 3b, we show a longitudinal mean across a slice of the Cascades from west to east. MPAS-CAM5 is biased high compared with SNODAS on the windward side and at lower elevations on the lee side. For the December–May mean, the pattern correlation for MPAS-CAM5 snowpack with SNODAS is 0.76. For CESM AMIP, the pattern correlation with SNODAS is only 0.48. Although SNODAS itself can be biased, especially where fewer stations are present, the greater realism gained from the 30-km resolution is clear, as overall MPAS-CAM5 simulates large-scale snowpack patterns far better than 1° CESM AMIP.
b. AR counts
In Figs. 4a and 4b, we show the number of ARs per season impacting the Northwest coast for the whole 30-yr time span in the three MPAS-CAM5 ensemble members and ERA-I. Counts of ARs impacting California are shown in an equivalent figure in the online supplemental material (Fig. S1). The model simulation of AR days is biased high compared with ERA-I in most years and in the mean. Variability from year to year is also greater in the model compared with ERA-I, especially in winter. Finally, we note a higher number of ARs occurs in winter than spring across all ensemble members and ERA-I. The monthly climatology of the number of AR days (Fig. 4c) further reveals that the model peaks in January and ERA-I peaks in November. It is unlikely that a longer historical record would change this discrepancy in the month of peak AR number, as the bias in AR climatology exceeds one standard deviation of the model ensemble in almost every month.
Despite the biases, the model considered here performs better than typical global climate models and previous versions of the MPAS model. Hagos et al. (2015) found that AR frequency was underestimated in AMIP-style simulations with MPAS-CAM4 (an older version of MPAS that is hydrostatic and uses CAM4 physics). AR frequency depends on the jet position and strength, which are biased equatorward and weaker in the hydrostatic version of the model in Hagos et al. (2015) and Sakaguchi et al. (2015). Zhao et al. (2016) anticipated that nonhydrostatic MPAS would also perform better because of its improved simulation of low-level westerlies.
The bias in the number of AR days depends on the detection method; for example, a threshold chosen relative to the climatological mean state instead of an absolute threshold can remove such a bias (e.g., Hagos et al. 2015). For our purposes, a constant moisture threshold is the simplest way to be consistent across datasets. By normalizing our AR indices, biases in absolute AR count do not influence our results.
To test how the AR indices relate to extreme precipitation, we correlate AR indices with the time series of the number of winter days exceeding the 95th percentile in winter-mean daily precipitation for all points (for the Northwest, see Fig. 5; for California, see Fig. S2). While the AR index is correlated with this metric of extreme precipitation frequency, the correlation is, unsurprisingly, not equal at every point in the region. The general southwesterly flow associated with ARs may explain why higher correlations are seen in the northern part of the region, especially farther inland from the coast. An AR that impacted the Oregon coast would need to approach with a more zonal orientation for the trajectory of elevated winds and moisture to be directed toward southeast Oregon and in that case would encounter the rain shadow of the Cascades.
c. AR evolution
The elevation and climatological temperature of a mountainous area in a given season influence the AR–snowpack relationship and so should the large-scale atmospheric circulation associated with AR days. A landfalling AR is generally associated with a warm 2-m temperature anomaly over land downstream (Fig. 6). There is a striking warm anomaly in winter, which peaks the day after the event, and persists at least three days later in both model and reanalysis. A similar, although weaker, anomaly is present in spring (Fig. S3).
In both seasons, a trough occurs just offshore and a ridge onshore, setting up southwest flow and channeling a band of high moisture into the warm sector, consistent with expectation. The magnitude of warm anomaly varies with ensemble member. Whether this influences the AR–snowpack relationship will be explored in section 4b.
The 2-m temperature composites for AR days in Fig. 6 are consistent with the composites in Warner et al. (2012) of 850-hPa temperature on days with the most extreme precipitation along the Northwest coast. Kim et al. (2013) found good agreement in 700-hPa temperature (and geopotential height) composites between AR days in ERA-I and WRF, with the latter at comparable resolution to the region in MPAS-CAM5. The magnitude of the warm anomaly associated with ARs may be a function of the AR-detection method: As discussed in Huning et al. (2017), defining ARs from total column precipitable water thresholds (as we have done) yields fewer but generally warmer AR days than a definition based on vapor flux.
Additionally, we examine how the freezing level evolves over the course of AR events in the Northwest and how this compares with the climatological freezing level (Fig. S4). The freezing level several days before an AR event is generally higher than climatology in winter but lower in spring. In both seasons, when an AR makes landfall, the freezing level moves up in elevation. In winter, the freezing level recedes entirely from the central Cascades but recovers to its former position three days later. In spring, the freezing level recovers only to approximately the climatological level. The patterns are similar across ensemble members.
d. Potential for AR predictability
The times series of indices for ARs reaching the Northwest are not significantly correlated across ensemble members (see the tables in the supplement), suggesting that SSTs are not a strong control on AR frequency in the Northwest, at least in this type of model configuration. This is consistent with Dong et al. (2018), who noted that SSTs only explain about 20% of year-to-year variability of extreme precipitation (mainly associated with ARs) in the western United States, based on analysis of multimodel simulations and large ensemble simulations. Guan and Waliser (2015) found a significant anomaly in AR frequency associated with positive Niño-3.4 index in the historical period from 1997 to 2014 for the northwestern United States and British Columbia but not California. However, we find no evidence that ENSO-related SST forcing is offering any AR predictability in the Northwest, otherwise we would expect the AR indices to be significantly correlated with one another across ensemble members. Correlations maps of Pacific SST with the Northwest AR indices vary across the MPAS-CAM5 ensemble members and ERA-I (Fig. 7). ARs reaching California, on the other hand, tend to be associated with an El Niño–like pattern of tropical Pacific SST (see Fig. S9). The question of predictability of AR frequency is important for water resource planners, especially in California, where a greater portion of the annual precipitation is associated with ARs, but we caution against drawing conclusions on the basis of short historical time series.
4. ARs and snowpack
Here we investigate how AR events modify snowpack and influence interannual snowpack variability in the Northwest and California.
a. Composites of snowpack centered on AR days
To understand the process by which ARs influence snowpack in winter, we examine the day-to-day variation in snowpack before and after AR events, defined as consecutive sequences of AR days. Most events (73%–84%) last a single day. We composite the snowpack change (the day-after minus day-before anomalies) for AR events (Fig. 8). In the SNODAS composites, the snowpack SWE increases with AR events over most of the mountainous regions of the Northwest and in the eastern plateau. The pattern is similar in the model but with generally smaller increases and with a decrease in much of the south Cascades, Olympic Mountains, and large volcanoes in Oregon. MPAS-CAM5 has insufficient horizontal resolution to capture the relevant small-scale relief of the highest peaks. SWE changes in the next several days after an AR event have similar patterns in both SNODAS and the model. In the SNODAS version, the SWE increase can be greater, with a maximum increase from one to three days after the event of 4 cm instead of 2 cm for MPAS-CAM5. The SWE increase occurs because snow is still falling when an AR event lasts more than one day or is followed by precipitation from postfrontal clouds even if AR conditions are no longer present.
b. Interannual effects
Beyond the immediate impact of an AR on snowpack, we wish to examine the cumulative impact of high-AR years on seasonal snowpack. To detect the spatial patterns of snowpack response to AR frequency, we regress monthly mean SWE on the AR indices for each of the months. Although AR frequency is climatologically high in November, there is little relationship with snowpack (see Fig. 9a) because there is too little accumulated snowpack (not shown). The lack of a relationship in November justifies focusing on the subsequent months in the analysis that follows.
A small but clear relationship is seen in December and January in Fig. 9. January, which is the coldest month, has the most widespread positive relationship between snowpack and AR count. ARs in December are more widely associated with lower monthly mean snowpack. The locations in February of positive linear regression slopes are confined to the (higher and colder) inland mountains, but February has negligible locations of significant linear regression slopes. In all winter months, there is an area of positive linear regression slope for snowpack in the highest part of the northeast Cascades. This enhanced snowpack is collocated with the coldest parts of the Cascades, but the regression is not significant there (see Figs. 9f–j). For the Cascades, the only month with negative linear regression slopes is December. The Rocky Mountain location’s positive coefficients are only significant in January. The lack of significant regressions in February could be due to its greater year-to-year variance in snowpack (Fig. S10), as the effect of all events that lead to accumulation, melt, or ablation of snowpack is cumulative up to the month in question.
Next, we investigate the influence of internal variability on the relationship between AR frequency and snowpack further by comparing seasonal regressions from the three ensemble members. The pattern of snowpack regressions on ARs in winter (Figs. 10a–c) is similar among ensemble members, although the negative relationship of snowpack and AR count is stronger in ensemble member 3. We note that ensemble member 3 also exhibits warmer anomalies in association with AR days in winter (Figs. 6k–o). The difference between the maps of daily changes and the seasonal-mean regressions of snowpack on AR indices (e.g., cf. Fig. 8 to Figs. 9e, 11a) could be due to the fact that more ARs occur in winters with overall warm temperature anomalies, where the warm anomaly is not limited to the AR days, as shown in the supplement (Figs. S5–S8), so more snowpack could melt on non-AR days, contributing to the seasonal effect.
In spring, there is no consistent synchronous relationship between Northwest snowpack and AR frequency in the model. The regression patterns in spring lack a common sign across ensemble members in much of the Cascades and in some parts of the Rocky Mountains. Ensemble member 3 does not stand out as the warmest in spring (Fig. S3k–o), but neither is it notably colder, which would be consistent with Fig. 10f. Thus, variations across ensemble members in 2-m air temperature composites do not explain the differences among ensemble members for snowpack regressions in spring, and another explanation is needed.
Next, we compare these relationships that occur in the model with observational estimates of recent years, using the ERA-I AR index combined with SNODAS and SNOTEL station data (Fig. 11). In winter, linear regression slopes are consistent between the model and SNOTEL in most of the locations for which we find significance in the model results (Fig. 9j). Notably, the linear regression slopes in the north Cascades have different signs between the two periods, and the signs also differ between ensemble members, suggesting the lack of a robust relationship between ARs and snowpack there. The northeast Cascades stations with positive linear regression slopes over the longer period (Fig. 11b) are outside the area in which the regressions are significant in the model (Fig. 9j). Even where the model ensemble members and observations over both periods are consistent, most Northwest stations have small linear regression slopes and the observational period is too short to be significant at P < 0.05 (Fig. 11b).
In the spring, when the MPAS-CAM5 ensemble members’ linear regression slopes have no consistent sign in the Northwest, in observations the sign is mostly negative (Figs. 11c,d) and, as in winter, not significant. As a caution, we note that stronger negative linear regression slopes (which are significant at P < 0.05) are found using only the period from 2003 onward available for SNODAS (Fig. 11c). SNOTEL is consistent with SNODAS over the same period. However, during the shorter period, April and May both have only three years in which a single AR impacts the Northwest and no years with more than one—not enough information upon which to base a robust index of interannual variance. Our analysis that follows suggests that spring snowpack is not influenced by spring AR frequency in a systemic fashion.
Snowpack reflects the cumulative influence of weather since the start of the water year, so by spring the synchronous relationship between ARs and snowpack is diluted. Additionally, the spring signal could vary on whether precipitation falls as rain or snow, whereas in winter it more consistently falls as snow. For the same reasons, the linear regression slopes in Fig. 9 in February are less significant than December. Indeed, unsurprisingly, there is no significance in the individual months of March–May in plots (not shown) similar to Fig. 9. In section 3c and in the supplement (Figs. S5–S8), we show that the model represents spring processes as faithfully as winter and that spring large-scale conditions are more influenced by internal variance. Thus, any synchronous relationship between ARs and monthly to seasonal snowpack observed in the real world in spring should be treated as the potential result of chance. In section 4d, we show that cumulative AR-indices have a clearer effect on the spring season.
c. California comparison
Next, we compare the relationship between ARs and snowpack in the Northwest to that in California, shown in Figs. 12 and 13. As in the Northwest, we see a consistent relationship between ARs and snowpack in California for the winter season in both the model and observations. The propensity for increased AR frequency to result in greater winter snowpack in the Sierra Nevada has been known since Neiman et al. (2008). It is reassuring that we reach the same conclusion using different methods. The positive linear regression slopes in the Sierra Nevada contrast with the negative regressions in much of the Cascades because the Sierra Nevada are higher and portions are climatologically colder than the Cascades. Additionally, the winter linear regression slopes in the Sierra Nevada are significant (P < 0.05) in the observations in addition to the model, and the result is more spatially coherent.
In the regressions of spring snowpack on synchronous spring AR index, we see that in the model (Figs. 12d–f), the sign of the relationship between AR frequency and Sierra Nevada snowpack in spring depends on the ensemble member, just as it did for the mountains in the Northwest. For SNODAS and SNOTEL data regressed on the ERA-I AR index, a higher frequency of spring ARs is also associated with more snow overall in the Sierra Nevada for the longer period (Fig. 13d) but less at most stations in the shorter, more recent period (Fig. 13c). Almost none of the stations have a significant linear regression slope, suggesting that little can be learned from synchronous spring relationships. Two ensemble members reproduce the observed pattern in the longer period in SNOTEL, while one does not. This is consistent with the sign of the linear regression slopes in spring being influenced by internal variance and an example of a situation where even a small ensemble can suggest caveats to conclusions based solely on observations over a relatively short historical period.
d. Cumulative ARs and snowpack
Finally, we test whether spring snowpack is influenced by the preceding winter ARs. Instead of testing the relationship between spring ARs and spring snowpack, we start with the cumulative AR count, beginning in the preceding winter up to the month in question [Eq. (3)]. We fit a linear relationship between a cumulative AR index and spring SWE at each grid cell, just as we did above for the synchronous AR indices. We obtain a regression that is significant over a greater geographic extent in the mountains for all months extending into the spring in the model (regression slopes shown in Fig. 14 and Fig. S11).
The month in which the magnitude of the linear regression slopes maximizes is March, which is also the peak month for monthly area-averaged SWE climatology. The sign of the linear regression slopes in all spring months tends to be the same as for those found earlier for the winter months. That we must account for ARs in the preceding winter in order to get a significant snowpack regression in the spring suggests that antecedent winter snowpack conditions, themselves influenced by the number of AR days, are more important to spring snowpack than the additional snow that accumulates or melts in association with spring ARs.
The cumulative AR-index-based linear regression slope pattern in the model is consistent with observations for the Northwest (Fig. 16), but most stations are still not significant at P < 0.05. In observations for California, we note that the sign of the coefficients for the spring months’ snowpack regressed on cumulative AR index is significant (Fig. S13), which is not true of the relationship based on synchronous spring ARs alone in the observations.
We have employed an experimental version of a state-of-the-art multiscale model (MPAS-CAM5) to explore relationships between atmospheric rivers and seasonal snowpack in the western United States. The model is uniquely suited to this problem because it allowed us to consider large-scale features and the resultant effects on snowpack in a seamless and self-consistent framework. We were able to resolve mountain snowpack with sufficient resolution to note important spatial variations across latitude and between windward and leeward sides of a range. We chose to simulate 30 years with three ensemble members to do a thorough study of variability, rather than simulate the higher spatial resolution necessary to better characterize differences from one watershed to another. At the scales we consider, the model climatology is in accord with historical datasets.
We found no significant correlation between the Northwest AR indices of the three MPAS-CAM5 ensemble members, indicating that SST (and hence ENSO as well) has little control on Northwest AR frequency. The inconsistent response of AR frequency to ENSO is not promising for improved AR seasonal prediction in the Northwest, but there may be more coherent correlations with tropical SST for the ARs impacting California. Further work is needed to establish the mechanisms by which large-scale variability influences AR frequency for distinct parts of the west coast of North America.
The structure of ARs and the environment that they encounter determine whether ARs produce rain or snow and their impact on the overall snowpack. Focusing on ARs landfalling in the Northwest, the model simulates the day-to-day changes in snowpack with AR events similarly to SNODAS, with the exception of some high-elevation peaks that are not captured in MPAS-CAM5. In both model and observations, the day-after response differs from the seasonal snowpack regression on AR index, likely because the seasonal large-scale conditions in high-AR years are anomalously warm on the whole over land.
In the Northwest, in both model and observations, years with more ARs are associated with less seasonal-mean snowpack in most of the Cascades in winter, with the exception of parts of the north Cascades where the signal is not robust. Compared to observational estimates, the model produces a snowier winter response in the Rockies to ARs making landfall along the Northwest coast; it can be seen in Fig. 3 for one part of the Cascades that the model has too much precipitation penetrating inland. In California’s Sierra Nevada in winter, a significant positive relationship between AR count and snowpack is present in both model and observations.
In contrast to winter, we found no significant synchronous relationship between snowpack and AR frequency in spring in the Northwest or California. We cannot rule out that the observed association between more frequent spring ARs and less spring snowpack is the result of chance over such a short period of observations. In fact, it seems likely that no significant relationship could be found for the spring season because the small effect of AR frequency on snowpack is overwhelmed by the cumulative effects of all sources of snowpack interannual variability. Further, a cumulative AR index better predicts spring snowpack in a given month or season than an AR index based on the month or season in question alone.
Finally, we consider the implications of our findings on variability in snowpack associated with AR frequency to potential future climates. In winter, the season for which a robust relationship currently exists between ARs and snowpack, we see that the biggest constraint on the sign of the snowpack response is geography. In high mountain locations, such as the Rockies and Sierra Nevada, that are climatologically colder and drier than the Cascades because of the combination of elevation, latitude, and distance from the coast, the result of anomalously high AR frequency is to enhance the winter snowpack.
In the Cascades, which are lower and closer to the coast in the path of the storm track—and thus climatologically warmer and wetter—the relationship is reversed, but the effect is small and significant in the model only when all three ensemble members are considered together. We considered the ensemble members both in aggregate and in comparison with one another to learn in which locations effects are robust and in which locations internal variance dominates and no consistent effect exists. High AR frequency in the southern Cascades is associated with a decrease in seasonal snowpack, but the north Cascades are more influenced by internal variance. An increase in ARs is projected with warming (Warner et al. 2015; Gao et al. 2015; Hagos et al. 2016), and loss of snowpack from warming in the Cascades would only be amplified by the increased ARs given a negative linear regression slope. As the climate warms, we also reason that, because of thermodynamic considerations, areas where greater AR frequency is associated with higher seasonal snowpack should shrink.
Changes related to atmospheric circulation, which are not consistently projected across models (Gao et al. 2015), would primarily affect the frequency of events. Any systematic shifts in the Hadley circulation or storm tracks would influence the latitude to which moist subtropical air extends and could change the relative number of ARs at any given latitude. The net impact of ARs on snowpack could be lower if there are systematically fewer events or higher if there are systematically more events in a given location.
The importance of internal variability that we have demonstrated for the spring has implications for the interpretation of future changes, as a different sign in the regressions could develop even in the absence of changes in forcing. The cumulative influence of winter and spring ARs on spring snowpack may remain more important than the synchronous spring signal in California, where this is clearly the current case, motivating a better understanding of future winter changes. The AR–snowpack relationship should be secondary to the warming trend itself when the warming signal in the western United States region becomes significant.
This study is supported by the U.S. Department of Energy Office of Science Biological and Environmental Research as part of the Regional and Global Climate Modeling program. This research used computational resources from the National Energy Research Scientific Computing Center (NERSC), a DOE User Facility supported by the Office of Science under contract DE-AC02-05CH11231. PNNL is operated for DOE by Battelle Memorial Institute under contract DE-AC05-76RL01830. All raw model simulation data were archived at the NERSC High Performance Storage System (HPSS).
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0268.s1.