## Abstract

Nonlinear principal component analysis (NLPCA), via a neural network (NN) approach, was applied to an ensemble of six 47-yr simulations conducted by the Canadian Centre for Climate Modelling and Analysis (CCCma) second-generation atmospheric general circulation model (AGCM2). Each simulation was forced with the observed sea surface temperature [from the Global Sea Ice and Sea Surface Temperature dataset (GISST)] from January 1948 to November 1994. The NLPCA modes reveal nonlinear structures in both the winter 500-mb geopotential height (Z500) anomalies and surface air temperature (SAT) anomalies over North America, with asymmetric spatial anomaly patterns during the opposite phases of an NLPCA mode. Only during its negative phase is the first NLPCA mode related to the El Niño–Southern Oscillation (ENSO); the positive phase is related to a weakened jet stream. Spatial patterns of the NLPCA mode for the Z500 anomalies generally agree with those for the SAT anomalies.

Nonlinear canonical correlation analysis (NLCCA), also via an NN approach, was then applied to the midlatitude winter GCM data and the observed SST of the tropical Pacific. Nonlinearity was detected in both the forcing field (SST) and the response field (Z500 or SAT) at zero time lag. The leading NLCCA mode for the SST anomalies is a nonlinear ENSO mode, with a 30°–40° eastward shift of the positive SST anomalies during El Niño relative to the negative SST anomalies during La Niña. The leading NLCCA mode for the Z500 anomaly field is a nonlinear Pacific–North American (PNA) teleconnection pattern. The ray path of the Rossby waves induced during El Niño is 10°–15° east of that induced during La Niña. The nonlinear atmospheric response to ENSO is also found in the leading NLCCA mode for the SAT anomalies.

## 1. Introduction

Substantial progress has been made during the past two decades toward understanding the wintertime extratropical atmospheric responses to the tropical forcings associated with El Niño sea surface temperature (SST) anomalies (Trenberth et al. 1998). The most prominent teleconnection pattern is the Pacific–North America (PNA) pattern (Horel and Wallace 1981), which is thought to link the changes in the extratropical circulation to the tropical SST through Rossby wave dynamics (Wallace and Gutzler 1981; Hoskins and Karoly 1981).

The simplest view of the atmospheric climate signal associated with the El Niño–Southern Oscillation (ENSO) phenomenon is that the atmosphere responds linearly, with anomalies during the El Niño phase being the reverse of those during the La Niña phase. However, recent evidence shows that the responses to warm and cold events are not exactly opposite. Sittel (1994) found the marginal probabilities of extreme rainfall and temperature over the southeastern United States to be highly nonlinear functions of the phase of the Southern Oscillation (SO). For example, it was found that, since 1946, the warm event–enhanced rainfall signal has been much larger than the cold event–suppressed signal in the southeastern United States. Richman and Montroy (1996) examined the composite January temperature and precipitation patterns over the United States and parts of Canada associated with El Niño and La Niña events. Their results suggest that El Niño and La Niña have their own unique characteristics in terms of temperature and precipitation. Asymmetric spatial patterns of Canadian surface air temperature and precipitation associated with the SO were also detected by Shabbar and Khandekar (1996) and Shabbar et al. (1997). Further evidence for a nonlinear response of the North America climate to ENSO was provided by Hoerling et al. (1997) with composite analysis and numerical experiments using an idealized atmospheric general circulation model, where a shift in the equatorial positions of the maximum rain responses and a phase shift of teleconnection patterns in the upper troposphere were found. The inherent nonlinearity in the tropical rain response may itself be responsible for the phase shift in the extratropical teleconnection patterns.

The robustness of the nonlinear climate response to ENSO's extreme phases has also been investigated with four GCM simulations (Hoerling et al. 2001), which were all found to have a 500-mb height response to extreme warm tropical Pacific SST that was twice as strong as the response to extreme cold SST. The longitudinal phase of the GCM's teleconnections also shifted eastward during warm events as compared with cold events, though this displacement is smaller than that observed. A nonlinear identification of the atmospheric response to ENSO was also addressed by Hannachi (2001) using general circulation models.

