Observed Feedback between Winter Sea Ice and the North Atlantic Oscillation

Feedback between the North Atlantic Oscillation (NAO) and winter sea ice variability is detected and quantiﬁed using approximately 30 years of observations, a vector autoregressive model (VAR), and testable deﬁnitions of Granger causality and feedback. Sea ice variability is quantiﬁed based on the leading empirical orthogonal function of sea ice concentration over the North Atlantic [the Greenland Sea ice dipole (GSD)], which, in its positive polarity, has anomalously high sea ice concentrations in the Labrador Sea region to the southwest of Greenland and low sea ice concentrations in the Barents Sea region to the northeast of Greenland. In weekly data for December through April, the VAR indicates that NAO index ( N ) anomalies cause like-signed anomalies of the standardized GSD index ( G ), and that G anomalies in turn cause oppositely signed anomalies of N . This negative feedback process operates explicitly on lags of up to four weeks in the VAR but can generate more persistent effects because of the autocorrelation of G . Synthetic data are generated with the VAR to quantify the effects of feedback following realistic local maxima of N and G , and alsoforsustainedhighvaluesof G . Feedbackcanchangetheexpectedvalueofevolvingsystemvariablesbyas much as a half a standard deviation, and the relevance of these results to intraseasonal and interannual NAO and sea ice variability is discussed.


Introduction
The North Atlantic Oscillation (NAO) strongly influences weather and climate over the North Atlantic and surrounding continents, and its forcing mechanisms and predictability have received considerable attention (Marshall et al. 2001;Trigo et al. 2002;Hurrell et al. 2003). Although the NAO is generally considered to be an internal mode of atmospheric variability, it has important interactions with other components of the climate system. The NAO strongly impacts North Atlantic sea surface temperatures (SSTs), and a weak but nonnegligible influence of tropical and extratropical SST on the NAO has been found (Marshall et al. 2001;Kushnir et al. 2002;Czaja et al. 2003;Mosedale et al. 2006).
There are also important interactions between the NAO and winter sea ice on synoptic to decadal time scales. As the leading mode of variability, the NAO might be expected to exert an influence on the spatial distribution of winter sea ice via wind-driven anomalies of sea ice velocity, surface vertical heat flux, and possibly horizontal oceanic heat flux. Observations show that positive NAO years tend to have anomalously high sea ice concentrations over the Labrador Sea and anomalously low sea ice concentrations over the Barents Sea, which is consistent with the NAO's circulation pattern (Deser 2000). Accordingly, there is strong observational evidence connecting Arctic sea ice decline with the increasing NAO trend from the 1960s to the early 1990s (e.g., Rothrock and Zhang 2005), although a signal of overall Arctic sea ice decline is emerging that is not directly attributable to a trend in the overlying atmospheric circulation (Deser and Teng 2008).
Modeling studies indicate that NAO-driven sea ice variability in turn feeds back onto the NAO Deser et al. 2004;Alexander et al. 2004;Kvamstø et al. 2004;Seierstad and Bader 2009). This is a negative feedback, meaning that the sea ice patterns associated with the positive polarity of the NAO tend to generate negative NAO-like atmospheric response patterns. Deser et al. (2007) recently examined the transient response of an atmospheric global climate model forced by a positive NAO-driven sea ice anomaly pattern. Their results revealed that the negative feedback process begins with a localized baroclinic response that reaches peak intensity in 5-10 days and persists for 2-3 weeks. If the ice forcing is maintained for several months, the atmosphere develops a larger-scale, equivalent barotropic response that resembles the negative polarity of the NAO, and this pattern is maintained primarily by nonlinear transient eddy fluxes of vorticity (Deser et al. 2007) related in part to changes in tropospheric Rossby wave breaking (Strong and Magnusdottir 2009b).
The research presented here was undertaken to test whether the negative feedback from sea ice to the NAO could be detected in observations and, if so, to quantify this feedback. We tested for the presence of negative feedback and investigated its effects within the framework of a linear stochastic model representing the evolution of the NAO index and sea ice concentration as a vector autoregressive process. Mosedale et al. (2006) recently used a similar statistical framework for examining the relationship between the NAO and Atlantic SSTs and found that the SST-to-NAO feedback is responsible for slower-than-exponential decay in lag autocorrelations of the NAO at lags longer than 10 days. We describe our methods in section 2 with additional details given in the appendix. The manuscript then proceeds with three results sections in which we illustrate the interaction between sea ice and the NAO, develop a linear stochastic model of the interaction, and perform numerical experiments with the model. The final section of the manuscript provides a summary and discussion.

