1 Introduction

Accurate seasonal precipitation prediction can help support water resource management; there has been an increasing call for better predictions, particularly over regions like the western United States (WUS) where a multi-year persistent drought recently adversely affected regional water supply, agriculture, ecology, and wildfire risks (Seager et al. 2015). The El Niño-Southern Oscillation (ENSO) phenomenon is the major driver for interannual variations of precipitation (Ropelewski and Halpert 1986, 1989, 1996; Kahya and Dracup 1993), as well as the major source of seasonal predictive skill over the WUS (Jia et al. 2015; Yang and DelSole 2012). During the 2015/16 winter season, a very strong El Niño event developed over the tropical Pacific, referred to as a “Godzilla El Niño” in the popular press (Kintisch 2016). Its intensity—based on NIÑO3.4 (170°W–120°W, 5°S–5°N) sea surface temperature (SST) anomalies—was comparable to the extreme event in 1997/98 (Fig. 1). However, the observed precipitation anomalies over the WUS in Winter 2015/16 differed markedly from the canonicalFootnote 1 El Niño pattern (Fig. 2a, b), with southern California notably maintaining a dry (instead of an expected wet) anomaly. The operational seasonal climate prediction models from the North American Multi-Model Ensemble (NMME) (Kirtman et al. 2014) successfully predicted the tropical Pacific SST aspects of this extreme El Niño event months in advance. However, the models predicted a canonical El Niño precipitation pattern with wet conditions over Southern California, thus failing to predict the observed precipitation anomalies over WUS (Wanders et al. 2017). In contrast, the observed precipitation anomalies during the winter of the other most recent major El Niño in 1997/98 largely resembled the typical pattern over the WUS (Fig. 2a, f), and gave anecdotal support for a repeat in the original 2015/16 prediction.

Fig. 1
figure 1