Nonlinear methods are required to investigate the nonlinear behavior of the North American climate variability and its relation to ENSO. Standard multivariate statistical techniques such as principal component analysis (PCA) [also known as empirical orthogonal function (EOF) analysis], and canonical correlation analysis (CCA) are linear methods. Composite analysis does not assume linearity, but is restricted to the analysis of differences between specific phases of the SO or equatorial SST indices. Recently, neural networks (NN; Hsieh and Tang 1998) have been used for nonlinear PCA (NLPCA; Kramer 1991) and nonlinear CCA (NLCCA; Hsieh 2000). NLPCA has been applied to the Lorenz three-component chaotic system (Monahan 2000), tropical Pacific SST, and sea level pressure fields (Monahan 2001). NLPCA was recently used to represent the quasi-biennial oscillation (QBO) in the equatorial stratospheric wind (Hamilton and Hsieh 2002), and to explore the nonlinear characteristics of Canadian surface air temperature (Wu et al. 2002). NLCCA has been applied to study the nonlinear relation between the tropical Pacific sea level pressure (SLP) and SST fields (Hsieh 2001a), as well as between the wind stress and SST fields (Wu and Hsieh 2002). Also, NLCCA has been applied to forecasting the tropical Pacific sea surface temperatures in the Experimental Long-Lead Forecast Bulletin (more information available online at www.iges.org/ellfb).

In this paper, the NLPCA (Hsieh 2001b) and NLCCA (Hsieh 2001a) models will be applied to study the wintertime climate variability over North America and its relation to ENSO as simulated in an ensemble of six 47-yr simulations produced with the Canadian Centre for Climate Modelling and Analysis (CCCma) second-generation atmospheric general circulation model (AGCM2). This paper is organized as follows. In section 2, the data and the two methods (NLPCA and NLCCA) are briefly introduced. Nonlinear modes of the winter 500-mb geopotential height and surface air temperature that were extracted by the NLPCA are described in section 3. Correlated nonlinear modes of atmospheric and tropical Pacific SST variability that were extracted by the NLCCA are described in section 4. Section 5 presents a summary and discussion.

## 2. Methodology and data

### a. NLPCA

A variable **x**, which consists of *l* spatial stations and *n* observations in time, can be expressed in the form **x**(*t*) = [*x*_{1}(*t*), … , *x*_{l}(*t*)], where for each *i* (*i* = 1, … , *l*), *x*_{i}(*t*) (*t* = 1, … , *n*) is a time series of length *n.* PCA is used to find a scalar variable *u* and an associated vector **a**, with

so that

where 〈⋯〉 denotes a sample or time mean. Here *u,* called the first principal component (PC), is a time series resulting from a linear combination of the original variables *x*_{i}, while **a**, the first eigenvector of the data covariance matrix (the first EOF), often describes a spatial pattern. The second PC can similarly be extracted from the residual **x** − **a***u,* and so on for the higher modes. In practice, the common algorithms for PCA extract all modes simultaneously by calculating the eigenvalues and eigenvectors of the data covariance matrix.

The fundamental difference between NLPCA and PCA is that NLPCA allows a nonlinear continuous mapping from **x** to *u* whereas PCA only allows a linear mapping. NLPCA is performed with a NN, such as that displayed in Fig. 1, which contains three “hidden” layers of variables (or “neurons”) between the input and output layers. These layers are called the encoding, bottleneck, and decoding layers, respectively. Four transfer functions *f*_{1}, *f*_{2}, *f*_{3}, *f*_{4} are successively used to map from the input layer to the output layer (**x** → **h**^{(x)} → *u* → **h**^{(u)} → **x**′), where **x**′ is the least square approximation of **x**. The mappings are

