## Abstract

A statistical method is developed to assess the full climate feedback of nonlocal climate feedbacks. The method is a multivariate generalization of the univariate equilibrium feedback assessment (EFA) method of Frankignoul et al. As a pilot study here, the generalized EFA (GEFA) is applied to the assessment of the feedback response of sea surface temperature (SST) on surface heat flux in a simple ocean–atmosphere model that includes atmospheric advection. It is shown that GEFA can capture major features of nonlocal climate feedback and sheds light on the dynamics of the atmospheric response, as long as the spatial resolution (or spatial degree of freedom) is not very high.

Given a sample size, sampling error tends to increase significantly with the spatial resolution of the data. As a result, useful estimates of the feedback can only be obtained at sufficiently low resolution. The sampling error is also found to increase significantly with the spatial scale of the atmospheric forcing and, in turn, the SST variability. This implies the potential difficulty in distinguishing the nonlocal feedbacks arising from small-scale SST variability. These deficiencies call for further improvements on the assessment methods for nonlocal climate feedbacks.

## 1. Introduction

One critical issue in climate dynamics is understanding the feedback response of the atmosphere to its lower boundary forcing over the ocean and land. This feedback response is usually difficult to assess, because of the overwhelming internal atmospheric variability that occurs independent of the boundary forcing. In particular, in the real world, with a single realization, climate feedback can be assessed only statistically. In comparison, climate feedbacks in a climate model can be assessed dynamically with ensemble experiments to suppress internal atmospheric variability.

Based on the stochastic climate theory of Frankignoul and Hasselmann (1977), Frankignoul et al. (1998) proposed a simple statistical method for the assessment of midlatitude oceanic feedback to the atmosphere. Their method takes advantage of the time-scale separation between the rapid (<week) atmosphere and slow (>month) sea surface temperature (SST). As such, at a slow climate time scale, the response of an atmospheric variable *x*(*t*) to the SST *y*(*t*) can be approximated as (Frankignoul 1985)

with *n*(*t*) being an internal atmospheric variability. Here, *by*(*t*) represents the quasi-equilibrium atmospheric response to SST with *b* as the feedback parameter. The word “feedback” here is used because the SST variability is always forced predominantly by the atmospheric variability in the first place, and therefore its impact on the atmosphere represents a feedback to the atmosphere. For this reason, the atmospheric response is also sometimes called feedback response or simply feedback. The feedback parameter *b* can be assessed from the SST-led ocean–atmosphere covariance, as opposed to the simultaneous covariance. A simultaneous covariance tends to mix the effects of the atmospheric forcing with the SST feedback because part of the SST variability is forced by, and therefore is correlated simultaneously with, internal atmospheric variability. However, SST variability is not correlated with atmospheric internal variability of sufficiently later times, because SST cannot be forced by internal atmospheric variability of later times. Therefore, we have

where 〈*p*, *q*〉 represents the covariance between *p* and *q*, and *τ* is a lag longer than the persistence time of the atmospheric internal variability *τ _{n}*.

^{1}With (1.2),

*b*can be estimated by first multiplying

*y*(

*t*−

*τ*) on both sides of (1.1) and then ensemble averaging as (Frankignoul et al. 1998)

This approach makes use of the SST-led covariance with the atmosphere in a quasi-equilibrium response state and therefore will be called the equilibrium feedback assessment (EFA). Although the scalar formula of EFA in (1.3) was originally used to assess the local feedback of SST on the overlying atmosphere (Frankignoul et al. 1998), it has more general applications. Indeed, it can be used to assess the nonlocal atmospheric feedback to any single SST index *y*(*t*), such as a regionally averaged SST (Liu and Wu 2004), or an SST EOF coefficient (Ferreira and Frankignoul 2005). This univariate EFA in (1.3) has been applied to the feedback of extratropical SST on atmospheric dynamic fields (Czaja and Frankignoul 2002; Liu and Wu 2004; Ferreira and Frankignoul 2005; Liu et al. 2007), and the feedback of vegetation on climate (Z. Liu et al. 2006; Notaro et al. 2006). For convenience, the lower boundary variable will hereafter simply be called the SST.

To further assess the full climate feedback between the entire SST field and an atmospheric field, including local and nonlocal feedbacks, one may formally extend the univariate EFA in (1.3) by performing the EFA assessment of the atmospheric variability on the SST variability at each point. However, as will be shown later, the interpretation of this assessment is problematic, because SST variability at different points could covary with each other and therefore the contribution of each individual SST point to the atmosphere becomes unclear. This motivated us to further develop a more general EFA that can assess the full climate feedback response. The univariate EFA will be generalized to a multivariate EFA in section 2. As a pilot study here, the generalized EFA and its sampling error are studied in a simple coupled ocean–atmosphere model in section 3. A summary and further discussion are given in section 4. The generalized EFA is shown to be able to identify some important features of nonlocal climate feedbacks. In the mean time, it is shown that sampling error tends to increase significantly with the resolution of the data and the spatial scale of the SST variability. As a result, given a finite sample size, useful assessment can only be made at a limited resolution. The potential challenges and problems associated with the assessment of nonlocal climate feedback are also discussed.

## 2. The nonlocal assessment

### a. The feedback response

Given the anomalous atmospheric and SST fields of *I* and *J* points, respectively, the atmosphere and SST at time *t* can be written in vectors **x**(*t*) and **y**(*t*), respectively, as

As a direct extension of the univariate response in (1.1), the atmospheric variability in region *i* is assumed to be forced by the entire SST field as

with *n _{i}*(

*t*) as the local atmospheric internal variability. The

*b*represents the impact of the

_{ij}*j*th SST component on the

*i*th atmosphere component that is

*independent of other SST variability*. The equilibrium atmospheric response in (2.1) at