a. Data
We use National Snow and Ice Data Center (NSIDC) sea ice concentrations derived from Nimbus-7 Scanning Multichannel Microwave Radiometer and Defense Meteorological Satellite Program Special Sensor Microwave Imager radiances (Cavalieri et al. 2008). These sea ice data are available on 25-km grids nominally once every two days for 1978-86 and once daily for 1987-2007. There is a gap in the data from 3 December 1987 to 13 January 1988, so we do not include winter 1987/88 in our study. For 2008, we use the NSIDC's preliminary version of the same sea ice data (Meier et al. 2008). We also use sea level pressure and surface 10-m wind ve-locity data from the National Centers for Environmental Prediction (NCEP)-National Center for Atmospheric Research (NCAR) Reanalysis for the same years.
Our study focuses on processes occurring at weekly time scales during winter to early spring when the NAO intensity and sea ice extent are at maximum. We thus analyze weekly mean data for the 21 weeks with starting dates ranging from 4 December through 23 April, where 4 December was chosen to give an integer number of weeks in December. Where we calculate lagged correlations for these 21 weeks, we incorporate weekly mean values for data prior to 4 December or after 23 April as needed. The data were detrended and, for each grid point and week, seasonality was removed by subtracting the 1978-2008 mean for that grid point and week. Data processed as described above are referred to as ''21week'' data (e.g., 21-week sea level pressure).

b. Indices
To define the indices used in the study, we use empirical orthogonal function (EOF) analysis. Specifically, the leading eigenvector of the spatial covariance matrix of the data is the first EOF, and the time series obtained by projection of the data onto that EOF is the leading principal component. Using 21-week data (defined in section 2a), we define the NAO index (N) as the leading principal component of sea level pressure for the period 4 December 1978 through 23 April 2008 over the domain used in Hurrell (1995): 208-808N and 908W-408E. A portion of the associated EOF is shown in Fig. 1a and features the well-known meridional seesaw with a positive (negative) center of action near 408N (658N). We define a sea ice index for the same time period as the leading principal component of 21-week sea ice concentration over the domain 408-908N and 908W-908E. The EOF associated with this sea ice index ( Fig. 1b) resembles that found in observations by Deser et al. (2000) and has a positive center of action to the southwest of Greenland and a negative center of action extending from the Greenland Sea to the Barents Sea. We refer to this sea ice variability pattern as the Greenland Sea ice dipole (GSD) and denote its principal component by G. As an example of the data, the N and G indices for 2005-08 are shown in Fig. 2.

c. Linear stochastic model
We develop a linear stochastic model of the relationship between N and G within the framework of a qth-order vector autoregression [VAR(q)]. The VAR has two equations, one for the value of G at week t, denoted G t , and one for the value of N at week t denoted N t . For the so-called primitive or structural form of the VAR, we write where G t and N t are assumed to be stationary, « Gt and « Nt are white noise disturbances with respective standard deviations s G and s N , and the associated disturbance time series « G and « N are uncorrelated. The f terms allow for the possibility of a contemporaneous influence of one variable on another, and the coefficients inside the summations provide a means for lagged relationships.
The system cannot be estimated directly unless one of the f terms is restricted to zero (e.g., Sims 1980). We choose the restriction f 2 5 0, meaning that the sea ice does not contemporaneously affect the NAO. This is a reasonable assumption since transient boundary forcing experiments show that it takes several weeks for the NAO-driven sea ice variability to feed back onto the NAO (Deser et al. 2007). At the same time, it is reasonable to allow f 1 to take on a nonzero value because the integrated effect of one week of NAO-related atmospheric variability can be expected to impact the spatial distribution of sea ice during that same week via, for example, wind-driven sea ice velocity anomalies. The VAR in (2.1) provides a means for exploring the relationship between N and G and invoking the testable definitions of causality and feedback pioneered by Granger (1969). More details of the methods associated with the VAR are provided as the model is fit and synthetic data are generated with the VAR in section 3b and the appendix.

d. Comparing models
We will at times compare the relative predictive skill or strength of two versions of a statistical model: un-restricted and restricted. The second model is a restricted version of the first in that one predictor has been eliminated by restricting its regression coefficient to zero or, in the case of an autoregressive model, the order of the restricted model is reduced by one. To compare the models, we used the log-likelihood ratio given by Sims (1980): where T is the number of usable observations, c is the maximum number of regressors in the longest equation, and jS u j and jS r j are the determinants of the covariance matrices of the unrestricted and restricted model's residuals, respectively. The test statistic L has an asymptotic x 2 distribution with degrees of freedom equal to the number of restrictions in the system; it was tested for significance at the 95% confidence level. We use L to determine the appropriate lag order for autoregressive models and also to compare models for the purpose of testing for Granger causality.