and **x** is the input column vector of length *l.* The encoding layer **h**^{(x)} is described by a column vector of length *m* (*m* is the number of the hidden neurons in the encoding layer). The parameters that control the transformation to this layer are 𝗪^{(x)}, which is an *m* × *l* weight matrix, and **b**^{(x)}, a column vector of length *m* containing the bias parameters. Index *k* ∈ [1, *m*]. The bottleneck layer contains a single neuron, which represents the nonlinear principal component *u.* The transformation between encoding and bottleneck layers is controlled by an *m*-element weight vector **w**^{(x)} and a bias parameter *b*^{(x)}. The decoding layer contains the same number of neurons *m* as the encoding layer, and the output layer is again a column vector of length *l.* The transformations between these layers are controlled by *m*-element weight and bias vectors **w**^{(u)} and **b**^{(u)} (decoding layer), and an *l* × *m* weight matrix 𝗪^{(u)} and an *l*-element bias vector **b**^{(u)} (output layer). Transfer functions *f*_{1} and *f*_{3} are generally nonlinear (e.g., the hyperbolic tangent function), while *f*_{2} and *f*_{4} are taken to be the identity function. See Hsieh (2001b) for more details.

Optimal values of the weight parameters for each mapping layer are found by minimizing the cost function *J* [see Eq (7)], rendering the outputs to be as close to the inputs as possible within constraints imposed by a weight penalty term (discussed below). Data compression is achieved at the bottleneck, with the bottleneck neuron representing a single degree of freedom, namely *u,* the nonlinear principal component (NLPC).

The nonlinear optimization was carried out with a quasi-Newton method. To avoid the local minima problem (Hsieh and Tang 1998, p. 1859), an ensemble of 30 NNs with random initial weights and bias parameters was run. Also, 20% of the data was randomly selected as testing data and withheld from the training of the NNs. For the NLPCA, runs where the mean-square error (mse) for the testing dataset was 10% larger than that for the training dataset were rejected to avoid overfitted solutions. Then the NN with the smallest mse was selected as the desired solution.

As noted above, the cost function used to identify the NLPCA model in this study has an extra weight penalty term,

where *p* ≥ 0 is called the weight penalty parameter. Increasing *p* increases the concavity of the cost function, thereby pushing the weights 𝗪^{(x)} to be smaller in magnitude and consequently yielding smoother and less nonlinear solutions than when *p* is small or zero. With a large enough *p,* the danger of overfitting is greatly reduced (Hsieh 2001b). The NLPCA was run repeatedly with various values of the penalty parameter (ranging from 0.001 to 0.02), and the solution with the smallest mse was chosen as the final solution.

After the first NLPCA mode has been subtracted from the data, the residual is again input into the NLPCA network to extract the second NLPCA mode. Because of the noisier conditions, the penalty parameter is increased to the range 0.01–0.2 for the second mode.

The weight penalty parameters in NLPCA (also in the NLCCA described below) are selected by means of a search. Future research will hopefully provide a more objective way to select these parameters.

### b. NLCCA

Given two sets of variables **x** and **y**, CCA is used to extract the correlated modes between **x** and **y** by looking for linear combinations

where the canonical variates *u* and *υ* have maximum correlation; that is, the coefficient vectors **a** and **b** are chosen such that the Pearson correlation coefficient between *u* and *υ* is maximized.

In NLCCA, we follow the same procedure as in CCA, except that the linear mappings in Eq. (8) are replaced by nonlinear mapping functions using two-layer feed-forward NNs. The mappings from **x** to *u* and **y** to *υ* are represented by the double-barreled NN on the left-hand side of Fig. 2. By minimizing the cost function *J* = −corr(*u, **υ*), one finds the parameters that maximize the correlation corr(*u, **υ*). After the forward mapping with the double-barreled NN has been solved, inverse mappings from the canonical variates *u* and *υ* to the original variables, as represented by the two standard feed-forward NNs on the right side of Fig. 2, are to be solved, where the cost function *J*_{1} is the mse of the output **x**′ relative to **x** (mse_{x}), and the cost function *J*_{2}, the mse of the output **y**′ relative to **y** (mse_{y}) are separately minimized to find the optimal parameters for these two NNs (see Hsieh 2001a for details).

An ensemble approach is also used for the NLCCA—runs where −corr(*u, **υ*), the mse_{x}, or the mse_{y} for the testing dataset were 10% larger than those for the training dataset were rejected. The NNs with the highest cor (*u, **υ*), and smallest mse_{x} and mse_{y} were selected as the desired solution. A weight penalty was again used (Hsieh 2001a) as in NLPCA.

