Four datasets were employed in this study, each spanning a period of 40°years from 1979 to 2018, as follows:

(1) The atmospheric reanalysis data used in this study are from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR), with a 2.5 horizonal resolution from 1,000 to 10 hPa comprising 17 pressure levels (Kalnay et al., 1996), including the monthly mean SAT, geopotential height, and zonal and meridional winds.

The multivariate empirical orthogonal function (MV-EOF) analysis method (Wang 1992) is adopted here to distinguish the leading variability modes that combine the variations of the December and following-January SATAs of CA during 1979–2018. In this study, an area-weighted correlation coefficient matrix is constructed for the combined SATA (weighted by the square root of the cosine of latitude) over CA in December and the following January during 1979–2018 to carry out the MV-EOF. The method has been widely applied in climate research. For example, Li et al. (2021) pointed out that the MV-EOF method can distinguish the mode of the reversal of December and January SATAs over China.

The wave activity flux (Takaya and Nakamura 2001) is employed to describe the stationary Rossby wave propagation. This wave flux is parallel to the local group velocity of a stationary Rossby wave train in the Wentzel–Karmers–Brillouin approximation. The horizontal flux is calculated as follows: W = P 2 | U → | { U ( v ′ 2 − ψ ′ v x ′ ) + V ( v ′ u ′ − ψ ′ u x ′ ) U ( − u ′ v ′ − ψ ′ u x ′ ) + V ( u ′ 2 − ψ ′ u y ′ ) (1) where P is the pressure scaled by 1,000 hPa, (U, V) denotes the horizonal mean winds, ψ ′ is the perturbation geostrophic streamfunction, and ( u ′ , v ′ ) denotes the perturbed winds. The climatological mean values are calculated during the period 1979–2018.

The strength of the synoptic-scale transient eddy activity (STEA) is measured by v ′ 2 ¯ at 300 hPa, where v ′ denotes the synoptic-scale daily meridional winds subject to a 2.5–6-day band-pass filter, and the overbar represents averaging over a month (Ren et al., 2010). The eddy-induced geopotential height tendency defined in Lau and Nath (2014) is calculated as follows: ( ∂ Z ∂ t ) e d d y = f g ∇ − 2 [ − ∇ ( v → ′ ζ ′ ) ¯ ] (2) where Z is the monthly mean geopotential height; g is the gravitational acceleration; v → ′ = ( u ′ , v ′ ) and ζ ′ are the synoptic-scale daily winds and relative vorticity at 300 hPa, respectively, subject to 2.5–6-day bandpass filtering.

The Student’s t-test was used to estimate levels of statistical significance, and the linear trends of all data used for regression and correlation were removed. The number of equivalent degrees of freedom for regression was calculated following Zwiers and von Storch (1995).

An atmospheric general circulation model, ECHAM5 (Roeckner et al., 2003), is employed to validate the mechanisms given in this study. The ECHAM5 simulations were performed with the spectral T63 horizontal resolution and 19 vertical levels (T63L19). The detailed experimental design is given in section 3d.

Figure 2 shows the spatial patterns of the first two leading modes (MV-EOF1 and MV-EOF2) and corresponding time series (PC1 and PC2). The first and the second leading modes explain 30.7 and 22.3% of the total covariance, respectively, and are both well separated in accordance with the method given by North et al. (1982). MV-EOF1 is an in-phase variation mode, which has similar spatial patterns in December and January, with strong amplitudes situated over the middle and northern part of CA (Figure 2A,B). MV-EOF2, meanwhile, shows maximum amplitudes distributed along the northeast–southwest of CA in December and exhibits nearly opposite anomalies in the same region in January (Figure 2C,D), which presents the out-of-phase nature of the SATAs between December and January. Similar in-phase and out-of-phase modes can be separated using the GHCN dataset, which have nearly the same spatial distributions as MV-EOF1 and MV-EOF2 (not shown), and the corresponding time series are also highly correlated with PC1 and PC2 (higher than 0.97). Thus, the reversal of the SATAs between December and January over CA can be represented by MV-EOF2.

To examine the atmospheric circulation anomalies associated with MV-EOF2, the PC2-related SATAs and the 850 and 300 hPa geopotential height and wind anomalies in December and January are analyzed (Figure 3). The positive PC2 corresponds to the remarkable December warming and January cooling (Figures 3A,B). In December, anomalous southeasterlies can be seen over CA at 850-hPa, with a southwest‒northeast distribution pattern (Figure 3C). The distribution pattern of the anomalous southeasterlies corresponds well with that of the SATAs (Figure 3A,C), which indicates that the anomalous warm advection brought by the southeasterlies may contribute to the SATAs in December. Besides, significant southwest‒northeast positive height anomalies appear at 300 hPa (Figure 3E), which can induce a decline in cloud cover. Thus, enhanced downward shortwave radiation may also contribute to the warm SATAs over CA in December. In January, significant positive geopotential height anomalies appear over the north side of CA at both 850 hPa and 300 hPa (Figure 3E), accompanied by anomalous northeasterlies over CA. The cold advection brought by these anomalous northeasterlies may cause the cold anomalies over northern CA. The contrasting atmospheric circulation anomalies between December and January play an essential role in the reversal of the SATAs over CA.

Because MV-EOF2 represents an out-of-phase relationship between the December and January SATAs, the signs of the mean SATAs of the 2 months are more likely to be opposite when the amplitude of MV-EOF2 is “relatively stronger” than that of MV-EOF1. Years with opposite signs are marked in Figure 2E,F. A higher PC2 is likely to appear in these years. To verify this, the index R p c 2 = | P C 2 | / ( | P C 1 | + | P C 2 | ) is defined to represent the “relative strength” of MV-EOF2, and the standardized time series of R p c 2 is given in Figure 2G. The years in which the December and January average SATAs are of opposite sign are marked in Figure 2G. It can be seen that most of these years show positive standardized values. The average value for the years with opposites signs (16 years) is 0.67, while for the years with same signs (24 years) is −0.46, the difference between the two values is statistically significant at the 95% level according to the Student’s t-test.

To reveal the remarkable remote SST forcing that modulates the reversal of the SATAs between December and January over CA, the PC2-related SSTAs and outgoing longwave radiation (OLR) in December and January are analyzed (Figure 4). Significant negative SSTAs appear over CTP within the Niño4 region (5°S‒5°N, 160°‒150°W) in both December and January (Figure 4A,B), and the pattern of the SSTAs resembles that of CP La Niña (Yuan and Yan, 2013, see their Figure 2B). Furthermore, positive OLR anomalies appear in the Niño4 region in both December and January, and these positive OLR anomalies correspond well to the SSTAs (Figure 4C,D). This indicates that the PC2-related diabatic heating is directly associated with the SSTAs in the 2 months. As the Niño4 index is defined as the area-averaged SSTA over the Niño4 region (Kug et al., 2009), the −Niño4 index is adopted to represent the variation of the CTP SSTAs. The correlation coefficients between the −Niño4 index in December and January and PC2 are 0.37 and 0.38 (both statistically significant at the 99% level), respectively. Thought the correlation coefficients are not too high, the result still suggests that a significant correlation exist between the between the cold CTP SSTAs and the MV-EOF2.

Moreover, the −Niño4-related SATAs and circulation anomalies over CA are analyzed (Figure 5). In December, weak positive SATAs are located over northern and western CA, corresponding to the anomalous southerlies and southwesterlies over western and northern CA at both 850 and 300 hPa (Figures C and E). The warm advection by the anomalous winds may induce the weak positive SATAs in December. In January, negative SATAs appear over CA, corresponding to the significant negative geopotential height anomalies at 300 hPa and anomalous northeasterlies at both 850 and 300 hPa (Figures 5B,F), which indicates the cold advection by the anomalous winds may cause the cold anomalies over CA. Besides, the decline of the downward shortwave radiation caused by the negative geopotential height anomalies at 300 hPa may also contribute to the cold anomalies. Though differences exist between the PC2 and −Niño4-related warm SATAs over CA in December (Figures 3A, 5A), similarity at some level still exists between the PC2 and −Niño4-related anomalous geopotential height and winds. Furthermore, the abrupt change in the influence of the cold CTP SSTAs on the SATAs in December and January may contribute to the reversal of the SATAs over CA. The results also imply that the cold CTP SSTAs may not be the single factor for the reversal of SATAs over CA, however, it is out of the scope of this study but deserves further investigations in the future.

To explore the possible physical mechanisms behind the impact of remote SST forcing on the reversal of SATAs between December and January over CA, the −Niño4-related 300 hPa geopotential height anomalies and the wave activity fluxes in December and January are analyzed (Figure 6). In December, a wave train extends from the North Pacific to CA, crossing North America and the North Atlantic. This wave train consequently results in the positive geopotential height anomalies over CA (Figure 6A). In contrast, in January, the wave train originates from the mid-latitude North Pacific and mainly propagates in the meridional direction (Figure 6B). This wave train cannot reach CA to influence the local SATAs. Interestingly, a wave train starting from the North Atlantic could reach CA and result in the local negative geopotential anomalies and SATAs. Therefore, two questions are raised regarding the possible mechanisms underlying the different roles of the wave trains associated with the −Niño4 between December and January: 1) why does the propagation of the wave train from the North Pacific differ between December and January; and 2) how do the cold CTP SSTAs in January cause the negative SATAs over CA in January?