*T*succeeding times can be rewritten in the vector form as

The matrix for the response coefficient,

is called the *feedback response matrix*, or simply, *the feedback matrix*. It effectively acts as a Green function for full feedback responses. Our objective is to assess this feedback matrix.

The feedback matrix can be estimated using a multivariate generalization of the scalar EFA in (1.3) as follows. We define the variable matrix as

where the subscript *t* is for the convenience of the representation of lagged covariance to be discussed next. The atmospheric response in (2.2) can be written as

Right multiplication of 𝗬^{T}_{t−τ} on (2.4) yields

where the lagged cross-covariance matrices are

with the superscript “T” standing for the transpose. Since SST variability cannot be forced by atmospheric variability of later times, for infinite samples, we have

This leads to a generalized EFA (GEFA) estimator for the feedback matrix as

In practice, we can take *τ _{n}* = 0, if the data are binned over a period (e.g., monthly) longer than the adjustment time of the atmospheric internal variability.

For a finite sample, the sampling error is usually not zero and varies with the lag [i.e., ɛ(*τ*) = 𝗖_{NY}(*τ*)𝗖^{−1}_{YY}(*τ*) ≠ 0] and therefore the estimator 𝗕(*τ*) also varies with the lag. The sampling error tends to increase with lag, because of the decreasing autocovariance 𝗖_{YY}(*τ*) and, in turn, increasing 𝗖_{YY}^{−1}(*τ*), with lag [as discussed for the scalar EFA by Z. Liu et al. (2006)]. Therefore, the first lag is usually preferred. An estimator much more stable than (2.8) with lag turns out to be an integral estimator

because the vanishing of sampling error in (2.7) is now replaced by a weaker condition of a diminishing integrated sampling error Σ^{τ}_{s=τn+1}𝗖_{NY}(*s*) = 0. Hereafter, unless otherwise specified, we only use the estimator at the first lag *τ* = *τ _{n}* + 1, with the lag subscription in (2.9) omitted. [For the first lag, (2.8) and (2.9) are identical.]

In a climate model, in principle, the matrix 𝗕 can also be obtained dynamically with ensemble experiments, independent of any statistical assessment. This is important because it provides an independent check against the statistical method, which is obtained under assumptions, such as linearity. For a specific *b _{ij}*, we can perform an ensemble of atmospheric model experiments in which the SST anomaly is prescribed to be nonzero only at the point

*j*as

*y*(

_{j}*t*).

^{2}With 〈

*n*(

_{i}*t*)〉 = 0, the ensemble mean of the atmospheric model response in (2.1) gives 〈

*x*(

_{i}*t*)〉 =

*b*(

_{ij}y_{j}*t*). Therefore, the nonlocal feedback response can be assessed dynamically as

It is now clear that the statistical method is also critical for model–data comparison through a combined dynamical–statistical assessment that is potentially powerful for the understanding of climate feedback (Liu and Wu 2004; Liu et al. 2007). First, the statistical assessment is applied to both the observation and the model to evaluate the model performance on the feedback. If the model feedback is reasonable, we can then perform further specific dynamic experiments in the model to validate the method in (2.8) itself, which is derived under assumptions, such as linearity. The dynamical experiments also directly shed light on the mechanism of the atmospheric response.

### b. The total feedback response

A useful feedback matrix can also be derived by directly applying the scalar EFA in (1.3) to each pair of the response (*x _{i}*) and forcing (

*y*) as

_{j}This matrix will be called the *total feedback response matrix*, or simply the *total matrix*, and the corresponding feedback will be called the *total feedback*, to be distinguished from the true feedback represented by the feedback matrix 𝗕. (As will become clear later, the word “total” here is reminiscent of the total derivative, in contrast to the partial derivative analogy of the feedback matrix.) The two matrices 𝗔 and 𝗕 can be shown^{3} to be related to each other through the *equivalence relation*:

Here, 𝗠 is a nondimensional matrix of lagged SST autocovariance, called the *SST mutual response matrix* or simply the *mutual matrix*:

By definition, the diagonal elements satisfy *m _{jj}* = 1. Although inaccurate, the name “mutual response” here is meant to remind one of the analogy to the feedback response among the SSTs themselves.

^{4}

The equivalence relation in (2.11) shows that the total matrix 𝗔, in general, differs from 𝗕. In some sense, *a _{ij}* can be thought as the

*i*th atmospheric response to the

*j*th SST,

*only if y*

_{j}is the sole forcing on x_{i}in the coupled climate system such thatIn a general coupled system, however, *x _{i}* is also affected by SSTs of other regions and therefore

*a*does not represent the true response of

_{ij}*x*to

_{i}*y*alone. Rather, since all the SSTs vary slowly, at short lags,

_{j}*m*can be approximated as the regression coefficient

_{kj}*r*as

_{kj}Therefore, (2.12) suggests that *a _{ij}* could be thought as the total feedback response on

*x*from all the SSTs that covary with the SST

_{i}*y*in the coupled system. This is in contrast to the feedback

_{j}*b*that identifies the SST forcing impact on

_{ij}*x*from the SST

_{i}*y*alone. In this sense, the total matrix can be thought of as the total derivatives of the feedback responses to an individual SST forcing, while the feedback matrix 𝗕 corresponds to the partial derivative.

_{j}It is also clear from the equivalence relation that if the SSTs are largely independent of each other, 𝗠 is close to the identity matrix 𝗠 = 𝗜, and therefore 𝗔 = 𝗕. Now, the scalar EFA and GEFA give the same results. This also highlights the point that the total matrix can be substantially different from the feedback matrix 𝗕 if the SST forcings are correlated.

It should be pointed out that, even if 𝗔 is not the same as 𝗕, 𝗔 still provides important information on the dynamics of the atmospheric response. This can be seen by rewriting the equivalence relation in (2.12) as *J* atmospheric response relations:

where

are the atmospheric response field and mutual response SST field, respectively. Therefore, in a climate model, 𝗔 can also be obtained dynamically using ensemble atmospheric experiments as follows. For a prescribed SST forcing field **y** = **m**_{j}, the ensemble mean of the atmospheric experiments [in (2.2)] gives the atmospheric response field of **a**_{j}, that is,

As such, both 𝗕 and 𝗔 provide information on the atmospheric dynamics through the dynamic assessments in (2.10) and (2.16), respectively. The feedback matrix 𝗕 contains the complete information of atmospheric dynamics of *I* × *J* dynamic relations in (2.9). Given 𝗕, we can predict the atmospheric response to any SST field **y**, as in (2.10). In comparison, the total matrix 𝗔 only provides a subset (of *J*) dynamic relations in (2.16), which are associated with the atmospheric model response to *J* SST fields **m**_{j} that are generated in the coupled system.

As discussed in the rest of the paper, most challenging for the nonlocal feedback is the assessment of 𝗕 for a finite sample size. The GEFA in (2.8) is much more sensitive to sampling errors than the scalar EFA used for 𝗔 in (2.13), because the sampling error can increase dramatically in the former by the cross covariance of SST variability. The covarying SST tends to make the SST covariance matrix 𝗖_{YY}(*τ*) singular, which leads to a large 𝗖_{YY}^{−1}(*τ*), and, in turn, a large sampling error in (2.8).^{5}

## 3. Feedback in a simple model

### a. A thermally coupled model

As a pilot study, here we test GEFA by assessing the SST feedback on surface heat flux in a simple thermally coupled ocean–atmosphere model. This model extends the stochastic climate model of Frankignoul and Hasselmann (1977) by including a temperature advection in the atmosphere model in the form of

where *T**_{a} is the surface air temperature, *T** is the SST, and *n** is a stochastic forcing associated with atmospheric internal variability. With a proper nondimensionalization, the coupled system can be rewritten in nondimensional variables as

where

is the downward heat flux, *d* is the oceanic damping, and the model domain is now 0 ≤ *x* ≤ 1. Since the SST is forced locally by the heat flux as shown in the oceanic heat budget equation in (3.1b), the only nonlocal process in the coupled system is the atmospheric advection ∂*T _{a}*/∂

*x*in the atmospheric heat budget equation in (3.1a). The detailed form of this teleconnection is not important for our study here. The important thing is that our model in (3.1) reflects the fact that, at short time scales (monthly to seasonal), remote climate teleconnection is dominated by atmospheric processes (Liu and Alexander 2007). The relative importance of local coupling verse nonlocal advection is measured by

*λ*, with a larger

*λ*representing a stronger local coupling, or a weaker advection. In the limit of strong local coupling (

*λ*→ ∞), the atmospheric heat budget is dominated by the local heat budget 0 ≈

*λ*(

*T*−

*T*) +

_{a}*n*(

*x*,

*t*). For convenience,

*λ*> 0 is used below such that the direction of increasing

*x*is downwind. For clarity of presentation, the oceanic damping is set to

*d*= 0, unless otherwise specified.

At the upstream boundary *x* = 0, the air temperature is set as

which can be thought as a stochastic forcing generated by atmospheric internal variability and land–atmospheric interaction upstream. With this boundary condition, the atmospheric equation in (3.1a) can be integrated as