In brief, the major advantage of NLPCA and NLCCA is that both are free from linear constraint. This allows us to identify principal components or canonical variates that vary along empirically derived curves through the data space instead of only straight lines. Consequently, a larger fraction of the variance can be explained (relative to PCA) or higher canonical correlation can be reached (relative to CCA), thus raising the possibility of improving the interpretation and forecasting climate using the NLCCA model. The main difficulty in the application of NLPCA or NLCCA is that the cost function often has multiple local minima (Hsieh and Tang 1998). It takes considerable computation to obtain a robust solution that approximates the “global” minimum by running the optimization many times from random initial parameters. Another drawback is overfitting, that is, fitting to the noise in the data. Using a weight penalty and reserving part of the data as validation data can alleviate overfitting.

### c. Data

The monthly mean 500-mb geopotential height (Z500) and surface air temperature (SAT) data used in this study were produced by the CCCma AGCM2, a spectral model with T32 resolution in the horizontal, 10 levels in the vertical, semi-implicit time stepping, and a full physics package (McFarlane et al. 1992; Boer et al. 1992). An ensemble of six 47-yr runs of AGCM2 were carried out, in which each integration was started from different initial conditions and forced by SSTs from the Global Sea Ice and Sea Surface Temperature (GISST) dataset (version 2.2; Rayner et al. 1996). Simulations were performed from January 1948 to November 1994, and were initially reported by Zwiers et al. (2000). For each run, anomalies were calculated by subtracting the monthly climatology based on the whole period. Monthly anomalies were then smoothed by taking a 3-month running mean and removing linear trends. Only the smoothed data for December–February (DJF) are analyzed, thus the total number of months used from the six AGCM2 runs is 840 [(47 × 3 − 1) × 6]. The domain of interest is 20°–76°N, 150°E–50°W covering the North Pacific and North America. For the SAT, only the data over land grids were used.

Monthly SST used in this study was from the reconstructed global historical SST datasets by Smith et al. (1996) for the period 1950–2000 with a resolution of 2° × 2°. As the AGCM2 runs end in 1994, the SST data was actually used up to November 1994 in the NLCCA {the total number of months used in NLCCA is 804 [(45 × 3 − 1) × 6]}. Similar data processing was done for the SST data, as for the Z500 and SAT data. The area of interest for the SST is restricted to the tropical Pacific (21°S–21°N, 123°E–71°W).

Prior to the NLPCA and NLCCA, ordinary PCA (i.e., EOF) analysis was conducted on the Z500, SAT, and the SST anomalies to compress the data into manageable dimensions and to insure that the estimated variance–covariance matrices that enter into CCA calculation can be inverted (Barnett and Preisendorfer 1987). EOFs of the three leading modes of Z500 and SAT anomalies are shown in Fig. 3. Variance contributions from these three modes of the Z500 anomalies are 32.8%, 17.8%, and 12.7%, respectively, and for the SAT anomalies, 26.8%, 16.1%, and 8.9%, respectively. The SAT anomaly pattern for each mode generally reflects the anomalous circulation implied by the corresponding Z500 mode, although the signs for EOF2 are opposite (Figs. 3b and 3e). The accumulated variance contribution of the five leading modes is 77.0% for the Z500 anomalies, 61.9% for the SAT anomalies, and 83.7% for the tropical Pacific SST anomalies.

The five leading PCs (i.e., the EOF time series) of the winter Z500 and SAT anomalies from 1948 to 1994 were used as the inputs to the NLPCA. For the NLCCA, only data after 1950 were used since the SST data were available from January 1950. Also for the NLCCA, SST anomalies (DJF) were repeated 6 times to pair with the atmospheric data (Z500 and SAT) from the six GCM ensemble runs. The five leading SST PCs, and the five leading Z500 (or SAT) PCs were used as the inputs to the NLCCA.

## 3. Nonlinear modes extracted by NLPCA

### a. Mode 1