The observed DJF SST anomalies for a 1997/98, b 2015/16 and c 2015/16 minus 1997/98. NINO3.4 index (SST anomalies averaged over the domain 5°S–5°N and 170°W–120°W, blue box in a and b (blue bars) and NINO3 index [SST anomalies averaged over the domain 5°S–5°N and 150°W–90°W, black box in a and b] (red bars) for 1997/98 and 2015/16 are shown in d

Fig. 2
figure 2

The winter (DJF) precipitation anomalies for a the canonical El Niño composite, b 2015/16 observations, c 2015/16 P1 hindcasts, d 2015/16 P2 hindcasts, f 1997/98 observations, g 1997/98 P1 hindcasts, h 1997/98 P2 hindcasts; e the pattern correlation coefficients between hindcasts and observations for each ensemble (whiskers) and ensemble mean (dots) during 2015/16 (blue) and 1997/98 (red). The stippling indicates regions where the amplitude of the anomalies exceeds one standard deviation in observations, and significant at 5% level for P1 and P2 hindcasts

The convergence of the NMME prediction system to a teleconnection pattern opposite to that observed provides a test case to explore sources of predictability. One hypothesis is that the distinct SST spatial pattern during 2015/16—which was different from that predicted by the NMME—drives a different teleconnected precipitation pattern over the WUS. Alternatively, the fundamental difference between the regional pattern of precipitation anomalies over the WUS during 2015/16 and the canonical El Niño pattern might be reconciled if atmospheric internal variability could overwhelm the ENSO-induced pattern in winter WUS precipitation, even for a strong ENSO year. Generally, the seasonal predictability of atmospheric internal variability is substantially smaller than the predictability arising from coupled ocean–atmosphere modes of variability (like ENSO). Thus, atmospheric noise poses a great challenge to the accurate prediction of regional precipitation, even though ENSO and its global teleconnection patterns during winter can be predicted up to 9 months in advance (Jia et al. 2015; Yang et al. 2015a). Nevertheless, skillful seasonal prediction arising from the atmospheric state has been shown in some cases of summer heat waves and winter Arctic Oscillation for a one-season horizon (e.g., Jia et al. 2016, 2017). At shorter timescales, the influence of land initial conditions can be felt for weeks to months (Koster et al. 2004, 2011; Guo et al. 2012; Paolino et al. 2012). The intraseasonal low-frequency atmospheric oscillations with periods of 20–70 days (Ghil and Mo 1991; Plaut and Vautard 1994; Marcus et al. 1994) can also modulate short-term seasonal climate prediction. Therefore, to investigate the predictability source of the 2015/16 WUS winter precipitation, we focus on seasonal retrospective forecasts for a one-season horizon (i.e., initialized on 1-December to predict December, January and February) using one of the forecast systems included in the NMME using both its original NMME configuration and a new one designed to explore potential improvements for the hydroclimate prediction.

The failure of individual models in the NMME might arise from different factors that may relate to the model itself (e.g., model physical parameterization, resolution) or the initialization method. Given our available resources, the focus of this paper is to provide an in-depth assessment of precipitation predictive skill for the 2015/16 winter as a case study using one of the Geophysical Fluid Dynamics Laboratory (GFDL) operational seasonal forecast models in the NMME (Vecchi et al. 2014; Kirtman et al. 2014). In this study, we use a suite of targeted high-resolution experimental retrospective forecasts to explore the roles of atmosphere and land initial states in addition to the ocean state in predicting the unusual 2015/16 western United States precipitation pattern, characterized by a “flipped” pattern of precipitation anomalies (i.e., wet anomalies occur in the northern part of the WUS in contrast to wet anomalies in the southern WUS that are seen in the canonical state associated with ENSO). Details of model simulations, observational datasets and methods are discussed in Sect. 2. The prediction results are summarized in Sect. 3. The physical mechanism of the flipped precipitation patterns between 1997/98 and 2015/16 will be discussed in Sect. 4. Conclusions and discussions are given in Sect. 5.

2 Model simulations, observational data and methods

2.1 Models

We used the GFDL Forecast-oriented Low Ocean Resolution model (FLOR, Vecchi et al. 2014). The atmosphere and land components of FLOR are taken from the Coupled Model version 2.5 (CM2.5, Delworth et al. 2012) developed at GFDL, whereas the ocean and sea ice components are based on the GFDL Coupled Model version 2.1 (CM2.1, Delworth et al. 2006; Gnanadesikan et al. 2006; Wittenberg et al. 2006). FLOR has a spatial resolution of ~ 50 km in the atmosphere and land, ~ 100 km in the ocean, and 32 (50) vertical levels in the atmosphere (ocean). In this study, we used the flux-adjusted version of FLOR (Vecchi et al. 2014), in which climatological adjustments were made to the model’s momentum, enthalpy, and freshwater fluxes from the atmosphere to the ocean to bring the model’s long-term climatology of sea surface temperature (SST) and surface wind stress closer to that of observations (Magnusson et al. 2013).

FLOR is one of GFDL’s operational seasonal forecast models in the NMME (Kirtman et al. 2014) and provides routine seasonal forecasts to the public every month. It shows skillful seasonal prediction of tropical cyclone activity (Vecchi et al. 2014), near-surface temperature and precipitation over land (Jia et al. 2015) as well as extratropical storm activity (Yang et al. 2015a). The high spatial resolution of FLOR offers the potential for predicting regional hydroclimate over the WUS, particularly over regions of complex topography as it more faithfully resolves orography and precipitation extremes than lower resolution models (van der Wiel et al. 2016).

2.2 Retrospective experiments

We conducted three different sets of seasonal retrospective forecasts using FLOR; experiments are summarized in Table 1. They are used to explore various sources of prediction skill and highlight when and how initial conditions matter for improving the seasonal prediction of precipitation during the 1997/98 and 2015/16 El Niño events.

Table 1 Experiments used in this study

2.2.1 AMIP simulations

The first set of experiments consists of a 15-member atmosphere/land coupled AMIP simulations using FLOR’s atmosphere/land components. The model was integrated from 1982 to 2016, forced with observed time-varying sea surface temperature (SST) and radiative forcing. The SST data were taken from a merged product (Hurrell et al. 2008) based on the monthly mean Hadley Centre sea ice and SST dataset version 1 (HadISST1, Rayner et al. 2003) and version 2 of the National Oceanic and Atmospheric Administration (NOAA) optimum interpolation SST (OISST) analysis (Reynolds et al. 2002). Before 2005, the radiative forcing is an observationally based estimate of changing concentrations of greenhouse gases, aerosols, land-use changes, solar irradiance variations, and volcanic aerosols. After 2005, the radiative forcing is based on the representative concentration pathway (RCP) 4.5 scenario (Meinshausen et al. 2011). These simulations follow the AMIP’s standard protocol (http://www.pcmdi.llnl.gov/projects/amip/NEWS/overview.php). The AMIP simulations are used for assessing the response of the precipitation anomalies over the WUS to perfect SST conditions. A 3-year spin-up was performed during 1979–1981 prior to the start of the simulations (1st January 1982). The initial conditions of 15-member ensembles were generated from different days (i.e., day 1–15) in January 1979.

2.2.2 Retrospective seasonal forecasts

In the second set, retrospective seasonal forecasts were conducted for FLOR (called phase-1 forecasts, hereafter referred to as P1) from 1982 to 2016. The initial conditions of the ocean and ice components for FLOR forecasts were taken from the GFDLs ensemble coupled data assimilation (ECDA) system developed originally for CM2.1 (Zhang et al. 2007; Chang et al. 2013) but used now for FLOR given the common oceanic/ice model composition. The initial conditions for the atmosphere and land components were taken from the set of AMIP simulations described in Sect. 2.2.1. The radiative forcing used in P1 is identical to that in AMIP simulations. The retrospective forecasts were initialized from 1st December each year from 1982 to 2015, and integrated forward for 12 months with 12 ensemble members, however this study only focuses on the first 3 months of integration (December–February). Each ensemble member was initialized with the ocean and sea ice initial conditions from each member of CM2.1 ECDA. The atmospheric and land initial conditions were taken from three separate AMIP simulations. The first AMIP member supplied atmosphere and land initial conditions for the first four ocean/ice ensemble members, the second AMIP member supplied atmosphere and land initial conditions ocean/ice members from 5 to 8, and the third AMIP member supplied atmosphere and land initial conditions ocean/ice members from 9 to 12 (Vecchi et al. 2014).

In the third set, ensemble retrospective forecasts (called phase-2 forecasts, hereafter referred to as P2) were conducted with identical configurations as P1, except that the initial conditions of atmospheric and land components were taken from a set of simulations with the FLOR coupled model in which surface pressure, three-dimensional temperature and horizontal winds were nudged (relaxed) toward a reanalysis data set on a 6-hourly time scale. Version 1 of the Modern-Era Retrospective Analysis for Research and Applications (MERRA) reanalysis (Rienecker et al. 2011) was used. In addition, sea surface temperature (SST) was restored to time-varying observations on a 5-day time scale. The relaxation/nudging process started from 1979 to the present, and the outputs from the nudging runs were used as atmospheric and land initial conditions for P2 forecasts. The atmospheric and land initial conditions from the relaxation process are identical for all 12 members, but the ocean and sea ice initial conditions are different for each member produced by ECDA.

During the course of this study, the MERRA-2 reanalysis became available (Gelaro et al. 2017), so we repeated the same initialization experiments of P2, but using MERRA-2 data instead of MERRA. It is called P2M2 hereafter. Due to the computational resource limitation, only two hindcasts initialized on 1st December 1997 and 2015 were conducted for P2M2. MERRA-2 is an upgraded version of its predecessor MERRA, and includes updates to the Goddard Earth Observing System (GEOS) model, analysis scheme and observational data. The improvements in MERRA-2 include a better representation of the stratosphere, a reduction of some spurious trends and jumps related to changes in the observing system, and reduced biases and imbalances in aspects of the water cycle (Gelaro et al. 2017).

These prediction experiments were developed to enhance our understanding of prediction system skill using different initialization techniques and a single global climate model. The methodology presented is useful for guiding future prediction system development. To summarize: (1) the P1 methodology is used for the initial condition setup for FLOR used in the NMME (Vecchi et al. 2014; Kirtman et al. 2014). (2) P2 has been developed for research purposes to explore sources of prediction skill for various regional climate variations and extremes (Jia et al. 2016, 2017). (3) P2M2 has been developed for research purposes to understand how upgrades in reanalysis products may improve hindcasts of the two case study years. (4) P1 and P2 both use the flux-adjusted version of FLOR, and they only differ in the initial conditions of atmosphere and land components.

2.3 Observational dataset

The observationally constrained estimates of 850-hPa meridional wind, 700-hPa zonal and meridional winds, humidity, and geopotential height were obtained from the ERA-Interim reanalysis data for 1982–2016 (Dee et al. 2011). Precipitation data for 1982–2016 were obtained from PRISM (Daly et al. 2008) and the North American Land Data Assimilation System project phase 2 (NLDAS-2, Xia et al. 2012). Soil moisture content data for 1982–2016 was obtained from the Global Land Data Assimilation Systems Noah V2.7 (Rodell et al. 2004). The outgoing longwave radiation (OLR) data was obtained from NOAA/the National Centers for Environmental Prediction (Liebmann and Smith 1996).

2.4 Extratropical storm activity measurement

The extratropical storm index (ETSI) is defined as the seasonal standard deviation of filtered daily 850-hPa meridional wind (v850), and the filter is a 24-h-difference filter (Wallace et al. 1988). The ETSI can be written as follows:

$$ETSI=\sqrt {\frac{1}{{N - 1}}\mathop \sum \limits_{{t=1}}^{N} {{\left[ {v850\left( {t+24{\text{h}}} \right) - v850(t)} \right]}^2}} ,$$
(1)

where N is the sample size of the December-February (DJF) for each year and t represents a time step of 24 h. This method of measuring ETS has been widely used in the future projections of ETS (Chang 2013), the seasonal climate prediction of ETS variations (Yang et al. 2015a), and for attribution studies of the 2013/14-winter ETS extreme event over North America (Yang et al. 2015b).

2.5 Rossby wave source calculation

The Rossby wave source (S) (e.g., Sardeshmukh and Hoskins 1988) is calculated by:

$$S= - \nabla \cdot \left( {{v_\chi }\zeta } \right)= - \left( {\zeta \nabla \cdot {v_\chi }+{v_\chi } \cdot \nabla \zeta } \right),$$
(2)

where ζ is the absolute vorticity, \({v_\chi }\) is the divergent component of the horizontal winds. Following Scaife et al. (2017), the monthly data was used to calculate the Rossby wave source for observations and model simulations at the 200-hPa-pressure level. The absolute vorticity, divergence, and divergent winds were calculated using the standard scripts from the NCAR Command Language (Version 6.4.0, https://doi.org/10.5065/D6WD3XH5).

3 The seasonal predictions of the WUS precipitation during the two large El Niño events

3.1 Observational circulation anomalies

The observed 2015/16 WUS precipitation anomalies display a dipole pattern with wet anomalies to the north and dry anomalies to the south (Fig. 2b). In contrast, all coastal states (including California) had uniformly wet anomalies during 1997/98 (Fig. 2f). In general, the precipitation pattern during 1997/98 resembles the canonical El Niño pattern (Fig. 2a, f) except the far Northwest, while the pattern of anomalous precipitation during 2015/16 is nearly opposite (Fig. 2a, b).

To further explore dynamics associated with this precipitation pattern, we examine the atmospheric circulation fields, moisture and extratropical storm (ETS) anomalies during the 2 years. During 2015/16 DJF, a persistent low-level high-pressure anomaly was formed over the southeast North Pacific, extending to southern California (Fig. 3a). Consequently, due to blocking by the high-pressure anomaly, the anomalous winds with moister (drier) air flow from the ocean (land) to land (ocean) over northern (southern) California. In addition, the North Pacific storms shift northward due to this blocking, thus enhancing (reducing) ETS over northern (southern) California (Fig. 4a). In summary, the atmospheric circulation (the blocking high) and ETS anomalies favor the formation of wet (dry) precipitation anomalies over northern (southern) California that were observed (Fig. 2b). In contrast, during 1997/98 DJF, the observed low-level circulation anomalies display a dipole over North Pacific and Pacific coastal region with high (low) pressure anomalies in the south (north) resulting in onshore wind anomalies with moister air over California and the surrounding region (Fig. 3d). Consequently, the ETS anomalies tend to be enhanced over the southeast North Pacific (Fig. 4d). Thus, the circulation and ETS anomalies favor increased precipitation over the Pacific coastal region, including California (Fig. 2f).

Fig. 3
figure 3

The anomalous 700-hPa horizontal winds (vectors, m/s), specific humidity (shadings, g/kg), heights (contours, m) during the DJF season for 2015/16 a observations, b 2015/16 P1 hindcasts, c 2015/16 P2 hindcasts, d 1997/98 observations, e 1997/98 P1 hindcasts, and f 1997/98 P2 hindcasts. The sector (130°W–60°W) zonal mean is also removed from the height anomalies

Fig. 4
figure 4

The anomalous extratropical storm activities measured by the standard deviation of 24-h-difference-filtered daily 850-hPa meridional winds during the DJF season for a 2015/16 observations, b 2015/16 P1 hindcasts, c 2015/16 P2 hindcasts, d 1997/98 observations, e 1997/98 P1 hindcasts, and f 1997/98 P2 hindcasts

3.2 Seasonal prediction and predictability source

One hypothesis for the failed prediction of precipitation anomalies in 2015/16 is that during the 2015/16 winter, the observed SST pattern differed from what the models predicted, resulting in models predicting an atmospheric teleconnection response to a misrepresented SST (canonical El Niño) instead of the actual pattern of SST anomalies. Observed SST anomalies contrasting 2015/16 and 1997/98 are shown in Fig. 1: in the eastern tropical Pacific SST anomalies are substantially colder (up to 2°) and in the western Pacific, the SST anomalies are warmer (about 0.5°–1.0°) in 2015/16 versus 1997/98. We test the hypothesis for incorrectly predicted SST patterns influencing teleconnections with the 15-member ensemble AMIP runs with “perfect” observed SST values, and find that they cannot explain the difference between the WUS precipitation anomalies during the 2015/16 winter and those of a canonical El Niño (Figs. 2b, 5a). During the 1997/98 winter, the simulated anomalous precipitation pattern with perfect SST (Fig. 5) agrees to a large extent with the observed pattern as well as the patterns of model predictions when only ocean state was initialized by observations (P1) and when the entire climate system was initialized by observations (P2) (Fig. 2f, g, h). For 1997/98 the pattern correlation between observed and predicted precipitation over the WUS from AMIP is 0.81 (Fig. 2e).

Fig. 5
figure 5

The simulated precipitation anomalies driven by the observed SST from AMIP simulations. a The 2015/16 winter. b The 1997/98 winter. Simulations made with the atmospheric/land coupled models without coupling to the ocean. The stippling indicates regions where the amplitude of the anomalies is significant at 5% significance level

The winter simulations for 2015/16 provide a contrasting story of sources of prediction skill. The 2015/16 winter simulated precipitation pattern with perfect SST (Fig. 5a) substantially differs from the observed pattern (Fig. 2b) as well as the patterns of P1 and P2, with a pattern correlation coefficient with observations of − 0.25 (Fig. 2e). Note that the AMIP experiments can reproduce the observed reduction of precipitation in the southern part of WUS during the 2015/16 relative to the 1997/98, but not the wet anomalies in the northern part (Figs. 2, 5). Thus, the AMIP experiments reveal that the oceanic conditions dominated the source of regional precipitation anomalies over the WUS in 1997/98, while the oceanic conditions alone, even when forced to be “perfect”, cannot explain the observed precipitation anomalies in 2015/16. Therefore, it is either the case that: (1) our model is physically unable to recover the precipitation pattern in 2015/16, or (2) other climate subcomponents of the model initialization (e.g. atmosphere, land) beyond the ocean state were behind the observed anomalies. We test the second hypothesis with the P2 experiments where “correct” atmospheric and land initial conditions (IC) are used in addition to the ocean IC to initialize the seasonal prediction simulations.

Differences in seasonal prediction skill using the two IC setups (P1-ocean vs. P2-ocean/atmosphere/land) are clearly illustrated in Fig. 2. In 2015/2016, P1 predicted a general dipole pattern with dry (wet) anomalies in the northern (southern) part. This dipole pattern (Fig. 2c) is very similar to the canonical pattern (Fig. 2a), but opposite to the observed pattern found in Fig. 2b (the pattern correlation coefficient between ensemble mean prediction and observation is − 0.53). In contrast, P2 largely predicted the observed precipitation pattern (Fig. 2d) with a pattern correlation coefficient of 0.58. The sharp contrast of predictions between P1 and P2 indicates that the atmospheric/land IC play key roles in predicting the observed 2015/16 precipitation anomalies (Fig. 2e), despite the fact that 2015/16 was a very strong El Niño year.

For 1997/98, the predicted WUS precipitation anomaly patterns from P1 and P2 both agree to a large extent with the observed pattern as well as the canonical pattern (Fig. 2f, g, h); the pattern correlation coefficient is 0.76 and 0.81 for P1 and P2 respectively. The large similarity of predictions between P1 and P2 suggest that the oceanic ICs dominate over the atmospheric/land ICs in predicting the observed 1997/98 precipitation anomalies. Note that while the ensemble mean pattern correlation is similar in magnitude, the atmospheric/land IC substantially reduces the ensemble spread in P2 relative to that in P1 (Fig. 2e).

To investigate the dynamics underlying the predictions, we further examine the atmospheric circulation, moisture and ETS anomalies in P1 and P2. Interestingly, for the 2015/16 winter, P2 largely predicted the observed low-level blocking high over the southeast North Pacific with an extension to southern California (Fig. 3a, c). Consequently, P2 also largely predicted the observed wind-moisture pattern that the anomalous winds with moister (dryer) air flowed from ocean (land) to land (ocean) over northern (southern) California, due to the accurate prediction of the blocking high. Consistently, P2 predicted the enhanced (reduced) ETS activities over northern (southern) California (Fig. 4a, c). In contrast, P1 was not able to predict the blocking high over the southeast North Pacific and thus failed to predict the observed wind-moisture-ETS anomalous patterns (Figs. 3b, 4b). Thus, the skillful prediction of the 2015/16 precipitation anomalies is attributable to the atmospheric/land initial conditions.

For 1997/98, P1 and P2 both largely predicted the observed south-north dipole of 700-hPa height anomalies over the North Pacific and surrounding coastal region where low-level wind anomalies with moister air were directed from the ocean to California (Fig. 3d–f). Consistently, the increased ETS over North Pacific and California were well predicted by both P1 and P2 (Fig. 4d–f). The large similarity of predicting the atmospheric circulation and ETS anomalies between P1 and P2 further indicates that the oceanic conditions determine the formation of the 1997/98 precipitation anomalies.

To further test if the predictive skill of the WUS winter precipitation anomalies during 1997/98 and 2015/16 is sensitive to the verification data, we repeat the above analysis using the precipitation data from the North American Land Data Assimilation System project phase 2 (NLDAS-2, Xia et al. 2012). The predictive skill verified against NLDAS-2 precipitation data is similar to that against PRISM precipitation (Figs. 2, 6), indicating that FLOR’s predictive skill of WUS precipitation is robust.

Fig. 6
figure 6

The winter (DJF) precipitation anomalies (mm/day) for NLDAS-2 observations. a The observed canonical El Nino composite; b 2015/16 observations; d 1997/98 observations. c The pattern correlation coefficients between forecasts and observations for each ensemble (whiskers) and ensemble mean (dots) during 2015/16 (blue) and 1997/98 (red). AM in c denotes AMIP experiments. The stippling indicates regions where the amplitude of the anomalies exceeds one standard deviation in observations

We have extended our analysis beyond the two case study years to further test whether improvements from the P2 system are realized throughout the full set of NMME hindcasts. P2 exhibits systematic superior performance over P1 in predicting WUS winter precipitation during 1982–2016 (Fig. 7), suggesting that appropriate initialization of atmosphere/land, in addition to the ocean, could substantially improve the seasonal prediction of WUS precipitation. Interestingly, P2 also shows higher correlation skill than P1 in predicting global SST during 1982–2016 (Fig. 8). The considerable improvement with the correlation skill increase of over 0.2 is found over the Southern Ocean, South Pacific, North Pacific and North Atlantic. The improvement of the SST predictive skill in P2 is consistent with the atmosphere–ocean coupling process that atmosphere forcing is the dominant driving force of ocean anomalies in the midlatitudes on seasonal timescales (Wu et al. 2006; Wu and Kirtman 2007; Infanti and Kirtman 2017). Thus the addition of observed atmospheric forcing during initialization can improve SST prediction. Therefore, the addition of atmosphere/land initial conditions using observations to the existing ocean initial conditions enhances the seasonal prediction skill for the FLOR coupled prediction system.

Fig. 7
figure 7

The anomaly correlation coefficients (ACC) between retrospective forecasts and observations for the winter precipitation anomalies during 1982–2016 from P1 (top) and P2 (bottom). The stippling indicates ACC significant at 5% level. The verification data is from PRISM. The ratio of the grid points with positive significant correlation to the total grid points is 17 and 55% for P1 and P2 respectively

Fig. 8
figure 8

The anomaly correlation coefficients (ACC) between retrospective forecasts and observations for the winter SST anomalies during 1982–2016 from a P1 and b P2. The ACC difference between P2 and P1 is shown in c. Only the ACC at 5% significance level shown in a and b. The verification data is from OISST

4 The mechanism of the flipped precipitation pattern

In this section, we use model simulations and seasonal hindcasts to understand the possible physical drivers of the flipped precipitation patterns between the 2015/16 and 1997/98 winter seasons. To demonstrate the robustness of P2’s skillful predictions for the WUS winter precipitation, we also examined the P2M2 hindcasts, in which the MERRA2 data sets were used as an alternative input for initialization (see Sect. 2.2.2).

Figure 9 shows the precipitation difference between 2015/16 and 1997/98 for observations, ensemble means of AMIP simulations and hindcasts (P1, P2 and P2M2). The observed differential precipitation pattern displays a dipole with dry (wet) anomalies in the southern (northern) part of the WUS. The AMIP simulation can reproduce the observed dry anomalies in the southern part, but the amplitudes are substantially weaker than the observations. P1 predicted dry anomalies in the central part of the WUS. Interestingly, P2 and P2M2 successfully predict the amplitudes and positions of the dry (wet) anomalies in the southern (northern) part of the WUS.

Fig. 9
figure 9

The difference of the winter (DJF) precipitation anomalies (mm/day) between 2015/16 and 1997/98 for a observations, b AMIP simulations, c P1 hindcasts, d P2 hindcasts, and e P2 hindcasts with MERRA2 data. The stippling indicates regions where the amplitude of the anomalies exceeds one standard deviation in observations, and significant at 5% level for AMIP, P1 and P2 hindcasts

The ensemble distributions for the regional mean precipitation in the southern (northern) part (defined in the boxes of Fig. 9a) are shown in Fig. 10 for each experiment. The ensemble mean values of AMIP, P1, P2 and P2M2 agree on the dry anomalies in the southern part, although there are more than 25% of the members showing wet anomalies in AMIP and P1. For the northern part of WUS, the ensemble distributions of AMIP and P1 are similar to the climatological distribution, indicating there is no skill of predicting the wet anomalies for both AMIP and P1. For P2 and P2M2, the ensemble mean values are very close to the observed ones for the two wet and dry regions, and the signs of over 75% of the members agree with that of the ensemble mean prediction, indicating that the predictions of the flipped pattern using P2 and P2M2 are statistically reliable. The high similarity of the predicted precipitation patterns between P2 and P2M2 demonstrates the robustness of the P2 forecast system.

Fig. 10
figure 10

The box and whisker plots for the differential precipitation anomalies between 2015/16 and 1997/98 winters for a the southern part of the WUS (marked as red box in Fig. 9) and b the northern part of the WUS (marked as green box in Fig. 9) from AMIP, P1, P2 and P2M2 ensemble simulations. The ensemble mean values are denoted as blue circles, and the observed values are denoted as red asterisks

To examine if the different SST patterns in the tropics could drive the teleconnection precipitation pattern changes in the WUS between the two large El Niño years, the OLR differences between 2015/16 and 1997/98 are shown in Fig. 11 for the observations and model predictions. In observations, the OLR values during 2015/16 are significantly larger (smaller) than 1997/98 in the eastern and southern (western and northern) equatorial Pacific, indicating that deep convection is substantially weaker (stronger) in the eastern and southern (western and northern) equatorial Pacific in 2015/16 relative to 1997/98. The differential OLR pattern is consistent with the observed differential SST pattern such that the cold (warm) SST anomalies suppress (enhance) deep convection in the eastern (western) tropical Pacific. P1 and P2 largely predict the observed OLR and SST patterns, although the ensemble mean amplitudes are generally weaker than observations. Note that the amplitudes of differential OLR values in P2 are qualitatively closer to observations than those in P1. However, the AMIP simulated OLR anomalies are opposite to observed anomalies in the far western tropical Pacific and southern central tropical Pacific. This is consistent with previous studies that AMIP simulations cannot simulate the observed SST and OLR relationships in the tropical western Pacific (Wu et al. 2007).

Fig. 11
figure 11

The difference of the winter (DJF) OLR (W/m2, contours) and SST (°C, shadings) anomalies between 2015/16 and 1997/98 for a observations, b P1 hindcasts, c P2 hindcasts, and d AMIP simulations

The observed differential OLR pattern with substantially suppressed (enhanced) deep convection in the eastern and southern (western and northern) equatorial Pacific indicates that tropical circulation and its teleconnection pattern would change as a response to changes in convection-driven diabetic heating. Figure 12 shows the anomalous 200-hPa zonal winds and stream functions during the two winters for observations, AMIP simulations and P1/P2 hindcasts, respectively. During the 1997/98 winter, the observed pattern shows strong easterly anomalies in the central and eastern Pacific and a pair of anomalous anticyclones in the southern and northern off-equatorial regions at the 200-hPa pressure level (Fig. 12a), resembling a Gill-type response to the increased heating anomalies in the eastern tropical Pacific (Gill 1980; Jin and Hoskins 1995). The P1, P2 hindcasts and AMIP simulations for 1997/98 can reproduce the strong easterly wind anomalies in the equatorial Pacific and the two anomalous off-equatorial anticyclones (Fig. 12c, e, g), although the northern anticyclone shifts towards eastward in comparison with observations in the AMIP simulation.

Fig. 12
figure 12

The winter (DJF) 200-hPa stream function (106 m2 s−1, contours) and zonal wind (m/s, shadings) anomalies for 1997/98 (left panels) a observations, c P1 hindcasts, e P2 hindcasts, and g AMIP simulations; for 2015/16 (right panels) b observations, d P1 hindcasts, f P2 hindcasts and h AMIP simulations

The 2015/16 winter shows an altered convection and flow pattern relative to 1997/98. During the 2015/16 winter, the center of the easterly anomalies shifts towards the central Pacific and the amplitudes tend to be weakened about 50% relative to those of 1997/98 winter (Fig. 12b). Consistently, the two anomalous anticyclones shrink towards the central Pacific and their amplitudes are also reduced to about 50% relative to those of 1997/98 winter. The P1, P2 hindcasts and AMIP simulations for 2015/16 can simulate the reduced amplitudes of the equatorial easterly wind anomalies and the off-equatorial anomalous anticyclones relative to those found in 1997/98 (Fig. 12d, f, h). Remarkably, the P2 hindcast outperforms P1 and AMIP in terms of simulating the observed spatial patterns of the equatorial easterly wind anomalies and off-equatorial anticyclones.

To make a qualitative comparison of the ENSO-teleconnection circulation patterns between the 2 years, the differential zonal winds and stream functions at the 200-hPa level are shown in Fig. 13. Corresponding to the suppressed (enhanced) deep convection in the eastern and southern (western and northern) equatorial Pacific (Fig. 11a), the strong westerly wind anomalies (over 10 m/s) were formed in the tropical eastern Pacific in association with the cyclonic flow anomalies in the subtropical eastern Pacific (Fig. 13a). Consequently, the subtropical jet extending from the north Pacific to North America was substantially weakened (over 10 m/s) during 2015/16 winter relative to 1997/98. The weakening of the subtropical jet leads to the reduction of the extratropical storms in the southeastern North Pacific and southern California (Fig. 4a, b) during 2015/16 relative to 1997/98. Thus, the relative dry anomalies between the two El Niño winters in the southern WUS can be attributed to the difference of the SST spatial pattern and the consequent deep convection pattern in the tropical Pacific. The reduced (enhanced) deep convections in the eastern and southern (western and northern) tropical Pacific leads to the weakening of the subtropical westerly jet over the southeast North Pacific and southern WUS, resulting in the reduction of extratropical storms and precipitation during the 2015/16 winter.

Fig. 13
figure 13

The difference of the winter (DJF) 200-hPa stream function (106 m2 s−1, contours) and zonal wind (m/s, shadings) anomalies between 2015/16 and 1997/98 for a observations, b P1 hindcasts, c P2 hindcasts, and d AMIP simulations

The three models can largely simulate the observed westerly wind anomalies at the 200-hPa level in the eastern tropical Pacific, cyclonic flow anomalies in the southern and northern off-equatorial regions and the weakening of the subtropical jet (Fig. 13b–d). The pattern correlation coefficient of the 200-hPa stream function between model and observation over the domain (10°S–50°N, 140°E–90°W) is 0.67, 0.81, and 0.68 for P1, P2 and AMIP respectively, thus P2 outperformed P1 and AMIP in terms of simulating the amplitudes and spatial patterns of the ENSO-teleconnected circulation anomalies. The three models agree on simulating the weakening of the subtropical jet over the southeastern North Pacific and California, indicating that the dry anomaly in the southern part of the WUS is mainly due to the tropical El Niño teleconnection change between the two winters.

Interestingly, only P2 can predict the observed anomalous anticyclone over the northeastern North Pacific, in association with westerly wind anomalies in the northern part of the United States (Fig. 13a, c). The anomalous anticyclonic circulation and the westerly wind anomalies favor moister flow and more storms into the northern WUS, thus resulting in the wet anomalies in the northern WUS during the 2015/16 winter. Although the predictability source of the wet anomalies in the northern WUS in P2 is attributable to the atmospheric/land initialization, the physical mechanism of predicting the wet anomalies requires further study and should be the focus of future work.

To further elucidate the dynamical relationship between the tropical Pacific deep convection anomalies and the circulation anomalies during the 1997/98 and 2015/16 winters, we examine the relative difference of observed 200-hPa divergent winds, divergence and Rossby wave source (RWS) between the two winters (Fig. 14). Consistent with the differential OLR pattern with substantially suppressed (enhanced) deep convection in the eastern and southern (western and northern) equatorial Pacific (Fig. 11a), a general dipole differential divergence pattern with convergent (divergent) flow in the eastern and southern (western and northern) equatorial Pacific is formed at the 200-hPa level (Fig. 14c). Consequently, the RWS was substantially reduced (enhanced) in the equatorial (polar) flank of the subtropical jet in the eastern Pacific-North America sector (10°N–40°N, 150°W–100°W) during 2015/16 relative to 1997/98 (Fig. 14d–f). The observed differential divergence pattern in the central and eastern tropical Pacific (Fig. 14c) was well predicted by P1 and P2. Remarkably, the differential RWS pattern in the eastern tropical Pacific-North America sector was well predicted by P2 in terms of the pattern loading amplitudes and signs (Figs. 14f, 15e). The predicted differential RWS pattern by P1, however, shifts northward and with amplified negative loadings in the equatorial part compared with observations (Figs. 14f, 15d). The AMIP simulation only captured the observed east–west contrast of the divergence field in the tropical Pacific, but failed to reproduce the observed south-north contrast of the divergence filed (Figs. 14c, 15c). Consequently, the AMIP model only simulated the observed negative RWS anomalies extending from the Baja California Peninsula coast towards the central tropical Pacific, but failed to capture the observed positive RWS anomalies in the Pacific west coast (Figs. 14f, 15f). Thus, the P2 hindcast outperforms P1 and AMIP in terms of simulating the observed differential patterns of the divergence and RWS anomalies between 2015/16 and 1997/98 winter in the tropical Pacific and North American sector.

Fig. 14
figure 14

The observed winter (DJF) 200-hPa divergence (10−6 m2 s−1, shading) and divergent winds (m/s, vector) anomalies (left panels) a 1997/98, b 2015/16 and c the difference between 2015/16 and 1997/98; the observed winter (DJF) 200-hPa Rossby wave source anomalies (10−11 s−2, shading) (right panels) for d 1997/98, e 2015/16 and f the difference between 2015/16 and 1997/98. The contours denote the zonal wind (m s−1) in d and e

Fig. 15
figure 15

The simulated difference of winter (DJF) 200-hPa divergence (10−6 m2 s−1, shading) and divergent winds (m/s, vector) between 1997/98 and 2015/16 (left panels) for a P1, b P2 and c AMIP; the simulated difference of winter (DJF) 200-hPa Rossby wave source (10−11 s−2, shading) between 1997/98 and 2015/16 (right panels) for d P1, e P2 and f AMIP

5 Conclusions and discussions

Our observational analysis shows that the presence of a blocking high over the southeast North Pacific favors the formation of a precipitation pattern over the WUS, which is nearly opposite to the canonical precipitation pattern associated with El Niño. We call this a “flipped” pattern in reference to the north–south reversal of precipitation anomalies relative the canonical El Niño pattern. When the model is only initialized with ocean observations (P1), it fails to predict the observed precipitation pattern over the WUS. Remarkably, when we initialize the model with ocean, atmosphere and land observations (P2), WUS precipitation is more accurately predicted in 2015/16 and for the entire hindcast period (1982–2016) relative to when only the ocean is initialized (P1). The contrasting outcomes in prediction highlight that during the 2015/16 winter, atmospheric internal variations were superposed with the broad influence of the global-scale ENSO-teleconnected wave pattern; hence the skillful prediction of such events has to be linked to both oceanic, atmospheric and land initial conditions. Note that weather and subseasonal climate predictability contributes to the seasonal mean prediction skill found in our prediction system in this study. However, the attribution of weather, individual events, and subseasonal predictability to the seasonal predictability is beyond the scope of this study. The analyzed source of seasonal prediction skill may depend on the forecast model and the initialization methodology used in this study.

The 2015/16 winter season provides an important test case for prediction systems and reveals that a canonical precipitation pattern may not occur over the WUS even during a strong El Niño. The 1997/98 winter was the last major El Niño with very strong teleconnection patterns that were robustly predicted in retrospective simulations—they succeeded since the precipitation pattern during 1997/98 winter was mainly controlled by the ocean state. When a similarly strong El Niño was predicted in late 2015 with canonical precipitation anomalies found in the NMME over the WUS, it provided the appearance of a repeat year with similar forecast patterns, giving false hope for a wet year when reported in the popular press (Kintisch 2016) to a region experiencing a multi-year drought. While opposite patterns could be found in individual ensemble members in our ocean-only experiment (P1), the ensemble mean also exhibited a canonical response. The ensemble mean precipitation anomaly flips however with the incorporation of the atmospheric/land states (P2), suggesting that the atmosphere/land states strongly interacted with the ocean state during the 2015/16 winter.

Observational comparison of the ENSO teleconnection patterns between the 2015/16 and 1997/98 winters reveals that the differential SST pattern with cold (warm) anomalies in the eastern (western) tropical Pacific drives a differential OLR pattern with suppressed (enhanced) deep convections in the eastern (western) tropical Pacific during the 2015/16 (1997/1998) winters. Consequently, the amplitudes of the El Niño-induced circulation anomalies (including the equatorial easterly wind anomalies and the off-equatorial anticyclone anomalies) during 2015/16 are weakened by 50% relative to those during 1997/98. This reduced teleconnection leads to a weakening of the subtropical westerly jet over the southeastern North Pacific and southern WUS, thus resulting in the reduction of extratropical storms and precipitation in the southern WUS, including California, during the 2015/16 winter. P2 outperforms both P1 and AMIP for predicting the difference between the two case study years in terms of the tropical circulation elements (e.g., divergence, stream function and Rossby wave source) and ENSO teleconnection in the tropical Pacific-North America sector, thus accurately predicting the dry anomalies in the southern WUS. The P1 and AMIP experiments also support the hypothesis that a physical driver of the observed relative dry anomalies in the southern WUS is due to tropical circulation changes between the two years. The relative wet anomalies in the northern WUS are associated with an anomalous anticyclone over the northeastern North Pacific and westerly wind anomalies in the northern United States. These circulation anomalies favor the enhancement of moist flow and more storms into the northern WUS, thus resulting in the observed 2015/16 winter wet anomalies. P2 predicted the observed anomalous anticyclone over the northeastern North Pacific, thus predicting the wet anomalies in the northern WUS and the whole flipped precipitation pattern in the WUS. The physical mechanism of predicting the wet anomalies requires further study. One possibility is that the predictability arises from low-frequency atmospheric internal variability, which interacts with the El Niño teleconnection in the WUS. Meanwhile, the predictability source could also originate from the tropics, such that the distinct 2015/16 teleconnection pattern (e.g., the deep convections suppressed in the eastern and southern tropical Pacific) can drive the flipped WUS precipitation pattern (evidence for this is found in the Rossby wave source analysis shown in Figs. 14, 15). In other words, the successful prediction of the wet anomalies in P2 could be due to the accurate prediction of the El Niño teleconnection in the tropics. It also can be speculated that the predictability source of the wet anomalies in the northern WUS in P2 is the combination of the low-frequency atmospheric internal variability and the tropical El Niño teleconnection.

It has previously been reported that there is little evidence of an improved forecast of precipitation over land due to initialization of the atmosphere/land (Paolino et al. 2012). However, this work shows that appropriate initialization of the atmosphere/land in addition to the ocean could substantially improve seasonal precipitation prediction. Incorporating the land and atmospheric states not only improves prediction during the 2015/16 winter, but has also been shown to improve precipitation and SST predictions over the entire period (1982–2016) over which hindcasts were performed (Figs. 7, 8). To attain better seasonal predictions of precipitation with broad societal benefits, prediction systems should be updated to incorporate the atmosphere/land states in addition to the ocean. The impact of the atmospheric and land initial conditions on predicting 2015/16 precipitation anomalies over the WUS is limited to only one season horizon, which is consistent with the atmospheric and land impact on the seasonal prediction of summer heat waves (Jia et al. 2016) and winter surface air temperature (Jia et al. 2017) over the United States. Although the skillful predictive window due to the atmospheric/land initial condition is within one season in advance, this will provide better climate forecasts to address water related issues, including water supply management and various water-related hazards.