Dividing the model domain into *I* intervals of the width Δ*x* = 1/*I*, (3.4) can be integrated with the “downwind” scheme as

where, *T _{ai}*(

*t*) =

*T*(

_{a}*x*,

_{i}*t*),

*T*(

_{i}*t*) =

*T*(

*x*,

_{i}*t*), and

*n*(

_{i}*t*) =

*n*(

*x*,

_{i}*t*) are the values at

*x*=

_{i}*i*Δ

*x*. [The result here remains similar with the trapezoid integration scheme (not shown).] The discrete atmospheric equation in (3.5) can be put in the vector form as

where

and a transformed stochastic forcing is introduced as

with

representing the downwind decaying effect. The feedback matrix for *T _{a}* is

with

Physically, the positive diagonal elements represent the positive local response of air temperature to SST. The positive off-diagonal elements represent the nonlocal response, with the *i*th row representing the nonlocal response of air temperature *T _{ai}* to all the upstream SSTs

*T*(

_{j}*j*≤

*i*), whose influence decays with distance as represented by an increasing power of

*q*. Alternatively, the

*j*th column represents the remote impacts of the SST forcing

*T*on all the air temperature downstream

_{j}*T*(

_{ai}*i*≥

*j*), also decaying with the distance.

The SST equation in (3.1b) can also be discretized as

Equations (3.6) and (3.10) form a coupled system of *I* atmospheric variables and *I*(=*J*) SST variables. Substituting (3.6) into (3.10), we have the equation for the coupled system in terms of the SST as

Here, the heat flux has been derived as

and the feedback matrix for heat flux is

Here,

is required for the numerical stability of the coupled system (due to the finite difference in space), because (3.11) has *I* eigenvalues of the same value as −*γ*. Noting 𝗕_{a} in (3.9), the feedback matrix for heat flux 𝗕 can be interpreted as follows. The negative diagonal elements reflect the negative local ocean–atmosphere thermal feedback. [Note that our definition of feedback parameter is of the opposite sign to that of Frankignoul et al. (1998).] A positive SST warms the local air temperature [by a factor of *μ* > 0, as in (3.9)], which reduces the air–sea temperature difference, and, in turn, the downward heat flux (to *γ* ≡ 1 − *μ*; Barsugli and Battisti 1998). The off-diagonal positive elements reflect the nonlocal positive feedback response: along the *j*th column, a positive SST *T _{j}* warms the air temperature downstream, and, in turn, enhances the downward heat flux in the downstream regions (

*i*>

*j*).

### b. The two-point model

It is instructive to start with the two-point model (*I* = *J* = 2), with *i* = 1 and 2 representing the upstream and downstream points, respectively. The feedback matrices for air temperature and heat flux are now, respectively,

The SSTs and, in turn, their cross-covariance matrix 𝗖_{TT}(*τ*) can be calculated analytically (see appendix A). At the limit of short lead (*τ* → 0^{+}), 𝗖_{TT}(*τ*) can be obtained from (A.8) as

where *σ*^{2}_{1}, *σ*^{2}_{2}, *σ*_{12} are proportional to the covariance of the internal variability as defined in (A.5). The mutual response matrix can therefore be obtained from its definition in (2.12) as

As previously discussed, the feedback matrices in (3.14) reflect the local response (diagonal), the downwind decaying of the atmospheric advection (*b*_{a21} = *b*_{21} = *μq*), and the absence of an upwind atmospheric response (*b*_{a12} = *b*_{12} = 0). Note, however, a nonzero upwind SST “mutual response” *m*_{12} ≠ 0, because the long persistence of SST can lead to a positive covariance even for the downwind *T*_{2} leading the upwind *T*_{1}. [As lag increases, however, the upwind influence decays faster (*C*_{12} ∼ *e*^{−2γτ}) than the downwind influence (*C*_{21} ∼ *e*^{−γτ}), as shown in (A.8), because of the preferred downwind atmospheric advection.] Usually, the cross covariance in (3.15), and, in turn, 𝗠 in (3.16), are positive [see the extreme cases discussed following (A.5)]. This reflects the SST tendency of covariability with each other due to the atmospheric advection.

The total response matrices for air temperature and heat flux can be derived with (2.11), (3.14), and (3.16) as

The difference between the total and feedback matrices can be seen directly by comparing (3.17) with (3.14). This difference can also be seen in the examples in Fig. 1, which show all the elements of 𝗕, 𝗕_{a}, 𝗔, 𝗔_{a} as a function of the local coupling parameter *μ* for air temperature (Figs. 1a,c,e,g,i) and heat flux (Figs. 1b,d,f,h,j), in the case of independent stochastic forcings of *n*_{1}, *n*_{2}, and *T*_{a0}. The local response at the upwind point is always identical in the total response and feedback response (*b*_{11} = *a*_{11}, *b*_{a11} = *a*_{a11}; Figs. 1a,b), because of the absence of upwind response for air temperature, and, in turn, heat flux (*b*_{12} = *b*_{a12} = 0). The nonlocal total responses, however, are always more positive than the feedback response for air temperature (*a*_{a12} > *b*_{a12}, *a*_{a21} > *b*_{a21}; Figs. 1c,e). This occurs because the nearby SSTs in the coupled model tend to covary with each other with the same sign (*m _{ij}* > 0; Fig. 1i), which reinforces each other’s warming effect on the air temperature. Interestingly, this positive bias also exists for the local response at the downwind point (

*a*

_{a22}>

*b*

_{a22}; Fig. 1g), reflecting the warm advection effect from the upstream. This implies that the scalar EFA may bias the local feedback in a region that is affected significantly by remote nonlocal feedbacks. For the heat flux, the negative local feedback at the downwind point is reduced (

*a*

_{22}>

*b*

_{22}; Fig. 1h), because the nonlocal warming response to the upstream SST increases the air temperature and in turn reduces the heat loss locally. In the meantime, the positive nonlocal total response are decreased relative to the feedback (

*a*

_{12}<

*b*

_{12},

*a*

_{21}<

*b*

_{21}; Figs. 1d,f), because of the interference of the nonlocal impact from the negative local feedback nearby.

Finally, some insight can be gained on sampling error by examining the magnitude of the SST covariance matrix 𝗖_{TT}. As discussed for (2.8), sampling error tends to increase when 𝗖_{TT} becomes more ill conditioned. The determinant of 𝗖_{TT} is therefore a useful measure of the potential sampling error. For the two-point model, (3.15) gives

This shows that sampling error tends to increase with a stronger advection or a more coherent forcing pattern. A stronger advection (relative to local coupling) increases the sampling error because it reduces det(**C**_{TT}), as seen in the example of Fig. 1j with an increasing *μ*. This occurs because an increased advection enhances the SST correlation downstream. A more spatially coherent stochastic forcing also increases the sampling error. Indeed, in the limit of a perfectly correlated stochastic forcing, one can show [from (A.5)] that *σ*^{2}_{1}*σ*^{2}_{2} − *σ*^{2}_{12} = 0. Therefore, det(𝗖_{TT}) reaches the minimum ∼(*μqσ*^{2}_{1}/2*γ*)^{2}. Furthermore, if the advection intensifies with *μ* → 0, we have det(𝗖_{TT}) → 0. In this limiting case, it is no longer possible to recover the feedback matrix statistically from a coupled simulation. Physically, with all the forcing perfectly correlated, and a strong advection, the SSTs over the entire domain become perfectly correlated. This made it impossible to distinguish the SST impacts from one region to another from a single coupled simulation. In contrast, the covariance of SST does not have obvious impact on the sampling error of 𝗔, because the local estimation in (2.10) does not involve the cross covariance of SST. This suggests that the estimation of 𝗔 is much more stable than that for 𝗕, as will be seen later.

### c. Multipoint models

Now, we study the multipoint models numerically, with the focus on the potential sampling error. With our application to the North Atlantic heat flux feedback in mind, we will use a sample size of *T* = 400, with the data binned in a nondimensional time interval of 0.5, which corresponds roughly to a monthly dataset of 30–40 yr. [The observed SST persistence time is about 2 months in the midlatitude, which corresponds to our model SST persistence time of 1/(1 − *μ*) ∼ *O*(1).]^{6}

We first study a six-point model with *λ* = 2. The stochastic forcings are chosen to be independent of each other, with a standard deviation of *σ*(*n _{i}*) = 10 (

*i*= 1, . . . , 6) in the interior and

*σ*(

*T*

_{a0}) = 1 on the boundary. Figure 2a shows the evolution of three SSTs (at upstream

*i*= 1, midbasin

*i*= 3, and downstream

*i*= 6) in one realization. In spite of the independent forcings, the SSTs show a strong correlation with each other, with the amplitude increasing downstream. The dominant slow variability is well captured by the first principal component (PC), which accounts for 80% of the total variance (Fig. 2b). This strong correlation of SST is caused by the nonlocal atmospheric advection. Otherwise, with independent forcings, the SST variability would have varied independently from each other, with each EOF explaining about 1/6 of the total variance.

The true feedback matrix for heat flux 𝗕 is shown in Fig. 2c (solid line): the elements *b _{ij}* are ordered in a 1D array

*b*with

_{k}*k*=

*i*+ 6(

*j*− 1), increasing first with row

*i*, then with column

*j*. The negative local feedbacks can be seen in the negative spikes, which represent the diagonal elements

*b*=

_{jj}*λ*Δ

*x*− 1 = 2 × (1/6) − 1 = −2/3; the positive nonlocal downwind impact is seen as the positive

*b*that decreases following each negative spike, representing the decay with distance; the absence of nonlocal impact in the upwind direction is represented by the zeros preceding each negative spike. The GEFA estimator 𝗕 with the full data (circles) shows a good agreement with the true 𝗕 in the dominant negative local feedback and the positive downwind impact. The absence of upwind impact is also captured by the small values that are statistically indistinguishable from zero (at the 90% level according to the Monte Carlo test). For comparison, the corresponding total matrix 𝗔 is also plotted (pluses). The strong negative local feedback in 𝗕 becomes weaker in 𝗔, and the nonlocal positive feedback in 𝗕 become significantly negative in 𝗔. These differences are caused by the positive covariability of SST as seen in the time series in Fig. 2a as well as the significant positive off-diagonal elements in

_{k}Given a finite sample, the sampling error of the GEFA estimator tends to increase with the resolution, because the SST field becomes more correlated and therefore the SST covariance matrix 𝗖_{TT} becomes more ill conditioned. One way to reduce the sampling error is to first estimate 𝗕 in the truncated SST EOF space, and then to recover it back into the physical space. This EOF truncation filters out small eigenvalues of 𝗖_{TT} and therefore allows for a better conditioned 𝗖_{TT} in the EOF subspace.

To illustrate the effect of the EOF truncation, we discuss four cases of resolutions *I* = 3, 6, 12, and 24, with the stochastic forcings independent between different points (Fig. 3).^{7} For each resolution, we perform a 20-member ensemble experiments. For each ensemble member experiment, 𝗕 is estimated *I* times using (2.8) as 𝗕_{f} ( *f* = 1, . . . , *I*), with 𝗕_{f} obtained with the leading *f* SST EOFs. The accuracy of each 𝗕_{f} is measured against the true 𝗕 using the ensemble mean of the pattern correlation cor(𝗕_{f} , 𝗕) (Fig. 3a) and the amplitude ratio *σ*(𝗕_{f} )/*σ*(𝗕)(Fig. 3b). For *I* = 3, 6, and 12, the estimator with the full data ( *f* = *I*) is the optimal estimator. The pattern correlation increases with the number of EOFs, peaking with all the EOFs at the value of about 0.8–0.9; the amplitude ratio of 𝗕_{f} also increases toward the true 𝗕 (ratio 1) when almost all the EOFs are retained. In contrast, for the high-resolution case of *I* = 24, the accuracy of 𝗕_{f} decreases after *f* = 15, as seen in both the pattern correlation and amplitude ratio. Therefore, given a sampling size *T*, for sufficiently high resolution, the optimal estimator for 𝗕 is obtained with a truncation of SST EOFs.

It, however, remains unclear to us how to determine the optimal EOF truncation if the true 𝗕 is unknown, as in the case of the observation or a complex coupled model. One possible measure is the successive convergence, based on the successive pattern correlation cor(𝗕_{f−1}, 𝗕_{f} ) (Fig. 3c) and amplitude ratio *σ*(𝗕_{f−1})/*σ*(𝗕_{f} ) (Fig. 3d). Overall, the successive pattern correlation and amplitude ratio appear to increase with the number of EOFs and converges toward 1. Although there is no clear indication of the optimal truncation, the successive pattern correlation and amplitude ratio seem to plateau when the EOFs are increased near the optimal truncation (in the case of *I* = 24).

It is also interesting to examine the corresponding total matrix 𝗔_{f} with *f* SST EOFs retained. As expected, 𝗔_{f} differ significantly from the true 𝗕 at all the resolutions, as seen in pattern correlation (Fig. 3e) and amplitude ratio (Fig. 3f). Nevertheless, 𝗔_{f} converges rapidly with the number of EOF. Indeed, 𝗔_{f} remains virtually the same after the first 3 EOFs, as seen in the successive pattern correlation (Fig. 3g) and amplitude ratio (Fig. 3h). As discussed before, this rapid convergence occurs because the estimation of 𝗔_{f} does not rely on the SST covariance matrix, and therefore, is insensitive to the correlation of SSTs, and, in turn, the resolution.

The accuracy of GEFA also depends on the spatial coherence of the stochastic forcing. This is potentially an important problem for ocean–atmosphere interaction, because of the large spatial scale of intrinsic atmospheric variability. As discussed for the two-point model case, a perfectly correlated stochastic forcing leads to a decreased det(𝗖_{TT}), and in turn a greater sampling error. Figure 2d shows an example of the six-point model similar to that discussed in Fig. 2c. In contrast to the case in Fig. 2c that is forced by independent stochastic forcings, however, this case is forced by a “tripole” stochastic forcing in the interior (with *n*_{1} = *n*_{2}, *n*_{3} = *n*_{4}, and *n*_{5} = *n*_{6}). Now the estimated 𝗕_{6} with the full data (circles in Fig. 2d) becomes much noisier than in the case of independent forcing (circles in Fig. 2c). The estimation, however, is improved significantly with a truncation to three EOFs (asterisk in Fig. 2d), although the EOF truncation seems to smooth the estimation, especially on the negative local feedback spikes.

To illustrate the effect of the pattern of stochastic forcing systematically, we show in Fig. 4 the response matrices forced by three patterns of stochastic forcing: an independent forcing (circle), in which the stochastic forcings are completely independent in the interior and boundary, a tripole forcing (square), (as in Fig. 2d), and a monopole forcing (asterisk), in which *n*_{1} = *n*_{2} = *n*_{3} = *n*_{4} = *n*_{5} = *n*_{6}. Each pattern of forcing is used to generate a 20-member ensemble simulations. For each ensemble member, six 𝗕_{f}s ( *f* = 1, . . . , 6) are estimated with successive truncations of the SST EOF. The accuracy of 𝗕_{f} is measured against the true 𝗕 using the ensemble mean of the pattern correlation (Fig. 4a) and amplitude ratio (Fig. 4b). With the independent forcing, 𝗕_{f} improves monotonically with the addition of EOF, converging toward 1 in both the pattern correlation and amplitude ratio (circle), as in the six-point model case in Fig. 3. In the meantime, 𝗕_{f} becomes more sensitive to the EOF truncation, as indicated by the increased ensemble spread of the pattern correlation. Most dramatically, 𝗕_{f} deteriorates significantly when the forcing becomes a tripole forcing (square): the maximum pattern correlation decreases from over 0.8 to below 0.6, and the optimal truncation for 𝗕_{f} changes from 6 to 3. The ensemble spread seems to also increase with the number of EOFs, significantly when the truncation is beyond the optimal truncation. With a monopole forcing, 𝗕_{f} further deteriorates (asterisk), with the maximum pattern correlation reduced to 0.4 and the optimal truncation is limited to only two EOFs. These examples show that 𝗕_{f} deteriorates when the pattern of the stochastic forcing becomes more spatially coherent. This is consistent with earlier discussions in that a more coherent pattern of forcing generates more coherent SST variability, and, in turn, a more ill-conditioned 𝗖_{TT}, and eventually a greater sampling error. (The increased spatial coherence of SST can be seen from the leading EOF1, which explains about 78%, 84%, and 95% of the total variance for the cases of independent forcing, tripole forcing, and monopole forcing, respectively.)

The convergence of 𝗕_{f} also deteriorates when the forcing becomes more spatially coherent as seen in the successive pattern correlation and amplitude ratio (Figs. 4c,d). The ensemble spreads increase significantly for both the correlation and amplitude ratio, from the case of independent forcing toward the case of tripole forcing and then the monopole forcing. In the latter two cases, the successive pattern correlation and amplitude ratio exhibit a maximum near the optimal truncation, indicating a minimum sensitivity of the GEFA estimator 𝗕 with respect to the EOF truncation. Therefore, the optimal truncation appears to be the case near the maximum successive pattern correlation. As discussed in Fig. 3, the total matrix 𝗔 differs substantially from the true 𝗕, but converges rapidly as EOFs are added (not shown).

In short, the simple model study shows that, for a finite sample size, GEFA provides a reasonable estimation of 𝗕 at low resolutions. With a sufficiently high resolution, or large-scale stochastic forcing, however, the accuracy of GEFA decreases significantly due to the nature of nonlocal estimation. Nevertheless, an optimal estimator seems to be available, in principle, with certain truncation of SST EOFs. In comparison, the total matrix 𝗔 is stable due to the nature of the local estimation, but it could differ from the feedback matrix 𝗕 substantially.

## 4. Summary and discussion

Due to nonlocal climate dynamics, a general atmospheric feedback response to SST (or other lower boundary) forcing consists of nonlocal as well as local responses. In the context of linear dynamics, this general feedback, in principle, can be represented by a feedback matrix 𝗕—effectively the response of the Green function. This matrix provides important information on the dynamics of the atmospheric response, and, in turn, the feedbacks in the coupled system. The major objective of this paper is to develop a statistical method to assess both nonlocal and local climate feedbacks, with the focus on the estimation of 𝗕. Here, the univariate EFA of Frankignoul et al. (1998) is generalized to a multivariate EFA as GEFA. GEFA is then used to assess ocean–atmosphere thermal feedback in a simple ocean–atmosphere model. It is shown that GEFA is able to extract the nonlocal feedback associated with the downstream atmospheric teleconnection. Furthermore, physical insight can also be gained by comparing the feedback matrix 𝗕 with the total matrix 𝗔. Unlike the feedback matrix that identifies the nonlocal feedback response to independent SST variability at different regions, 𝗔 represents the total feedback impact from all the covarying SSTs, and therefore its implication on nonlocal dynamics needs to be treated with caution.

Given a finite sample size, it remains challenging to obtain the optimal estimator for the feedback matrix 𝗕 in a coupled climate system and the observation. In general, sampling error increases significantly with the spatial resolution of the data. Therefore, in a coupled system with highly correlated SST variability, the accuracy of the GEFA estimator can be limited intrinsically. Of course, for the application to a complex system such as the real world or a general circulation climate model, the method itself is also limited by its assumptions, notably linearity.

The GEFA has also been applied to the assessment of nonlocal feedback between turbulent heat flux and SST for the observed North Atlantic. At low resolutions (three–six regions over the North Atlantic), GEFA confirms the dominant local negative SST feedback on heat flux, and, furthermore, identifies a nonlocal feedback, with a warm Gulf Stream SST enhancing the downward heat flux downstream in the subpolar region. The detailed results will be presented elsewhere.

It is important to point out that the nonlocal assessment can be applied to much more general issues to understand the climate feedback to multiple factors, if the boundary variable is replaced by a mixture of variable factors. In particular, when these factors are correlated with each other, the scalar EFA is no longer valid. For example, we can study the atmospheric response to Niño-3 SST, North Pacific SST, Indian Ocean SST, and Eurasian land vegetation cover. Since the North Pacific SST and tropical Indian Ocean SST include some responses to, and therefore are correlated with, Niño-3 SST, a scalar EFA may not produce the correct estimate of nonlocal feedbacks. Instead, GEFA provides a promising method. It is also sometimes more relevant to perform the nonlocal feedback analysis in the EOF space itself, instead of the physical space, because some major climate modes are better defined with EOFs. In these cases, even with a small number of degrees of freedom, the feedback matrix will provide great insight into the roles of forcing of different climate factors. All these suggest the need for further improvement on the assessment methods for nonlocal climate feedbacks in the future.

Finally, it is important to compare and evaluate GEFA with other relevant statistical methods for nonlocal feedbacks, such as the maximum covariance analysis (MCA; Czaja and Frankignoul 2002), and especially the linear inverse method (LIN; Penland and Sardeshmukh 1995), or more generally, the fluctuation–dissipation theorem (FDT; Leith 1975; von Storch 2004; Kirk-Davidoff 2005, manuscript submitted to *J. Climate*; Gritsun and Branstator 2007). These works are in progress. Here, a brief discussion based on our preliminary results is given in appendix B.

## Acknowledgments

We are indebted to Dr. C. Frankignoul for his thoughtful comments on an earlier version of the manuscript. We thank Dr. Zhiguang Qian for helpful comments on the paper. We also thank Drs. B. Kirk-Davidoff, A. Mjada, G. Branstator, and A. Gritsun for informing us (including their unpublished manuscripts) of the fluctuation–dissipation theorem at the end stage of the work. Comments from three anonymous reviewers have helped the presentation of the paper. This work is supported by DOE, NSF, and NOAA. This work is also partially supported by Chinese NSF and the Ocean University of China.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### The Two-Point Coupled Model Solution

With two points, the coupled system in (3.11) is reduced to

where the transformed stochastic forcing is related to the stochastic forcing in the interior and upstream boundary as

Denote the cross covariance of the stochastic forcings as

The cross covariance of the transformed stochastic forcing can be derived as

with

In the extreme case of independent stochastic forcings, *ω*_{12} = *ω*_{1a} = *ω*_{2a} = 0, (A.5) is simplified to *σ*^{2}_{1} = *ω*^{2} + *q*^{2}*ω*^{2}_{a}, *σ*^{2}_{2} = (1 + *q*^{2})*ω*^{2} + *q*^{4}*ω*^{2}_{a}, and *σ*_{12} = *qω*^{2} + *q*^{3}*ω*^{2}_{a}. In another extreme case with all the stochastic forcings perfectly correlated, *ω*_{12} = *ω*^{2} and *ω*_{1a} = *ω*_{2a} = *ωω _{a}*, (A.5) degenerates to

*σ*

^{2}

_{1}= (

*ω*+

*qω*)

_{a}^{2},

*σ*

_{12}= (

*ω*+

*qω*)[(1 +

_{a}*q*)

*ω*+

*q*

^{2}

*ω*], and

_{a}*σ*

^{2}

_{2}= [(1 +

*q*)

*ω*+

*qω*]

_{a}^{2}.

### APPENDIX B

#### Relation of GEFA with Other Methods

Although a detailed comparison is beyond the scope of this paper, here, we briefly discuss the relation between the GEFA and other relevant statistical methods.

##### Relation with multivariate regression

Following Z. Liu et al. (2006), the feedback matrix estimated using GEFA in (2.8) can be shown to be related to the regression as

where 𝗥(*τ*) = 𝗖_{XY}(*τ*)𝗖^{−1}_{YY}(0) is the multivariate regression matrix of the predictive equation 𝗫_{t} = 𝗥(*τ*)𝗬_{t−τ} + 𝗘 and Γ(*τ*) = 𝗖_{YY}(*τ*)𝗖^{−1}_{YY}(0) is the autocorrelation matrix of the SST, all at the same lag *τ*. Physically, (B.1) states that the predictability of the atmospheric response to SST with a lead *τ*equals the instantaneous atmospheric response multiplying its persistence of lag *τ* following the SST. The relation in (B.1) can also be obtained by assuming the SST as a first-order autoregressive process (AR-1). This relation may be useful for studying the statistical properties of the feedback matrix.

##### Relation with MCA

Czaja and Frankignoul (2002) studied the feedback of SST on the atmospheric geopotential height field using an MCA analysis, in which the singular value decomposition (SVD; Strang 1976; Bretherton et al. 1992) is applied to the lagged cross-covariance matrix 𝗖_{XY}(*τ*) between the atmosphere and SST. The MCA analysis provides feedback information very different from the feedback matrix 𝗕(*τ*) = 𝗖_{XY}(*τ*)𝗖^{−1}_{YY}(*τ*). The former is specific to a coupled system because the SST variability is determined by the coupled system, somewhat similar to the total matrix 𝗔, while the latter only concerns the atmospheric dynamics, and is discussed below.

##### Optimal feedback response

One can show that the SVD decomposition of the feedback matrix 𝗕 gives the optimal SST forcing fields (right vectors) and the corresponding optimal atmospheric response (left vectors). Therefore, with 𝗕, one can identify the optimal SST pattern that generates the maximum atmospheric response. Furthermore, this optimal feedback decomposition provides another possible filter of the noise of the feedback matrix 𝗕. Our preliminary analysis suggests that, at high resolution, while the 𝗕 itself becomes noisy, its leading SVD modes remain relatively stable. Therefore, in spite of the noisy 𝗕, large-scale feedback responses may still be extracted from it with reasonable accuracy.

##### Relation with LIN/FDT

In principle, the feedback matrix can be derived directly from the coupled system

At climate time scales longer than the rapid atmospheric adjustment, the atmosphere equation in (B.2a) can be approximated as a quasi-equilibrium state

which then leads to the feedback matrix as 𝗕 = −𝗚^{−1}_{XX}𝗚_{XY}. The matrix for the coupled system 𝗚, in principle, can be obtained using the LIN (Penland and Sardeshmukh 1995), or more generally, the FDT (Leith 1975; von Storch 2004; Kirk-Davidoff 2005, manuscript submitted to *J. Climate*; Gritsun and Branstator 2007). It is therefore interesting to compare GEFA with LIN/FDT. Our preliminary study suggests that, for the same accuracy, the LIN/FDT requires data of higher temporal resolution (e.g., daily) that resolves the fast temporal dynamics of the atmosphere [such as (B.2a)]. GEFA, however, only makes use of the low-frequency data (e.g., monthly) as it only deals with the quasi-equilibrium response [such as (B.3)]. In the meantime, GEFA cannot assess the full temporal dynamics of the atmosphere. Rewriting the full atmospheric equation in (B.2a) as

the LIN/FDT estimates both 𝗚_{XX} and 𝗕, while GEFA is only able to assess 𝗕.

Finally, we point out an interesting difference between the feedback matrix 𝗕 and the total matrix 𝗔 from the perspective of the full coupled equations in (B.2a, b): 𝗕 is determined by the atmospheric dynamics, 𝗚_{XX} and 𝗚_{XY}, while 𝗔 usually depends on the entire coupled system (i.e., the entire 𝗚). Given an atmosphere (𝗚_{XX} and 𝗚_{XY}), **B** is determined uniquely, but 𝗔 may be different for different ocean dynamics (𝗚_{YX} and 𝗚_{YY}).

## Footnotes

*Corresponding author address:* Z. Liu, CCR, University of Wisconsin—Madison, 1225 W. Dayton St., Madison, WI 53706. Email: zliu3@wisc.edu

* Center for Climatic Research Contribution Number 919.

^{1}

For application to the observed climate variability such as the North Atlantic Oscillation (NAO) and Pacific–North America (PNA), the long persistence in the final total variability *x* is caused dominantly by the response to the slow SST [by (*t*) in (1.1)] instead of the internal variability [*n*(*t*) in (1.1)], which is unknown. Therefore, a lag-1 month is usually sufficient to allow (1.2) to be valid.

^{2}

In practice, in atmospheric general circulation models, one chooses a small region near point *j*, instead of the single point *j* (e.g., Gritsun and Branstator 2007).

^{3}

Since 〈*n _{i}*(

*t*),

*y*(

_{j}*t*−

*τ*)〉 = 0, multiplying

*y*(

_{j}*t*−

*τ*) in (2.1) and ensemble averaging yield 〈

*x*(

_{i}*t*),

*y*(

_{j}*t*−

*τ*)〉 = Σ

^{J}

_{k=1}

*b*〈

_{ik}*y*(

_{k}*t*),

*y*(

_{j}*t*−

*τ*)〉. Dividing both sides by 〈

*y*(

_{j}*t*),

*y*(

_{j}*t*−

*τ*)〉, and noticing (2.10) and (2.12), we have

*a*= Σ

_{ij}^{J}

_{k=1}

*b*, or (2.11).

_{ik}m_{kj}^{4}

The *m _{kj}* may appear to be derived using the scalar EFA as a mutual response between the SST pair

*k*and

*j*, with

*y*responding to

_{k}*y*, from an equilibrium response

_{j}*y*(

_{k}*t*) =

*m*(

_{kj}y_{j}*t*) +

*n̂*(

_{kj}*t*). Strictly speaking, however, it is not an SST response, because the SSTs all vary at the slow time scale and therefore the equilibrium response is not valid.

^{5}

Similar difficulties arise if one uses the equivalence relation in (2.11) to estimate 𝗕 as 𝗕 = 𝗔𝗠^{−1}. Both 𝗔 and 𝗠 are obtained with the local method. However, since the determinants of 𝗖_{YY}(*τ*) and 𝗠 are related to each other as det(𝗖_{YY}) × [Π^{J}_{j=1}〈*y _{j}*(

*t*),

*y*(

_{j}*t*−

*τ*)〉] = det(𝗠), an ill-conditioned 𝗖

_{YY}(

*τ*) also corresponds to an ill-conditioned 𝗠, and in turn a greater sampling error in 𝗕.

^{6}

Nevertheless, one should be cautious about the direct implication of the simple model results here to the observed North Atlantic in the next section, because the parameter regimes and the physical processes may differ significantly.

^{7}

An increase in resolution also reduces the accuracy of the estimation of 𝗕 by reducing the persistence time of the SST (Liu et al. 2006), which is proportional to the reciprocal of *γ* = 1 − *λ*Δ*x*. To separate this effect from the truncation of SST EOFs, we set the effective persistence time the same as in the case of a three-point model by introducing a negative ocean damping in (3.1b) as *d* = *λ*(Δ*x* − 1/3).