Figures 4a,b, which show the first NLPCA mode for the Z500 and SAT anomalies, respectively, reveal nonlinear structures in both datasets. The Z500 NLPCA mode 1 explains 35.7% of the total variance, versus 32.8% explained by PCA mode 1 (straight line in Fig. 4a). For SAT, the NLPCA mode 1 explains 30.6% of the variance, versus 26.8% by the linear mode. A measure of the degree of nonlinearity is the ratio (*r*) between the MSE of the NLPCA mode 1 and that of the corresponding PCA mode. Smaller *r* means stronger nonlinearity. When *r* = 1, the nonlinear mode is reduced to the linear mode. Here *r* is 0.921 for the Z500, and 0.901 for the SAT, suggesting moderate nonlinearity in both datasets, with SAT being somewhat more nonlinear.

Changing the value of a PC (i.e., selecting a point on the straight line in Fig. 4a) has the effect of changing the amplitude but not the spatial structure of the corresponding EOF pattern (Fig. 3a). In contrast, both the structure and amplitude of the spatial pattern of the NLPCA mode change smoothly as the NLPC *u* (which traces the curve in Fig. 4a) changes value. The NLPCA maps the bottleneck neuron *u* to the output layer of the NN used to represent the data. This produces values for the first five PCs that, in turn, can be combined with the corresponding PCA spatial patterns (i.e., the EOFs) to yield the spatial pattern corresponding to *u.* When *u* takes on its minimum value, three large anomalies appear in the Z500 anomaly field over the North Pacific and the North American continent, resembling a negative PNA pattern (Fig. 5a), with negative height anomalies over Canada and the United States except western Alaska and the southeastern United States, where there are positive height anomalies. The SAT anomaly pattern associated with the minimum *u* is roughly in agreement with the Z500 anomaly pattern, with positive anomalies (1°–2°C) over the southeastern United States and negative anomalies over the rest of North America (Fig. 5c).

When *u* takes on its maximum value, the Z500 anomalies (Fig. 5b) show no closed contours north of 50°N, with the United States covered by negative height anomalies, and Canada by positive anomalies, implying a weakened jet stream because the geopotential gradient is reduced in the region that is usually occupied by the jet stream. With maximum *u,* the SAT anomalies (Fig. 5d) show that the United States and southern Canada are cooler, while other areas of Canada and Alaska are warmer.

The NLPCA Z500 anomaly pattern for maximum *u* (Fig. 5b) resembles EOF1 pattern (Fig. 3a) while the pattern for minimum *u* (Fig. 5a) resembles EOF2 (Fig. 3b). This means that the weakened jet stream state associated with maximum *u* does not have an equally strong negative counterpart (i.e., an enhanced jet stream state). Instead the enhanced jet steam state is overshadowed by the negative PNA state found during minimum *u.*

Comparing Figs. 5a with 5c, and 5b with 5d, we can see some spatial correspondence between the anomalies of Z500 and SAT, where positive anomalies of Z500 tend to occur roughly together with positive anomalies of SAT, and similarly for the negative anomalies, thereby revealing some consistency in the nonlinear structures found by the NLPCA approach.

### b. Mode 2

After removing the NLPCA mode 1 from the data, the residuals were again input into the NLPCA network to extract the second NLPCA mode, which was also found to contain notable nonlinearity in both the Z500 and SAT fields. The second Z500 NLPCA mode explains 16.5% of the total variance, which is slightly more than the 15.2% explained by the corresponding linear mode—here the linear mode is not the same as the PCA mode 2, as the latter is extracted from the residual with the PCA mode 1 (not the NLPCA mode 1) subtracted from the original data. The mse ratio between the second NLPCA mode and the corresponding linear mode for Z500 is 0.953. There is stronger nonlinearity in the second SAT NLPCA mode, which explains 12.5% of the total variance, versus the 9.3% explained by the corresponding linear mode, with an mse ratio of 0.874.

The second Z500 NLPCA mode reveals mainly east–west differences over North America and the adjacent North Pacific. Positive anomalies lie in the west and negative anomalies to the east at minimum *u* (Fig. 6a). The anomalies reverse sign and shift southwestward for maximum *u* (Fig. 6b). The second SAT NLPCA mode shows even larger changes in the anomaly pattern when comparing minimum *u* (Fig. 6c) and maximum *u* (Fig. 6d).