For the first question, the different Gill-type responses to the cold CTP SSTA between December and January may play a key role. The Gill-type response is the stationary Rossby wave response to the tropical SST anomalies (Gill, 1980). The Gill-type response corresponding to the cooling SST anomalies of La Niña winters can be characterized by a pair of low-level anticyclonic anomalies residing at the northern and southern sides of the equator, along with low-level easterly wind anomalies and cyclonic anomalies symmetrical about the equator at the upper level of the troposphere (Yuan and Yan, 2013; Zhang et al., 2015). In Figure 7, the Gill-type response to the CTP SSTAs is remarkable: on the one hand, anomalous easterlies and symmetric anticyclones at 1,000 hPa appear around the region (20°S‒20°N, 170°E‒120°W) in both December and January (Figures 7A,B), and the easterly anomalies are relatively stronger in January than in December; whilst on the other hand, in both December and January, a pair of cyclonic anomalies appears over the CTP at 300 hPa, and these anomalous cyclones in January are much stronger than they are in December (Figure 7C,D). The result implies that the Gill-type response to the cold CTP SSTAs in January is stronger than in December.

With the enhanced Gill-type response in January, the influence of cold CTP SSTAs on the extratropical atmosphere over the North Pacific are strengthened. In January, corresponding to the enhanced 300 hPa cyclonic anomalies in the region (0°‒20°N, 160°‒130°W), stronger positive geopotential height anomalies appear in the region (30°‒60°N, 180°‒150°W) (Figure 7C,D). Consequently, intensified wave activity fluxes are driven by the positive geopotential height anomalies, as seen in Figure 6. Also, the more remarkable meridional dipole structure of the geopotential height anomalies over the North Pacific might cause the meridional propagation of the Rossby waves.