Results
In the first results section (3a), we use correlations and composite time series to illustrate the relationship between the observed NAO and sea ice variability. Then, in section 3b, we use the VAR to model the observed relationships and test for causality and feedback. Section 3c presents the results of some experiments performed with the VAR.

a. Illustrating the relationship between sea ice and the NAO
The spatial patterns of the EOFs associated with N and G (Fig. 1) suggest a physical relationship between the two underlying phenomena. The circulation anomalies around the NAO's centers of action can induce sea ice concentration anomalies via wind-driven sea ice advection anomalies as discussed in, for example, Deser et al. (2000). Analysis of coupled model output further indicates that NAO circulation anomalies induce sea ice melting or freezing anomalies in the Labrador and Barents Seas via wind-driven oceanic heat transport (Koenigk et al. 2009;Strong and Magnusdottir 2009a). When the NAO is in its positive polarity, sea ice concentrations tend to be anomalously high to the southwest of Greenland and anomalously low over the region extending from the Greenland Sea northeast into the Barents Sea.
The principal components N and G are autocorrelated and cross correlated on a variety of time scales and lags. The autocorrelation functions in Fig. 3a show that G has a longer time scale than N. The partial autocorrelation functions in Fig. 3b indicate the amount of autocorrelation not accounted for by lower lags and suggest that a model of this system might include lags of at least three weeks. The contemporaneous correlation between G and N is positive and accounts for 6% of the variance (r 5 0.24 at lag zero; Fig. 4a). This positive r(N, G) at zero lag motivates the inclusion of a f 1 term in the VAR and reflects the ability of NAO-driven circulations to generate sea ice advection anomalies on the time scale of days to a week. When N leads G by one or several weeks (left side of Fig. 4a), the correlation r(N, G) strengthens over the first week and remains positive, which is consistent with the concept (e.g., Zhang et al. 2000) that thermodynamics are essential to the connection between sea ice and the NAO. More specifically, NAO-driven circulation anomalies may need to integrate over several weeks in order to advect sufficient heat in the ocean and atmosphere to generate the melting anomaly patterns associated with the leading sea ice EOF.
When G leads N (right side of Fig. 4a), correlations rapidly drop toward zero and then become increasingly negative after three weeks. The correlation pattern when G leads N may be observational evidence of the negative feedback from sea ice to the NAO posited in previous studies-a hypothesis we will test in section 3b. When G and N are averaged over the 21 weeks making up each extended winter (Fig. 4b), the overall correlation is positive, which is consistent with the dominance of positive correlation on within-season time scales seen in Fig. 4a.
We now use a set of composite time series (Fig. 5) to further explore the lagged relationships between N and G. Lag zero in each composite corresponds to the time of a locally extreme value of either N or G as indicated in the caption for each panel in Fig. 5. We selected 29 distinct cases for each composite to have an average of one case per year. These composite time series provide a complementary view of the lagged relationships suggested by Fig. 4a.
The composite results indicate that locally extreme values of N tend to be followed by like-signed locally extreme values of G. Based on the shaded interquartile range in Fig. 5, more than 75% of the G values were positive two weeks after the local maximum of N (Fig. 5a), and approximately 75% of the G values were negative one week after the local minimum of N (Fig. 5c). Examining the behavior of N about locally extreme values of G suggests the presence of a negative feedback. Just after G reaches a local maximum, the composite N time series undergoes a strikingly rapid decrease, and approximately 75% of the N values in the composite are negative one week after the composite G maximum (Fig. 5b). The composite local minimum of G is followed by a composite local maximum of N, and approximately 75% of the N values in the composite are positive two weeks after the G minimum (Fig. 5d). in the relationship between observed sea ice and the NAO. In this section, we formally test these hypothesized interactions by fitting the VAR to our data and invoking the testable definitions of causality and feedback introduced by Granger (1969). Specifically, one variable X is said to Granger cause another variable Y, denoted X 0 Y, if past values of X contain information that improves the prediction of Y above and beyond the information contained in past values of Y alone. If X Granger causes Y and also Y Granger causes X, there is a Granger feedback between X and Y, denoted X * 0 Y.
To determine if Granger causality is present in a system, we compare the strength of a null model with Granger feedback to the strength of a model with one direction of causality eliminated. The relative strength of models is assessed using a log-likelihood ratio test of the model residuals at the 95% confidence level as described in section 2d. If the model with feedback is stronger than the null model with a direction of causality eliminated, we conclude that Granger causality is present in that direction.
For our application to N and G, we defined three models:  .1), and the process of transforming this system for the purpose of data fitting is described in the appendix. The other two models are developed from M[G * 0 N] by restricting certain parameters to zero-a process also described in the appendix. As an example, M[G * L N] is the system (2.1) with the additional restriction g i 5 0 for i 5 1, . . . , q, meaning that N evolves unaffected by past values of G. To determine the model order, we used the log-likelihood ratio test described in section 2d to incrementally compare M[G * 0 N] at some arbitrarily high-order q to itself at order q 2 1. The model order was the first q for which the q 2 1 model had a significantly larger residual error. Beginning this process at q 5 7, the parsimonious model order was found to be q 5 4.
The correlations in Fig. 4 and composites in Fig Table 1 and its caption. The first-order influence of N on G is positive (b 1 ) and the first-order influence of G on N is negative (g 1 ), consistent with the correlations in Fig. 4 and the composite results in Fig. 5. The white noise residual term to N (« N ) has a standard deviation of s N 5 0.95, which is almost 3 times as large as the corresponding noise term to G (« G with standard deviation s G 5 0.32). This reflects the short time scale of N relative to G. The negative sign on f 1 is consistent with the positive contemporaneous correlation r(N, G) in Fig. 4a. Here M[G * 0 N] accounts for 32% of N and 92% of G, whereas M[G * L N] accounts for 30% of N and 90% of G. These somewhat small differences in performance are discernible as significant in part because of the large sample size. For longer averaging periods, the effect of shorter-term fluctuations becomes less pronounced and the influence of feedback becomes a more important contributor to the total variance as discussed in Mosedale et al. (2006).
The presence of feedback in the system also changes the spread of the N and G distributions. To illustrate this, we generated synthetic, perpetual wintertime series N and G of length one million weeks using M[G * 0 N]. We then generated control time series N* and G* from M[G * L N]. These control time series represent the behavior of the system with feedback turned off via elimination of the Granger causality G 0 N (i.e., g i 5 0 for i 5 1, . . . , 4). Turning off the negative feedback from sea ice to the NAO increased the variance of N by a factor of 1.08 and increased the variance of G by a factor of 1.60. These are both statistically significant variance increases under bootstrapping (e.g., Boos and Brownie 1989), meaning that sea ice and NAO-driven weather and climate over the North Atlantic are significantly dampened by the presence of this negative feedback.