The Z500 and SAT anomaly patterns corresponding to minimum *u* (Figs. 6a and 6c) resemble EOF3 (Figs. 3c and 3f), with positive height and SAT anomalies centered over western Canada. Spatial patterns corresponding to maximum *u* (Figs. 6b and 6d) resemble EOF2 for both Z500 and SAT (Figs. 3b and 3e) with positive height and SAT anomalies over much of the United States and Canada. We can see some similarity between the spatial structure of the second Z500 and SAT NLPCA modes, even though NLPCA was performed separately on the two datasets.

To investigate the relations between the NLPCA modes and ENSO, composites of the winter (DJF) Z500 and SAT anomalies during warm (El Niño) event years (1958, 1966, 1969, 1973, 1983, 1987, 1988, and 1992) and cold (La Niña) event years (1950, 1951, 1955, 1956, 1965, 1971, 1974, 1976, and 1989) were computed. The years used for the composites are the same as those used by Hoerling et al. (1997). The Z500 composites during El Niño years and La Niña years show positive and negative PNA patterns, respectively (Figs. 7b and 7a). When El Niño takes place, positive height anomalies are dominant over North America except over the southeastern United States. The corresponding composite SAT anomaly field has negative SAT anomalies over the southeastern United States and positive SAT anomalies over the rest of North America (Fig. 7d). The composite Z500 and SAT anomalies during La Niña years are basically opposite to those during El Niño years. Thus, the asymmetries between El Niño and La Niña in the midlatitudes simulated by the CCCma AGCM are weaker than observed (Hoerling et al. 1997). Patterns displayed in Figs. 7a and 7c are similar to those shown in Figs. 5a and 5c, but with weaker amplitudes. However, patterns in Figs. 7b and 7d are quite different from those shown in Figs. 5b and 5d. Apparently, NLPCA mode 1 is related to ENSO only when *u* is negative. To relate the tropical Pacific SST to the North American Z500 and SAT, we turn to the NLCCA approach.

## 4. Nonlinear modes extracted by NLCCA

### a. NLCCA of SST and Z500

NLCCA between winter tropical Pacific SST anomalies and simultaneous Z500 anomalies (Fig. 8) reveals nonlinearity when compared to the linear CCA solution, which is shown as a straight line. For SST (Fig. 8a), nonlinearity in the PC_{1}–PC_{2} plane is manifested by a curve that links La Niña states at the left end to El Niño states at the right end. For Z500 (Fig. 8b), stronger nonlinearity appears in the PC_{1}–PC_{3} and PC_{2}–PC_{3} planes than in the PC_{1}–PC_{2} plane. The first NLCCA mode for SST explains 64.6% of SST variance, versus 62.3% explained by the first CCA mode. The mse ratio *r* is 0.856. The corresponding NLCCA mode for Z500 explains 23.6% of the variance, versus 21.5% explained by the first CCA mode, with a *r* value of 0.953. The correlation between the canonical variates (*u* and *υ*) is 0.702 for the nonlinear mode and 0.675 for the linear mode.

As in NLPCA, one can map values of canonical variate *u* and *υ* onto SST and Z500 anomaly patterns, respectively. Here, minimum and maximum *u* are chosen to present the La Niña and El Niño states, respectively, and Z500 anomaly patterns are considered for the values of *υ* that correspond to minimum and maximum *u.* The SST field that corresponds to minimum *u* presents a La Niña with negative anomalies (about −2.0°C) over the central-western equatorial Pacific. The corresponding Z500 field has a negative PNA pattern with a positive anomaly center over the North Pacific, a negative center over western Canada and a positive center over the eastern United States (Fig. 9a).

