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.
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 ni(t) as the local atmospheric internal variability. The bij represents the impact of the jth SST component on the ith 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 𝗬Tt−τ 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(τ)𝗖−1YY(τ) ≠ 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 bij, 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 yj(t).2 With 〈ni(t)〉 = 0, the ensemble mean of the atmospheric model response in (2.1) gives 〈xi(t)〉 = bijyj(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 (xi) and forcing (yj) as
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 shown3 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 mjj = 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, aij can be thought as the ith atmospheric response to the jth SST, only if yj is the sole forcing on xi in the coupled climate system such that
In a general coupled system, however, xi is also affected by SSTs of other regions and therefore aij does not represent the true response of xi to yj alone. Rather, since all the SSTs vary slowly, at short lags, mkj can be approximated as the regression coefficient rkj as
Therefore, (2.12) suggests that aij could be thought as the total feedback response on xi from all the SSTs that covary with the SST yj in the coupled system. This is in contrast to the feedback bij that identifies the SST forcing impact on xi from the SST yj 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.
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:
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 = mj, the ensemble mean of the atmospheric experiments [in (2.2)] gives the atmospheric response field of aj, 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 mj 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
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 ∂Ta/∂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 − Ta) + 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, Tai(t) = Ta(xi, t), Ti(t) = T(xi, t), and ni(t) = n(xi, t) are the values at xi = 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
and a transformed stochastic forcing is introduced as
representing the downwind decaying effect. The feedback matrix for Ta is
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 ith row representing the nonlocal response of air temperature Tai to all the upstream SSTs Tj (j ≤ i), whose influence decays with distance as represented by an increasing power of q. Alternatively, the jth column represents the remote impacts of the SST forcing Tj on all the air temperature downstream Tai (i ≥ j), also decaying with the distance.
The SST equation in (3.1b) can also be discretized as
Here, the heat flux has been derived as
and the feedback matrix for heat flux is
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 jth column, a positive SST Tj 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,
As previously discussed, the feedback matrices in (3.14) reflect the local response (diagonal), the downwind decaying of the atmospheric advection (ba21 = b21 = μq), and the absence of an upwind atmospheric response (ba12 = b12 = 0). Note, however, a nonzero upwind SST “mutual response” m12 ≠ 0, because the long persistence of SST can lead to a positive covariance even for the downwind T2 leading the upwind T1. [As lag increases, however, the upwind influence decays faster (C12 ∼ e−2γτ) than the downwind influence (C21 ∼ 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 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 n1, n2, and Ta0. The local response at the upwind point is always identical in the total response and feedback response (b11 = a11, ba11 = aa11; Figs. 1a,b), because of the absence of upwind response for air temperature, and, in turn, heat flux (b12 = ba12 = 0). The nonlocal total responses, however, are always more positive than the feedback response for air temperature (aa12 > ba12, aa21 > ba21; Figs. 1c,e). This occurs because the nearby SSTs in the coupled model tend to covary with each other with the same sign (mij > 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 (aa22 > ba22; 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 (a22 > b22; 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 (a12 < b12, a21 < b21; 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(CTT), 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 σ21σ22 − σ212 = 0. Therefore, det(𝗖TT) reaches the minimum ∼(μqσ21/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 σ(ni) = 10 (i = 1, . . . , 6) in the interior and σ(Ta0) = 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 bij are ordered in a 1D array bk with 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 bjj = λΔx − 1 = 2 × (1/6) − 1 = −2/3; the positive nonlocal downwind impact is seen as the positive bk 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
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 n1 = n2, n3 = n4, and n5 = n6). 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 n1 = n2 = n3 = n4 = n5 = n6. Each pattern of forcing is used to generate a 20-member ensemble simulations. For each ensemble member, six 𝗕fs ( 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.
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.
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
In the extreme case of independent stochastic forcings, ω12 = ω1a = ω2a = 0, (A.5) is simplified to σ21 = ω2 + q2ω2a, σ22 = (1 + q2)ω2 + q4ω2a, and σ12 = qω2 + q3ω2a. In another extreme case with all the stochastic forcings perfectly correlated, ω12 = ω2 and ω1a = ω2a = ωωa, (A.5) degenerates to σ21 = (ω + qωa)2, σ12 = (ω + qωa)[(1 + q)ω + q2ωa], and σ22 = [(1 + q)ω + qωa]2.
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
where 𝗥(τ) = 𝗖XY(τ)𝗖−1YY(0) is the multivariate regression matrix of the predictive equation 𝗫t = 𝗥(τ)𝗬t−τ + 𝗘 and Γ(τ) = 𝗖YY(τ)𝗖−1YY(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(τ)𝗖−1YY(τ). 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 𝗕 = −𝗚−1XX𝗚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).
Corresponding author address: Z. Liu, CCR, University of Wisconsin—Madison, 1225 W. Dayton St., Madison, WI 53706. Email: firstname.lastname@example.org
* Center for Climatic Research Contribution Number 919.
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.
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).
The mkj may appear to be derived using the scalar EFA as a mutual response between the SST pair k and j, with yk responding to yj, from an equilibrium response yk(t) = mkjyj(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.
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) × [ΠJj=1〈yj(t), yj(t − τ)〉] = det(𝗠), an ill-conditioned 𝗖YY(τ) also corresponds to an ill-conditioned 𝗠, and in turn a greater sampling error in 𝗕.
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.
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).