Additionally, the reason why the Gill-type response to the cold CTP SSTAs in January is stronger than that in December is explored. Previous studies have revealed that a strengthening of the tropical heating could lead to a stronger Gill-type response (Xing et al., 2014). Besides, Lee et al. (2009) suggested that a stronger and more equatorward jet stream could result in a stronger Gill-type response in the tropical atmosphere, along with stronger responses in the extratropical atmosphere. Figure 8 shows the CTP SSTI-related SSTAs and the climatological 300 hPa zonal winds in December and January. The cold SSTAs in the CTP region in January are slightly stronger than these in December, but their difference is not statistically significant (Figure 8A–C). Additionally, the subtropical jet stream over the North Pacific in January is stronger and situated more towards the equator (Figure 8D–F), which may lead to the strengthened Gill-type response in the tropics, as indicated by Lee et al. (2009). Thus, the CTP SSTAs in December and January may excite different Gill-type responses in the tropical atmosphere due to the difference in the climatological jet stream. Subsequently, the different Gill-type and extratropical atmospheric responses finally lead to a different propagation of Rossby waves between December and January.

In terms of the second question, the strong 300 hPa zonal wind anomalies associated with the CTP SSTAs may play a crucial role. In January, significant positive and negative 300 hPa zonal wind anomalies appear in the regions (40°‒60°N, 180°‒120°W) and (20°‒40°N, 180°‒120°W), respectively (Figure 9A). The meridional structure of the anomalous zonal wind can be linked to the dipole structure of the geopotential height anomalies in the January (Figure 7D). Corresponding to the wind anomalies, significant positive and negative STEA anomalies appear over the eastern North Pacific (Figure 9A,B). These STEA anomalies can travel eastwards to the North Atlantic via the downstream development of wave packets (Rivière and Orlanski, 2007; Li and Lau 2012a, Li and Lau 2012b). Accordingly, the strength of STEA over the North Atlantic shows non-significant positive anomalies over (40°‒60°N, 60°‒30°W) and significant negative anomalies over (20°‒40°N, 60°‒30°W) (Figure 9B). Such a dipole structure of STEA anomalies facilitates a positive geopotential height tendency over the region (40°‒50°N, 60°‒30°W) (Lau and Holopainen 1984), which may lead to the positive geopotential height anomalies associated with the −Niño4 index (Figure 9C). Hence, the CTP SSTAs may influence the STEA over the North Atlantic, by which the positive geopotential height over the North Atlantic can be forced, and then the Rossby waves excited by the positive geopotential height can propagate to CA to cause the local cold SATAs in January. Moreover, the positive geopotential height over northwestern CA forced by the anomalous STEA may also contribute to the driving of the Rossby waves and the cold SATAs.