When *u* takes on its maximum value, the SST field presents a fairly strong El Niño with positive anomalies (about 2.5°–3.0°C) over the central-eastern Pacific (Fig. 9b). The SST warming center shifts eastward by 30°–40° longitude relative to the cooling center in Fig. 9a. This asymmetric SST variation between El Niño and La Niña states has also been found by Hsieh (2001a) and Wu and Hsieh (2002). Note that the warming in Fig. 9b does not display a maximum off Peru, in contrast to the first NLPCA mode of SST (Fig. 10d of Hsieh 2001b). This difference between the first NLPCA and NLCCA modes suggests that warming confined to the eastern equatorial Pacific does not have a strong midlatitude atmospheric response, in agreement with Hamilton (1988). The Z500 response field contains a PNA pattern (Fig. 9b) that is roughly opposite to that shown in Fig. 9a, but with a notable eastward shift. The zero contour surrounding the North Pacific anomaly lies close to the western coastline of North America during El Niño (Fig. 9b), while it is about 10°–15° farther west during La Niña (Fig. 9a). The positive anomaly over eastern Canada and the United States in Fig. 9a becomes a negative anomaly shifted southeastward in Fig. 9b. Evidently, the Rossby wave train linking the extratropical atmospheric response to the tropical source changes between La Niña and El Niño events in the CCCma model. As the equatorial SST anomalies shift eastward from La Niña to El Niño, the Rossby wave train also shifts eastward. However, the atmosphere's eastward shift is not as large as that in the SST anomalies. This is probably because the atmospheric response is a direct result of tropical heating anomalies (the response to tropical convection) instead of the SST anomalies, although the former is caused by the latter. The amplitude of the Z500 anomaly over the North Pacific is stronger during El Niño than La Niña, but the anomaly over western Canada and the United States is weaker during El Niño than La Niña (Figs. 9a and 9b).

For comparison, CCA mode 1 describes a pattern of atmospheric response to La Niña that, apart from magnitude, is opposite to the response pattern for El Niño (Figs. 9c and 9d). Note, however, that the El Niño response does have somewhat stronger amplitudes. The SST anomaly patterns extracted with CCA are also completely symmetric between the two extremes. The anomaly centers are located at about the average positions between those shown in Figs. 9a and 9b.

### b. NLCCA of SST and SAT

NLCCA was also applied to analyze the covariability of SST and SAT. The first NLCCA mode for SST (shown in Fig. 10a) is very similar to that shown in Fig. 8a. For SAT, considerable nonlinearity occurs between the PC_{1} and PC_{2} (Fig. 10b). The first NLCCA mode for SST explains 65.5% of the total variance, versus 62.3% explained by the first CCA mode, with *r* being 0.850. The first NLCCA mode for SAT explains 23.2% of the total variance, versus 22.1% by the first CCA mode, with an mse ratio *r* of 0.967. The canonical correlation is 0.605 for NLCCA and 0.588 for CCA.

Figures 11a,b show the spatial anomaly patterns for both SST and SAT associated with La Niña and El Niño, respectively. When *u* takes on its minimum value, positive SAT anomalies (about 1°C) appear over the southeastern United States, while much of Canada and the northwestern United States are dominated by negative SAT anomalies. The maximum cooling center (−4°C) is located over northwestern Canada and Alaska (Fig. 11a). When *u* takes on its maximum value (Fig. 11b), the warming center (3°C) is shifted to the southeast of the cooling center in Fig. 11a, with warming over almost all of North America except the southeastern United States. The SAT anomaly patterns here are roughly consistent with the composite analysis from observed data by Hoerling et al. (1997). The first CCA mode for SAT and SST is shown for reference in Figs. 11c,d.

We note that the nonlinear SAT response is more strongly captured by NLCCA than by composite analysis (cf. Figs. 9a,b with 7a,b and Figs. 11a,b with 7c,d). This is not surprising because the averaging in the composite method mixes the strong and weak events, thereby producing weaker features than NLCCA. Another disadvantage with the composite approach is that one must decide, a priori, which events to include in the composites.

## 5. Summary and discussion

Two nonlinear multivariate statistical techniques, NLPCA and NLCCA, developed from neural networks, were applied to investigate the nonlinear behaviour of the North American winter climate as simulated by the CCCma AGCM2. NLPCA was used to find the nonlinear modes that could account for the maximum variance in a dataset, while NLCCA was used to detect the most strongly correlated nonlinear modes relating the atmospheric (Z500 or SAT) and tropical Pacific SST variability. These two methods provided valuable nonlinear diagnostics of a GCM simulation.