c. Model experiments
To further explore the effect of feedback in the sea ice-NAO system, we performed two additional sets of experiments with the VAR. The first set, in section 3c(1), illustrates how feedback affects the evolution of N and G following a realistic locally extreme value of G. The second set, in section 3c(2), illustrates how feedback affects the evolution of N and G following a sustained high or sustained low value of G.

1) EFFECT OF LOCAL MAXIMUM OF G
For these experiments, we used as initial conditions the composite values from 23 # t # 0 in Fig. 5b and then advanced M[G * 0 N] and M[G * L N] forward for 20 weeks (Fig. 6a). Consistent with the convention in the previous section, we denote the output from M[G * 0 N] by N and G and denote the output from M[G * L N] by N* and G*. We set the white noise terms « N and « G to zero and produced single realizations, which is equivalent to developing large ensembles of time series and then averaging since the white noise has zero mean. Although the initial conditions were based on a composite around local maxima of G, the simulation results are also informative regarding local minima of G since the VAR output reflects across the time axis if all the initial conditions are reflected across the time axis (i.e., if the initial conditions are scaled by 21).
In the control model (M[G * L N]), the initial N* and G* anomalies decayed monotonically toward zero (dashed curve and curve with open circles in Fig. 6a). In the simulation with feedback (M[G * 0 N]), the decay of N was accelerated relative to the control because of the negative feedback from the initially high G values; N became negative after t 5 3 and then decayed toward zero after t 5 10 (curve with filled circles in Fig. 6a). The lower N values in turn forced a faster decay of G toward zero, consistent with the first-order positive influence of N on G.
We stated above that reflection of all the initial conditions about the time axis reflects the results across the time axis. It is informative to reflect only the G initial conditions across the time axis (Fig. 6b) to illustrate what happens when N and G are initially of opposite sign. This can occur in observations despite the positive contemporaneous correlation between N and G because of the large autocorrelation of G and volatility of N. Reflecting the initial G conditions across the time axis of course had no effect on N* (curve with open circles, Fig. 6b). Affected by the negative G, N decayed toward zero monotonically and less rapidly than N*. These results for like-signed versus oppositely signed initial conditions illustrate two points that will be treated more generally in the remainder of this subsection. First, the largest absolute difference between N and N* was greater for like-signed initial conditions than oppositely signed initial conditions (cf. Figs. 6b and 6a). Second, large likesigned initial conditions tend to result in N changing sign prior to decaying toward zero as in Fig. 6a, whereas large oppositely signed initial conditions tend to result in a monotonic decay of N toward zero.
The shifts in G and N due to feedback scale approximately linearly with the values of the initial conditions. To illustrate this, we simplify the initial conditions by setting them equal to constant values G 0 and N 0 from t 5 23, . . . , 0 and then find D N , which is the signed difference N t 2 N t * that has the largest absolute value over the interval 1 # t # 20. In other words, D N is the largest positive or negative change in the NAO index attributable to feedback in the model. We likewise define D G for the signed difference G 2 G*. Figures 6c  and 6d show surfaces of D G and D N for the twodimensional space of initial conditions 23.0 # G # 3.0 and 23.0 # N # 3.0. These two surfaces are well approximated by the following planes through the origin (r 2 . 0.99 for each): From these equations and Figs. 6c,d, we see that the responses are strongest (D G 5 60.48 and D N 5 60.33) when the initial conditions are like-signed, consistent with the comparison of Figs. 6a and 6b made above. We also see that D G and D N are both more sensitive to G 0 than to N 0 .
The D G and D N results in Figs. 6c,d depart somewhat from planar along an axis approximately paralleling G 0 5 0, with the departures for D N being more pronounced. The abrupt changes in D N in Fig. 6d separate two sets of results: 1) results in which the evolving N monotonically approached zero and 2) results in which the evolving N underwent a sign change and then decayed toward zero. In the vicinity of the arrows, the initial conditions were like-signed, as in Fig. 6a, and N underwent a sign change before approaching zero. At a critical value of G 0 , the N results transitioned to monotonic decays as in Fig. 6b.