In this section, three experiments were designed to validate the mechanisms given above. The first one forces the model with the climatological seasonal cycle of global SST and sea ice (referred to as Ctrl). Ctrl was integrated for 50°years, and the first 20°years were used for model spin-up and the remaining years were used for analysis. In the second experiment (referred to as Run1), the boundary forcing was changed to the climatological seasonal cycle of global SST with additional SSTAs over the tropical Pacific region (10°S‒10°N, 150°E‒90°W) from December to the following January. The SSTAs used here were the regressions of the December (January) SST onto the simultaneous −Niño4 index during 1979–2017 (1980–2018), and keeping the same value in December (January). Run1 included 30 members. Each member was integrated from 30 November to 31 January of the following year with the initial conditions taken from Ctrl. The third experiment (referred to as Run2) was similar to Run1, except that the SSTAs added to the climatological seasonal cycle of the global SST in December and January were both from the regression of the December SST onto the −Niño4 index during 1979–2017.

The composite differences in the SAT and 300 hPa geopotential height between Run1 and the control run are given in Figure 10. In December, although the SATAs over CA in Run1 are weak and non-significant (Figure 10A), the geopotential height in December is similar to the −Niño4 index-related geopotential height anomalies (Figure 10C). The Gill-type geopotential height anomalies appear over the tropical Pacific and match with the obvious result. Furthermore, the positive anomalies over the Bering Strait, the negative anomalies over North America, and the weak negative anomalies over western Europe in Run1 correspond to the −Niño4 index-related geopotential height anomalies. Thus, the numerical experiments roughly reproduce the mechanism of the −Niño4 index’s impacts on CA in December. In January, strong and significant negative SATAs appear over CA (Figure 10B), and the geopotential height anomalies are similar to the observational anomalies (Figure 10D). Interestingly, the Gill-type geopotential height anomalies over the tropical Pacific are stronger than in December, as expected. The results from Run1 in December and January support the proposed mechanism. On the one hand, in December, there is a similarity between the geopotential height anomalies and the −Niño4-related geopotential height anomalies, which validates the suggestion that an eastward-propagating Rossby wave train may be excited by the anomalous cold CTP SSTA in December. On the other hand, in January, a stronger Gill-type response in the tropical Pacific is excited by −Niño4-related SSTAs, and may cause significant negative SATAs over CA.

Since the SST forcings used in Run1 in December and January were different, as shown in Figure 8C, it is possible that the different responses in Run1 are caused by the difference in the SST forcings. Thus, in Run2, the same SST forcing [the −Niño4-related SSTAs over the region (10°S‒10°N, 150°E‒90°W) in December] was added in December and January. The composite differences in the SAT and 300 hPa geopotential height between Run2 and Ctrl are given in Figure 11. The results closely resemble those of Run1. The similarity between Run1 and Run2 provides further support to the notion that changing climatological background winds might play an important role in the reversal of the SATAs over CA.