Nonlinearity is detected by NLPCA in both the winter Z500 and SAT anomaly fields, which display asymmetric spatial patterns of variability during the opposite extremes of the nonlinear principal component (NLPC) *u.* The first NLPCA mode for Z500 describes a negative PNA pattern during the minimum *u,* and a weakened jet stream pattern during maximum *u.* Although NLPCA was conducted on the Z500 and SAT anomalies separately, the derived NLPCA modes have similar spatial patterns.

The first NLCCA mode shows there is nonlinearity in both the SST anomaly field and the related atmospheric response fields (Z500 and SAT). The leading NLCCA mode for SST is a nonlinear ENSO mode, showing asymmetry between the warm El Niño states and the cool La Niña states with a 30°–40° longitude eastward shift of the SST anomalies during El Niño relative to the anomalies during La Niña. The leading NLCCA mode for Z500 is a nonlinear PNA teleconnection pattern, in which the positive PNA pattern associated with El Niño is situated about 10°–15° east of the negative PNA pattern associated with La Niña. As with NLPCA, the spatial anomaly pattern of the first NLCCA mode for the SAT generally corroborates with the corresponding Z500 anomaly pattern. The NLCCA was able to extract a stronger nonlinear atmospheric response to ENSO than the composite method.

Even with NLCCA, the nonlinearity of the atmospheric responses to ENSO reported in this paper is weaker than that found by Hoerling et al. (1997), where an appreciable 30° longitude phase shift between the warm and cold event circulation composites was detected from observational data. However, we should note that, the atmospheric sample available here from an ensemble of six runs is much bigger than is available for an observational study, so the opportunities for overfitting and consequently finding apparent strong nonlinearity are somewhat reduced. It might therefore be more appropriate to compare our results with other models. Hoerling et al. (2001) surveyed the response to ENSO in four GCMs. Nonlinearity was found in all four models, existing in both the strength of the midlatitude response and its spatial phase. Similar to our results, the spatial phase displacements of the four GCMs' teleconnections during warm events and cold events are smaller than observed. The degree of nonlinear atmospheric responses to ENSO is model dependent.

Hoerling et al. (2001), used one-sided linear regression to analyze the nonlinearity. The procedure involves calculating linear regressions between all occurrences of one sign of an SST index (e.g., the leading EOF PC of tropical Pacific SST anomalies) and the response fields (e.g., Z500 or SAT anomalies). For comparison, this method was also applied to CCCma AGCM2's simulation and the results are shown in Fig. 12. The regression patterns for Z500 and SAT anomalies are generally consistent with the spatial patterns given by NLCCA (Figs. 9a,b and 11a,b), with the positive PNA teleconnection shifted 10°–15° eastward of the negative PNA teleconnection (Figs. 12a,b). The phase difference in the CCCma AGCM2, as indicated by the placement of zero contours, is comparable to that in version 3 of the National Center for Atmospheric Research (NCAR) Community Climate Model (CCM3) and the ECHAM-3 models, but weaker than that in the Geophysical Fluid Dynamics Laboratory (GFDL) and Medium-Range Forecast (MRF) models [see Hoerling et al. (2001) for details of the four models]. The negative Z500 anomalies over the North Pacific are much stronger during warm events than the positive anomalies that occur during cold events. However, the anomalies over the North American continent have similar amplitudes during warm and cold events (Figs. 12a,b), which can also be seen in the NLCCA results (Figs. 9a,b).

Hence the CCCma AGCM2 is basically capable of simulating the nonlinear responses of North American climate to ENSO, although the nonlinearity is somewhat weaker than that observed. NLCCA, which successfully extracted the nonlinear mode, offers a fully nonlinear diagnostic tool to study GCM output.

## Acknowledgments

This work was supported by the Natural Sciences and Engineering Research Council of Canada via a strategic grant to Hsieh, Shabbar, and Zwiers, and a research grant to Hsieh. This paper was improved by constructive comments on an earlier draft provided by Greg Flato and John Fyfe.

## REFERENCES

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

## Footnotes

*Corresponding author address:* Dr. Aiming Wu, Dept. of Earth and Ocean Sciences, University of British Columbia, Vancouver BC V6T 1Z4, Canada. Email: awu@eos.ubc.ca