2) EFFECT OF SUSTAINED VALUES OF G
The second set of experiments with the VAR were motivated by previously published boundary forcing experiments in which the atmosphere was allowed to evolve over an approximately constant sea ice anomaly pattern representing an extreme value of G. To duplicate this experimental setup in the VAR, we used the same initial N and G conditions from Figs. 6a,b but then maintained the sea ice at a constant value G h equal to its initial condition at time t 5 0. To model this system with and without feedback, we wrote an equation for N t affected by G h and an equation for N t * not affected by G h : where parameters remain as given in Table 1. We calculated (3.3) and (3.4) with « Nt 5 0, which is equivalent to averaging over a large ensemble of realizations. The resulting N and N* time series are shown in Figs. 7a,b. The sustained positive and negative sea ice values increased the maximal absolute response in N by approximately 50% over the response to realistic local G extrema (see Figs. 7a,b versus Figs. 6a,b). For the simulations shown in Figs. 7a,b, the difference d t [ N t 2 N t * asymptotically approached 60.21, which is 20.14G h (open circles in Figs. 7c,d). Further examples would illustrate that the asymptotic approach of d t to 20.14G h as t / ' holds generally and is independent of the initial conditions of N. To show this analytically, we write the equation for d t , which is (3.3) minus (3.4): (3.5) where d tÀi 5 N tÀi À N tÀi * . Equation (3.5) is a fourthorder, deterministic, linear difference equation with a particular solution: which is 20.14G h for the parameter values in Table 1.
To obtain an expression for how d t asymptotically approaches d, we consider the general solution to (3.5), which is the sum of the particular solution (3.6) and four homogeneous solutions: (3.7) In (3.7), i is the imaginary unit, the numbers in parentheses are the roots of the characteristic polynomial of (3.5), and d 1, . . . , 4 are arbitrary coefficients determined by four initial conditions. Since the absolute value of each root is less than unity, (3.7) expresses an asymptotic approach of d t to its particular solution. Fitting (3.7) to the initial conditions in Figs. 7a and 7b yields the curves in, respectively, Figs. 7c and 7d. These analytical results match closely the differences N t 2 N t * obtained from the model output (open circles in Figs. 7c,d).

Conclusions
Using observations from 1978-2008, we calculated two weekly indices: an NAO index, N, and an index G representing the seesaw of winter sea ice concentration between areas to the northeast and southwest of Greenland. We found that N and G are correlated on a variety of time scales and lags and determined that a fourth-order vector autoregressive model was appropriate for capturing these relationships in weekly data. Based on the VAR(4) model, we concluded that N Granger causes G and also that G Granger causes N, meaning there is Granger feedback in the winter NAO-sea ice system, and found that N anomalies Granger cause like-signed anomalies of G, consistent with the observational analysis in Deser et al. (2000). Our finding that G anomalies Granger cause oppositely signed anomalies of N has not, to our knowledge, been shown in observations before. We quantified this negative feedback for the range of initial conditions jGj, jNj # 3.0 by comparing a model with feedback to a control model wherein Granger feedback was turned off. With feedback turned off, the variance of N was inflated by a factor of 1.06 and the variance of G was inflated by a factor of 1.60, meaning that the presence of negative feedback dampens the variability of sea ice concentration over the North Atlantic and, to some degree, NAO-driven weather and climate as well. In the following, the phrase ''feedback response'' will be used to refer to results from the model with feedback (M[G * 0 N]) minus results from the control model in which feedback was turned off (M[G * L N]). The largest feedback responses for N and G occurred when the initial conditions were likesigned; also, N and G were more sensitive to the initial value of G than to the initial value of N. Using local extrema as initial conditions, the largest feedback responses were 60.33 for N and 60.48 for G. When G was held constant at a value G h for the entire experiment rather than allowing G to evolve from the local maximum, the N response was approximately 50% larger for realistic values of G h . The N response to the sustained G h value asymptotically approached the value 20.14G h independent of the initial N value-a result we showed analytically.
These findings provide an observational perspective for interpreting boundary forcing experiments in which the atmosphere is exposed to sustained sea ice anomalies corresponding to large positive values of G. Sustained, positive G-like sea ice anomalies generate a localized baroclinic response that lasts several weeks, after which wavemean flow interaction maintains a pattern of large-scale height anomalies resembling the negative polarity of the NAO (Deser et al. 2007;Strong and Magnusdottir 2009b). Our observational analysis indicates that G-like sea ice anomalies tend to decay after a few weeks based on autocorrelation and composite time series. The lack of sustained G anomalies in observations may appear to limit the time frame during which G anomalies can generate large-scale, negative NAO-like height anomalies via wave-mean flow interaction. However, we showed 21-week winters such as 2007/08 during which G was consistently positive (Fig. 2) with a mean greater than 1.0 (Fig. 4b). We can think of such winters as naturally occurring boundary forcing events wherein the atmosphere was exposed to unusually large G for 21 weeks, resulting in a reduction in the variance and value of the NAO index.
Multi-month sustained sea ice anomaly patterns may be of increasing relevance over the next century. The G value for an entire extended winter or series of extended winters could be markedly above or below average as a result of projected changes in sea ice transport or melting. Even if the energy spectrum of the associated weekly values of G resembles the contemporary spectrum, the shift in the seasonal mean G could impact the central tendency of the N distribution via the negative feedback shown here. It would be difficult to detect such a process operating in observations since the time scale of the transient feedback process is within season, whereas the relevant lags for the associated VAR and Granger causality tests would be interannual and longer. Further experimentation with coupled climate models will thus be useful for investigating the interaction between NAO and winter sea ice as each evolves under projected climate change.

APPENDIX
To fit the VAR model to data (e.g., Enders 2004), we concisely write (2.1) as where y t [ [G t , N t ] T , F and A i are the 2 3 2 matrices on the left-and right-hand sides, respectively, of (2.1), and e t [ [« Gt , « Nt ] T . Solving for y t yields where B i 5 F 21 A i and e t 5 F 21 e t . The system (A2) is fit to data using seemingly unrelated regression (Zellner 1962), yielding values for B i and the residuals e t . To recover the 3 1 4q parameters of the original model (2.1) from the results of the fitting, we write out (A2) and solve the system A i 5 FB i , i 5 1, . . . , q, and (A4) where S « is the covariance matrix of « t and S e is the covariance matrix of e t .
To define the null model where only N Granger causes G (G * L N), we fit (A3) imposing the restriction g i 5 0. To define the null model where only G Granger causes N (G G 0 N) we impose the restriction b i 2 f 1 h i 5 0. Note that neither of these restrictions eliminates the possibility of a contemporaneous effect of N on G since f 1 is unrestricted. The complete set of restrictions used in the full and null models is given in Table 2.