Sensitivity of Surface Temperature to Oceanic Forcing via q-Flux Green’s Function Experiments. Part I: Linear Response Function (original) (raw)

1. Introduction

Despite the spatial homogeneity of the increase of greenhouse gas (GHG) concentration, the surface temperature response to GHG forcing shows rich structures in space. A number of studies have been dedicated to understanding the spatial pattern of the sea surface temperature (SST) response to the GHG forcing through examining the surface energy budget (Xie et al. 2010; DiNezio et al. 2009; Lu and Zhao 2012; Luo et al. 2015; Liu et al. 2017), a perspective that might fall into the so-called surface energy budget fallacy. As noted by Pierrehumbert (2010), “if one is in a regime where the surface fluxes tightly couple the surface temperature to the overlying air temperature as is in the low-frequency SST evolution, there is no need to explicitly consider the surface balance in determining how much the surface warms or cools.” For slowly evolving SST, the surface energy budget terms and the oceanic convergence or divergence always conspire to ensure that the SST is approximately in equilibrium, making it difficult to assign causality based solely on which energy budget term contributes the most to the tendency. A more constructive approach to understanding the SST response to external forcing may be through the linear response function (LRF) of the SST to the surface energy flux. Considering the key feature of the coupled ocean–atmosphere system that the low-frequency (decadal and longer) variations of SST and the intrinsic atmospheric time scales from days to weeks are well separated, one can have a significant simplification in the treatment of the slow SST response or variations (see Palmer 1999; Alexeev 2003; Winton et al. 2010; Liu et al. 2012). The SST (δT) can be thought of as in equilibrium, resulting from the constant balance between the oceanic convergence or divergence of heat (δf; hereinafter referred to as q flux) and the surface feedbacks:

e1

e1

where the surface energy flux feedbacks are expressed to be linearly proportional to the SST perturbation. Note that matrix

is the LRF that represents all the relevant physical processes such as evaporative damping effect (Knutson and Manabe 1995), wind–evaporation–SST feedback mechanism (e.g., Xie and Philander 1994), Planck feedback, cloud–radiation feedbacks (e.g., Ramanathan and Collins 1991), and so on. Since

has taken into account all the (linear) relationships between the SST anomalies and the feedback processes and is a function of the mean climate, one can obtain the linear SST response to a given q flux without resorting to energy budget analysis. Conceptually,

can be understood as the linear Jacobian of the full operator evaluated at a reference state, and the nonlinear component in the operator is either ignored or treated as a nonlinear noise forcing. Note also that the atmospheric dynamics can spread the heat from a local q flux via teleconnection mechanisms, so a local SST response is the product of the _q_-flux forcing everywhere, and the mutual relationship between SST and the q flux at low enough frequency can be all encapsulated in the matrix

.

Constructing the LRF of the SST to q flux is the first goal of Part I of this study. Through this LRF, we will be able to predict the global SST response to a given arbitrary q flux without running expensive model experiments. Moreover, the LRF of the SST can provide valuable information through its neutral vector (Marshall and Molteni 1993; Goodman and Marshall 2002, 2003; Watanabe and Jin 2004; Hassanzadeh and Kuang 2016a). The neutral vector is the right singular vector of the LRF with the smallest singular value, and hence it is most likely to manifest in response to a random _q_-flux forcing. Therefore, the patterns of neutral vectors or the combination thereof are expected to be the most excitable under climate change forcings; under stochastic atmospheric forcing, they may also emerge as the most prevailing patterns of the low-frequency internal variability. Several studies (e.g., Watanabe and Jin 2004; Ring and Plumb 2008; Hassanzadeh and Kuang 2016a) have used the neutral vector analysis and identified that the annular mode–like patterns are the most excitable dynamical mode in the atmosphere. This prompted us to examine the neutral vector of the SST in both an AGCM with a slab ocean and a coupled climate model with ocean dynamics to see if it can be identified with the established modes of the SST variability. In Part II (Liu et al. 2018, manuscript submitted to J. Climate), the global surface temperature response to a localized oceanic forcing can be interpreted as the sum of the excited leading neutral modes, which, in turn, organizes the pattern of the radiative feedbacks in such a way to account for the distinct global warming sensitivity for different forcing locations.

What equally motivates this study is the urge to construct the sensitivity map of any given climate index and/or phenomenon of interest to the q flux in the same spirit as Barsugli et al. (2006, hereinafter B2006). Through the Green’s function (GRF) approach, B2006 identified the areas in the tropical Indo-Pacific region where the SST anomalies can most effectively drive the global mean surface warming and global increase of precipitation. Here we follow suit from B2006, except that now the inquiry is about the sensitivity to q flux instead of SST. This is also motivated by the latest proposition that the ocean heat uptake (represented by negative q flux) in the Southern Ocean and North Atlantic is key to delaying global surface temperature warming under a sudden increase of CO2, owing to the greater effectiveness of the Southern Ocean heat uptake in cooling the global climate than the warming effect of CO2 per unit of atmospheric energy perturbation (e.g., Winton et al. 2010; Armour et al. 2013; Rose and Rayborn 2016). It is the second goal of the current study to investigate the global warming or cooling effectiveness of the ocean heat uptake everywhere in the global ocean by constructing a global surface warming sensitivity map; this will be achieved via the LRF of surface temperature to q flux. In Part II (Liu et al. 2018, manuscript submitted to J. Climate), we will delve into the specific feedback processes responsible for the spatial dependence of the global warming sensitivity using a radiative kernel.

There are at least two ways to construct LRF: fluctuation–dissipation theorem (FDT) approach and GRF approach, with their corresponding LRFs denoted as and , respectively. FDT states that by observing the statistical characteristics of the unperturbed system, an LRF can be estimated that gives the response of the quasi-Gaussian system to weak external forcing. It has a relatively long history of use as a powerful tool in the climate science research, beginning with the seminal work of Leith (1975), followed, among others, by Bell (1980), North et al. (1993), Gritsun and Branstator (2007), Liu et al. (2012), Lutsko et al. (2015), Fuchs et al. (2015), and Hassanzadeh and Kuang (2016b). Some recent studies have further proposed a time-periodic (seasonal) operator with more accuracy based on larger data samples (e.g., Majda and Wang 2010; Gritsun 2010; Fuchs and Sherwood 2016). In the same spirit as Langen and Alexeev (2005), who found some success in constructing an LRF of SST using the FDT approach for a zonally symmetric aquaplanet AGCM, we here apply the similar method to a realistic AGCM.

A more straightforward approach to construct the LRF is the GRF approach (e.g., Branstator 1985; Kuang 2010; Hassanzadeh and Kuang 2016a). Here, we perform a large set of simulations with a state-of-the-art AGCM coupled to a slab ocean model using an array of localized _q_-flux anomaly patches that cover most of the global ocean surface (Fig. 1). Thus, the SST response to an individual _q_-flux patch can be considered akin to a Green’s function for the _q_-flux forcing in that location, and the SST responses to all imposed _q_-flux patches can then be used to construct .

Fig. 1.

Fig. 1.

Fig. 1.

Configuration of _q_-flux perturbation patches, each being illustrated by the 6 W m−2 contours. Note that the size of the patch is actually larger than the contoured area.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

Hypothetically, if data samples are long enough for the FDT approach, and the forcing for the GRF approach is small enough, the two methods should converge to the same result for a Gaussian system. However, they could depart from each other because of reasons like (i) departure from Gaussianity of the SST statistics (Loikith and Neelin 2015; Berner and Branstator 2007), or (ii) sampling errors resulting from limited data length. Thus, as a part of the tasks of this study, we also evaluate the LRFs calculated from FDT and GRF, and explore the possible sources of errors and the repercussions in their applications. Finally, we restrict our focus of investigation to the SST–q flux relationship in this study, leaving the problem of air–sea interaction, which is often framed as generalized equilibrium feedback system (e.g., Frankignoul et al. 1998; Liu et al. 2008, 2012) or linear inverse model (e.g., Penland 1989; Penland and Sardeshmukh 1995; Winkler et al. 2001), to future investigation.

The rest of the paper is structured as follows. Section 2 describes a long control slab model simulation for constructing and the design of the _q_-flux patches for the estimate of , as well as a test case to assess the accuracy of the LRFs. Sections 3 and 4 present the detailed procedure to construct and assess the two LRFs. Section 5 discusses the sources of the inaccuracy in the LRF built from the FDT approach (). In section 6, we demonstrate two applications of the estimated LRF: (i) creating a sensitivity map of global SST response and the related feedbacks and (ii) calculating the neutral vector of the SST variations with possible implications on the leading multidecadal modes of the SST variability. Finally, the paper will conclude with a summary and discussion.

2. Model and experiment

The model we use is Community Earth System Model, version 1.1 (CESM1.1), coupled with a slab ocean model (CESM–SOM), in which the active Community Atmosphere Model, version 5 (CAM5), the Community Land Model, version 4 (CLM4), and the Community Ice Code (CICE) are included with the CESM–SOM. The horizontal resolution of CAM5 and CLM4 is 2.5° longitude × 1.9° latitude, with the atmospheric component having 30 vertical levels. The horizontal resolution of the CICE and SOM is at a nominal 1°, telescoped meridionally to approximately 0.3° at the equator. In this model, the ocean and atmosphere are only thermodynamically coupled and SST is computed from surface heat flux and q flux that accounts for the missing ocean dynamics.

The construction of requires a long unperturbed simulation, so we make use of a long, 900-yr CESM–SOM control simulation (CTRL) from NCAR, which is forced with preindustrial carbon dioxide levels and solar insolation. Both the mixed layer depth and the q flux used in CTRL are derived from the climatology of a fully coupled CESM1.1 control simulation. We also make use of a 1000-yr coupled CESM1.1 control simulation to validate the significance of the neutral vector obtained with CESM–SOM in section 6b.

The creation of

needs a large number of simulations. Therefore, starting from January 1st of an arbitrary year from the equilibrated CTRL run, we perform 106 pairs of “warm patch” and “cold patch” simulations, in which the _q-_flux patch is added to or subtracted from the slab ocean. The locations of the 106 _q_-flux patches are illustrated in Fig. 1. Following Barsugli and Sardeshmukh (2002), the _q_-flux patch is designed as a localized cosine hump:

e2

e2

where Q is the maximum _q_-flux anomaly of ±12 W m−2 at the center of the rectangular patch (ϕ k ± ϕ w, λ k ± λ w), and ϕ w = 30° and λ w = 12° are half-widths of the rectangle in zonal and meridional directions, respectively. Note that what is contoured in Fig. 1 is half of the _q_-flux anomaly maximum (Q/2), and the adjacent patches actually overlap. The sum of the q flux from all the patches amounts to a uniform heating (cooling) of 12 W m−2 over the global ocean except near the edges of the oceans. Given that approximately 20 years are needed for the CESM–SOM to approach equilibrium, each of the patch experiments is integrated for 40 years, and only the data of the last 20 years are averaged and compared to that from the CTRL to give the anomalous SST response. The linear and nonlinear parts of the SST anomalies can be estimated as

and

, where the subscript + (−) denotes the response in +12 (−12) W m−2 simulation.

It should be noted that the choice of the forcing amplitude of 12 W m−2 and the simulation length of 40 years is a trade-off between the requirement for generating a robust response signal (which entails long simulations and large forcing perturbation) and the need to stay within the linear regime (which entails long simulations and small forcing perturbation) within the computational affordability. As it turns out, the averaged magnitude of local SST response at each patch is in the order of 0.5 K (see Fig. 4c), which is comparable to B2006, who used SST patches with a mean of 0.667 K. However, we also notice that the 12 W m−2 amplitude of the _q_-flux perturbation may not be small enough to ensure the linear regime.

To test the accuracy of the constructed and , we further perform a pair of AGCM–SOM experiments (labeled CPX) forced by a _q_-flux perturbation with more complexity (see Fig. 6a), which was originally designed to add heat to the equatorial Pacific to correct the cold tongue bias there (e.g., Mechoso et al. 1995). (The in CPX averaged over years 21–40 is shown in Fig. 6b.) In addition, the patch experiments themselves can also serve as test cases for the LRFs.

3. Construction of the LRFs

a. FDT approach

Following previous studies (Gritsun and Branstator 2007; Lutsko et al. 2015; Fuchs et al. 2015; Hassanzadeh and Kuang 2016b) we use the most common Gaussian formulation of FDT to estimate

as follows:

e3

e3

where

represents covariance matrix, τ is the time lag, and

(0) is the autocovariance matrix. Equation (3) could also be applied to a quasi-Gaussian system (Majda and Wang 2010; Gritsun and Branstator 2007), which is more closely representative of the climate system. Still, the significant departure of the SST distribution from Gaussianity, as found in both observations (e.g., Loikith and Neelin 2015) and model simulations (e.g., Berner and Branstator 2007), may impair the accuracy of the FDT approach. Moreover, although we aim to estimate the equilibrium SST response, it is impractical to set the upper bound of the integral infinite, so we use a finite but sufficiently large upper boundary _N_Inf = 120 months to capture the time scale of the response. Further tests show that increasing _N_Inf does not appreciably improve the performance of the FDT operator (not shown).

To construct for atmospheric variability, daily data is needed, as the dominant atmosphere modes manifest themselves mostly at intraseasonal time scales (e.g., Feldstein, 2000; Hoerling and Kumar 2003). On the other hand, the SST modes such as El Niño–Southern Oscillation (ENSO; Philander 1990), interdecadal Pacific oscillation (IPO; Zhang et al. 1997), and Atlantic multidecadal oscillation (AMO; Kerr 2000), are all characteristic of interannual to multidecadal time scales, so monthly data are sufficient for constructing the SST . To test the sensitivity to sampling frequency, we repeat the construction using daily SST, and the result shows little difference from that using monthly data, whereas the computing cost of the former approach is an order of magnitude greater. Thus, all the results reported here are based on monthly data.

b. Dimension reduction

Although our model resolution is already relatively coarse, it still has over 8000 ocean grid points globally, so it is impractical to calculate the covariance matrix of that size, and dimension reduction is necessary. As conventionally practiced (Gritsun and Branstator 2007; Lutsko et al. 2015; Fuchs et al. 2015; Hassanzadeh and Kuang 2016b), we project the raw SST data onto the empirical orthogonal functions (EOFs) of the SST in the long CNTRL simulation, truncated based on North’s criterion (North et al. 1982) for well-resolved EOFs. The resultant operator and LRF are all expressed based on the EOFs with reduced dimensions. In addition, dimension reduction can largely reduce the effective number of degrees of freedom _N_eff and reduce the likelihood that (0) becomes ill-conditioned (Martynov and Nechepurenko 2006).

Some studies (e.g., Gritsun and Branstator 2007; Fuchs et al. 2015; Fuchs and Sherwood 2016) simply used the leading EOFs that explain a certain proportion (varies between 90% and 97.5%) of the total variance without a rigorous assessment of the number of well-resolved EOFs. Here, we adopt North’s criterion that demands the ratio of sampling error of a particular eigenvalue λ i to the spacing between neighboring eigenvalues to be less than 1:

e4

e4

where

is the sampling error ratio. The effective degree of freedom is estimated following Trenberth (1984):

e5

e5

where ρ(τ) is the autocorrelation of each principal component (PC). According to Fig. 2, the first few EOFs easily satisfy the C r < 1 condition with well-separated eigenvalues despite relatively small _N_eff. However, C r increases gradually with EOF modes as the separation between the eigenvalues decreases much faster than the increase of _N_eff. A power law function fit (orange line in Fig. 2) to the raw C r curve (black line in Fig. 2) and its intercept with line C r = 1 are used to identify the reasonably resolved EOFs. As a result, the first 62 EOFs, which explain approximately 80% of the total variance, are retained to construct the

. Sensitivity tests using either more or fewer than 62 EOFs for the base of

(not shown) all result in deteriorate results (measured against the actual response to the _q_-flux patch experiments). All in all, both North’s criterion and empirical tests point to 62 EOFs as the optimal number for dimension reduction.

Fig. 2.

Fig. 2.

Fig. 2.

North’s criterion C r (black line) as a function of the number of EOFs, overlaid with its power law function fit (orange line).

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

c. Green’s function approach

The mean SST response to an individual patch shown in Fig. 1 can be considered to be the GRF for the forcing in that location. Therefore, the equilibrium SST response to an arbitrary _q_-flux forcing pattern δf could be approximated as a weighted combination of the GRF responses to the local _q_-flux patches:

e6

e6

where

is a rectangular matrix comprising all the Green’s function. Apparently, −

(hereinafter referred to as

) is the LRF derived from GRF approach. Following Barsugli and Sardeshmukh (2002), each column vector Gk can be calculated as

e7

e7

where k, ranging from 1 to 106, is the index of the _q_-flux patches,

is the linear parts of global SST response to the _k_th forcing, and x j represents the grids within the _k_th patch (see Barsugli and Sardeshmukh 2002; Barsugli et al. 2006 for details). Since

is not a square matrix and can only be calculated from pseudoinversion, we project the SST response and _q_-flux forcing onto the selected EOFs to obtain a square matrix for

. Employing this dimension reduction will result in two major benefits. First, the comparison between the performances of

and

will show if dimension reduction is a source of error in constructing LRF. Second, dimension reduction enables

to have the same matrix dimension (62 × 62) as

, thus facilitating the comparison of their different eigen properties.

d. Comparison of eigenvalues

Figure 3 shows the eigenvalues λ of and , with the decay rate (real part of λ) of each eigenmode on the ordinate, and the frequency (imaginary part of λ) on the abscissa. Note that the growing eigenmodes (λ has positive real part) of LRF should be corrected, as all eigenmodes are expected to decay in a steady system (Kuang 2010). For all the real parts of eigenvalues are negative, and thus the eigenmodes are decaying. However, 8 of the 62 eigenmodes of are growing, and they are likely unphysical. Following Hassanzadeh and Kuang (2016a), we simply set the 8 eigenvalues with positive real part to zero and then reconstruct . Notwithstanding, setting the positive eigenvalues to zero or not does not impact the results substantively in both the forward and inverse applications of the linear response matrix in our analysis. Furthermore, both the calculations of and involve the matrix inversion operation. According to Hassanzadeh and Kuang (2016a), the eigenvalues with smaller magnitude tend to be calculated with higher accuracy, as they are less influenced by the errors through the matrix inversion. Figure 3 shows that the decay rates of eigenmodes λ of are much greater than those of , and thus is more susceptible to the uncertainties in the estimate of eigenvalues and therefore less accurate. As shown later in the next section, indeed fails to estimate the _q_-flux forcing for a given SST response pattern via Eq. (1).

4. Validation of the two LRFs

In this section, we evaluate the performance of the two LRFs by testing 1) whether they can capture the true (modeled) time-mean SST response to _q_-flux forcing (referred to as forward problem), and 2) whether they can give accurate estimates of the _q_-flux forcing for a given SST response pattern of interest (referred to as inverse problem). Both the patch experiments and CPX experiment described in section 2 will serve as the test cases.

a. The q-flux patches

We first evaluate the performance of the two LRFs in solving the forward problem to each localized _q_-flux patch. More specifically, the 106 _q_-flux patches shown in Fig. 1 are first projected onto the selected 62 EOFs, and the SST responses are then directly calculated from δT = −δf. The spatial pattern correlation between the true global SST response and the values linearly constructed through the LRFs is calculated for each _q_-flux patch. The pattern correlation maps shown in Figs. 4a and 4b serve as metrics for the skill of and , respectively.

Fig. 4.

Fig. 4.

Fig. 4.

Pattern correlations between the modeled and constructed SST responses with (a) and (b) for each _q_-flux patch. The color of each dot indicates the pattern correlation coefficient between the constructed and modeled SST response. (c) The local SST response (K) averaged over each patch is shown.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

Here shows limited skill in capturing the SST response to _q_-flux forcing with the pattern correlations (with the modeled pattern) around 0.5 in average. For _q_-flux forcings over the Atlantic Ocean and tropical Indian Ocean, as well as in the 20°S and 20°N bands in the Pacific, the predicted patterns exhibit weak and even negative correlations with the simulated counterparts. A further inspection of the correlation map and the local SST response (Fig. 4c) suggests that they are somehow related. In general, the weaker is the local SST response, the lower is the pattern correlation. Here we select two representative cases in the northern tropical Indian Ocean (centered at 12°N, 60°E) and eastern equatorial Pacific (centered at 1°N, 90°W), regions with the lowest and highest pattern correlations, respectively. In the former case, the constructed SST can hardly capture any remote true responses beyond the regions of forcing (comparing Figs. 5c and 5a, the pattern correlation is −0.45). There are myriad reasons behind this case of failure. One may be related to the weak and locally confined SST response, as the q flux fails to produce strong enough local convection anomalies, and hence wave source for linear teleconnections, and thus the remote response is dominated by the nonlinear or chaotic processes (this can be seen from the very limited stippled regions in Fig. 5a). In contrast, the eastern Pacific case shows rather strong and broad warming near the _q_-flux perturbation, and most of the global ocean is featured with significant linear SST responses (denoted by stippling in Figs. 5b,d,f). As a result, shows sizable skill (pattern correlation is 0.69) in predicting the SST responses for that case. Overall, the accuracy of prediction is poor when the local SST response is weak and confined locally near the _q_-flux forcing. Similar spatial dependence has also been noticed in Fuchs and Sherwood (2016), who used FDT to study the atmospheric responses to a thermal heating profile.

Fig. 5.

Fig. 5.

Fig. 5.

(a),(b) Modeled and predicted SST response to the _q_-flux patch by (c),(d) and (e),(f) . (left) The _q_-flux perturbation centered at (12°N, 60°E) and (right) the q flux at (0°, 150°E). The location of the _q_-flux patch is highlighted by a black ellipse in each panel. In this and the following figures, the crosses denote the regions of linearity (measured by > 1).

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

Figure 4b shows that the skill of is much greater than in replicating the SST response to _q-_flux forcings. The pattern correlations are greater than 0.6 over most of the global ocean, even in regions where fails to capture the sign of the true SST response. In addition, approach appears to capture the magnitude of the SST response better than (not shown), while both tend to overestimate the magnitude. The better skill of the latter is exemplified by the comparison between the predictions by (Fig. 5c; pattern correlation is −0.45) and (Fig. 5e; pattern correlation is 0.87) for the case of the northern tropical Indian Ocean.

b. CPX test case

One might argue that the above comparison is not fair, as is constructed with the patch experiments themselves and should work better if tested against the patch experiments. For a fairer test, we further make use of experiment CPX, which is forced by a complex _q_-flux perturbation, to test the ability of the two LRFs in solving both forward and inverse problems.

In response to the CPX q flux in Fig. 6a, the target SST pattern for the forward problem is presented in Fig. 6b, and the corresponding predictions constructed with and are presented in Figs. 6d and 6f, respectively. The prediction exhibits relatively low skill not only in spatial pattern (pattern correlation is 0.42), but also in the magnitude, with a severe overestimation. On the other hand, the predictions are superior in all aspects in capturing both pattern and magnitude.

Fig. 6.

Fig. 6.

Fig. 6.

The _q_-flux perturbation used to force (a) the CPX experiment and (b) its corresponding SST response; (c),(d) constructed _q_-flux perturbation and SST response with ; (e),(f) constructed _q_-flux perturbation and SST response with ; and (g),(h) constructed _q_-flux perturbation and SST response with .

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

The application to the inverse problem is also evaluated using the same test experiment. The constructed q fluxes calculated via f = −T are compared with the actual q flux shown in Fig. 6a. The completely fails (Fig. 6c), and this may be due to the too-large eigenvalues of compared to , which tends to unrealistically amplify the uncertain eigenmodes in predicting the q flux. By contrast, the prediction with can capture the major feature of the target q flux, although much smoother and with much weaker amplitude. Further analysis suggests the underestimation of the amplitude is attributable to the dimension reduction, as only 45% of the variance in the target _q_-flux forcing is retained after projecting onto the first 62 EOFs. Comparing to the target _q_-flux forcing represented by the retained EOFs (Fig. 7), the prediction (Fig. 6e) actually agrees reasonably well with the projected target q flux. Thus, for a given anomalous SST pattern, is somewhat skillful in retrieving the source of the oceanic forcing.

5. Sources of inaccuracy in

The poor performance of the FDT method may be blamed on (i) the non-Gaussian statistics in SST variability, (ii) insufficient sample size, and (iii) nonnormality of the linear operator extracted from the slab AGCM.

In appendix A, we assess the non-Gaussianity of the leading PCs of the SST from the CTRL simulation. Indeed, 12 PCs of the kept EOF basis, including the leading two, deviate significantly from Gaussian distribution. Insufficient sample size is always a limiting factor for the number of well-resolved EOFs, hence the robustness and the effective spatial degrees of freedom of the LRF. For a reference, Gritsun and Branstator (2007) used 106 day-long data to achieve an accurate atmospheric operator for the daily atmospheric variability. Extrapolating to the much slowly varying SST, simulation of length in the order of 105 years might be needed to reach a similar accuracy for an SST operator.

What further degrades the accuracy of is its nonnormality, because errors are introduced by dimension reduction, as the included and excluded EOFs are strongly coupled through a nonnormal LRF. This is illustrated in appendix B by comparing how the EOF modes interact through a normal versus a nonnormal operator.

To circumvent the problem of EOF interactions, we use Eq. (7) to construct on the model grid without dimension reduction, and the predictions for the inverse and forward problems for the test case CPX are shown in Figs. 6g,h. The prediction skill indeed improves compared to the approach: For the forward problem, the predicted SST pattern shows improved pattern correlation (increasing from 0.69 to 0.89), and for the inverse problem, the predicted q flux captures better the magnitude of the target forcing.

6. Applications

In view of the superior performance of compared to , we hereby apply the LRF derived from the GRF approach to the following applications: (i) global surface temperature sensitivity and (ii) neutral vector for the most predictable SST mode.

a. Global surface temperature sensitivity map

The sensitivity Sk to q flux of any scalar function of the model variables [e.g., temperature, precipitation, and top-of-the-atmosphere (TOA) radiative forcing] is a special case of the Green’s function (Barsugli and Sardeshmukh 2002), and it can be used to evaluate the effectiveness of a _q_-flux forcing at different locations in exciting global changes. An immediate example is the global surface air temperature (TS) sensitivity map, which is calculated as the ratio of the global mean TS anomaly divided by the area-integrated _q_-flux anomaly over each patch:

e8

e8

where angle brackets represent global area-weighted mean. The equation for the sensitivity is very similar to Eq. (7) for the GRF matrix except that the numerator here represents the global mean TS anomalies in response to the _k_th _q_-flux forcing, and thus Sk can be alternatively calculated as the global mean of Gk. Figure 8a shows the sensitivity map Sk, with high sensitivity mainly concentrated in higher latitudes, especially in the Southern Hemisphere south of 30°S. This is consistent with previous studies (Rose et al. 2014; Marshall et al. 2015) that highlight the role of heat uptake over the Southern Ocean in affecting climate response to increased GHG concentrations. By contrast, the tropics overall feature relatively low sensitivity values (roughly one-third of the sensitivity in higher latitudes), especially in the Indian Ocean and Atlantic Ocean. B2006 arrived at a similar sensitivity pattern to the SST forcing in the tropical Pacific. Our sensitivity map shows that the _q_-flux forcing in higher latitudes are 3–4 times more effective in causing global warming or cooling compared to the tropics. To understand this meridional structure in the sensitivity map, we decompose the net TOA radiative forcing into longwave (LW) and shortwave (SW) components from clear and cloudy skies (all defined as positive downward) and composite them for the q flux at each latitudinal band. The resultant sensitivities of these fluxes per unit (PW) q flux from each latitude are shown in Fig. 8b, from which one can see that the enhanced TS sensitivity (Fig. 8b, dashed black line) to high-latitude q flux stems mostly from cloud shortwave feedback (Fig. 8b, magenta line) and clear-sky shortwave feedback (Fig. 8b, red line). While the importance of cloud shortwave feedback in enhancing global TS sensitivity to high-latitude forcing is consistent with the earlier conclusion derived from idealized aquaplanet experiments (e.g., Rose et al. 2014; Rugenstein et al. 2016), the greater importance of the clear-sky shortwave feedback in the large TS sensitivity to Northern Hemisphere high-latitude q flux seems to be a feature in more realistic climate. In-depth investigation into the feedback mechanisms behind the interhemispheric asymmetry will be reported in Part II (Liu et al. 2018, manuscript submitted to J. Climate).

Fig. 8.

Fig. 8.

Fig. 8.

(a) Map of global TS sensitivity (K PW−1). (b) The zonal mean of the global TS sensitivity (dashed black line; K PW−1) and the corresponding radiative feedbacks (W m−2 PW−1): SW clear-sky radiative flux (red line), SW cloud radiative flux (magenta line), LW clear-sky flux (blue line), and LW cloud radiative flux (green line). A spatial smoother derived from the “inpaintn” algorithm (Garcia 2010) is employed here to ensure that only the statistically robust features were retained.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

b. Neutral vector

Neutral vectors are obtained as the right singular vectors of the LRF with the smallest singular values, representing the most excitable patterns of the low-frequency internal variability (see more details in the introduction). Here we will focus the discussion on the first neutral vector of calculated using singular value decomposition (Fig. 9), as it has important bearing on the recent debate regarding the origin of the Pacific decadal variability (Clement et al. 2011).

Fig. 9.

Fig. 9.

Fig. 9.

(a) The first neutral vector (shading) computed as the first right singular vector of , rescaled to have an amplitude of 1 K. Overlaid also are the regression of wind stress (vectors; 10−2 N m−2) and SLP (black contours at intervals of 20 Pa; zero contours are thickened, and negative values are dashed) upon the neutral vector time series; (b) the first optimal forcing computed as the first left singular vector of , and rescaled to have an amplitude of 1 W m−2. (c),(d) As in (a),(b), but for the neutral vector and optimal forcing of derived from a fully coupled CESM1.1.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

The leading neutral vector of displays an IPO-like pattern in the Pacific basin, with horseshoe patterns in midlatitudes and ENSO-like patterns along the equator. Note that large SST anomalies in the tropical Pacific are not symmetric about the equator compared to the canonical IPO; this is common to the simulated SST variability in both CESM–SOM and fully coupled models (e.g., Clement et al. 2011; Okumura 2013). In particular, the warm anomalies tend to be confined to south of the intertropical convergence zone (ITCZ). Because of the northward placement of the ITCZ from the equator, only the internal atmospheric variability of Southern Hemisphere origin can reach the equator and exert significant impact on equatorial eastern Pacific variability through surface turbulent heat flux (Liu and Xie 1994; Xie and Philander 1994; Okumura 2013). As a result, the tropical Pacific tends to be more tightly coupled to the South Pacific than the North Pacific in slab ocean AGCM (e.g., Zhang et al. 2014a,b). In the North Atlantic basin, the neutral vector shows a tripolar pattern resembling the AMO. Since the neutral vector here represents the most excitable mode in a slab model without ocean variability, this SST variability pattern can only be of atmospheric origin, lending some support to the notion that the IPO- and AMO-like variability can exist in a slab ocean under stochastic atmospheric forcing (Clement et al. 2011, 2015). The dominant patterns of the wind stress and sea level pressure (SLP), obtained by regressing onto the neutral vector time series, also correspond well with their observed counterparts (Clement et al. 2011; Okumura 2013; Dong and McPhaden 2017). However, we caution against equating this result to the evidence for the atmospheric origin for the multidecadal oscillation of the SST, as the power spectrum of the first neutral vector shows no significant peak at decadal–multidecadal scales (Fig. 10a, red line). Interestingly, a striking similarity is found between our leading neutral vectors and the most predictable modes of internal variability in Srivastava and DelSole (2017): Our first neutral vector corresponds to their first most predictable mode, and our second and fourth neutral vectors (Fig. 11), characterized by interhemispheric asymmetric pattern and ENSO pattern, closely resemble their second and third most predictable modes, respectively.

Fig. 10.

Fig. 10.

Fig. 10.

(a) The power spectra and (b) autocorrelation functions of the time series resulting from the projection of the CESM–SOM and fully coupled CESM1.1 datasets onto the first neutral vector of CESM–SOM. The red (blue) lines are for the projection of the CESM–SOM (CESM1.1) dataset. The dashed lines in (a) give the 99% confidence intervals for a red noise background.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

To check if the first neutral vector is of any relevance to the most excitable mode in a climate system with active ocean dynamics, we construct using the FDT approach based on a 1000-yr-long control simulation of the fully coupled CESM1.1. The resultant coupled first neutral vector is shown in Fig. 9c, which exhibits considerable spatial similarity to its CESM–SOM counterpart, especially in the IPO- and AMO-like patterns. Therefore, the first neutral vector revealed by our Green’s function experiments with CESM–SOM indeed represents the most excitable mode of SST, irrespective of coupling to the ocean dynamics or not.

The first left singular vector of the LRF, representing the _q_-flux forcing that can optimally excite the corresponding neutral vector, is also estimated from the LRF derived from both the Green’s function experiments with CESM–SOM and the long coupled CESM simulation (Figs. 9b,d). Given the different methods used for the estimate of the LRF,1 the conspicuous difference between the two left singular vectors is no surprise. It is not necessary either that the left singular vector should resemble spatially the corresponding neutral vector, as the SST anomalies are often not locally forced. It is noteworthy that the optimal forcing of the leading neutral vector is characterized by a tripolar pattern in the North Atlantic basin as well, and this is the case for both LRFs from CESM–SOM and fully coupled CESM1.1. Therefore, to the extent that the ocean circulation can organize the tripolar ocean heat divergence/convergence and drive the similarly shaped SST pattern, a corresponding low-frequency dynamically coupled mode is conceivable in the North Atlantic.

To check how ocean dynamics might modulate the temporal characteristics of the leading neutral mode, we project both the 900-yr CESM–SOM and 1000-yr fully coupled datasets onto the neutral vector of CESM–SOM (Fig. 10a) and examine the resultant respective time series. Their power spectra, as well as the autocorrelation functions, are shown in Fig. 10. The spectrum of the CESM–SOM data is overall characteristic of a red noise with an _e_-folding time scale of about 11 months (red lines in Fig. 10). Interestingly, the inclusion of interactive ocean dynamics has little impact on the decorrelation time scale (~12 months), consistent with the findings of Srivastava and DelSole (2017). Yet, coupling to ocean dynamics boosts the power on the interannual time scales with some significant interannual peaks standing out from the background red noise (blue line in Fig. 10a). The oscillatory behavior of the coupled neutral vector can also be discernible from the autocorrelation function in Fig. 10b (blue line). However, no significant decadal–multidecadal periodicities emerge in both spectra.

7. Summary

Motivated by the need to understand the role of the ocean dynamics in the climate response to GHG-induced warming, in this study, we construct the linear response function of SST to the q flux of an AGCM coupled to a slab ocean using both FDT and GRF approaches. The former is predicated on the hypothesis that the SST can be approximated as a Gaussian variable and that its response to external forcing behaves the same way as the variability arising from internal noise, while the latter is achieved through a large set of numerical simulations with patches of q flux specified over the open ocean surface everywhere. The LRF from FDT exhibits some modest skill in capturing the SST response to prescribed q fluxes, especially those from the equatorial Pacific and Southern Ocean, but in general is inferior to the LRF from GRF, epitomizing the challenges in the quantitative understanding of the SST variability and response. In addition to the limited data length used for the construction, suffers also from the non-Gaussianity in the internal variability of the SST and the nonnormality of the LRF itself. Despite sharing the latter two difficulties, the GRF method exhibits considerable skill in both predicting the SST response to a given q flux (forward problem) and retrieving the q flux for a target SST pattern (inverse problem).

The immediate application of the GRF is the sensitivity of the global mean surface temperature to the q flux from different geographic locations. Consistent with an earlier study using an idealized aquaplanet model and idealized q flux (e.g., Rose et al. 2014), we find that high latitudes are 3–4 times more effective in driving global surface temperature change. Somewhat differing from the previous result, we also find an interesting interhemispheric asymmetry in the feedbacks to the amplified global warming sensitivity to q flux from high latitudes: The TOA clear-sky shortwave feedback plays a more important role than the shortwave cloud feedback in response to the NH high-latitude forcing, while the positive cloud feedback is the leading positive feedback for the enhanced global warming response to SH high-latitude forcing. Further detailed analysis using a radiative kernel to dissect the specific feedbacks will be reported in Part II (Liu et al. 2018, manuscript submitted to J. Climate).

The neutral vector analysis of the derived LRF reveals that the SST pattern of IPO corresponds to the leading neutral vector of . Thus, IPO is the most excitable mode of the global SST variability by the stochastic weather noise, corroborating the notion that AGCM coupled to a slab ocean can develop IPO-like variability in the absence of ocean dynamical variability (Clement et al. 2011, 2015). The leading neutral vector is also almost identical to the most predictable SST mode found by Srivastava and DelSole (2017) via average predictability time analysis. Comparison with the neutral vector of a model coupling with ocean dynamics indicates that oceanic feedback can interact with this SST pattern and bring oscillatory behavior to it. But in the CESM modeling system examined, this only occurs on interannual time scales, and we find no significant periodicities at decadal–multidecadal range in the spectra of both the neutral vectors. Thus, one should refrain from extrapolating this result to longer time scales to argue for an oceanic origin for the multidecadal oscillation of the IPO. In Part II (Liu et al. 2018, manuscript submitted to J. Climate), evidence will be presented to show that neutral modes can organize the global radiative feedback to give rise to distinct global temperature response despite the same _q_-flux forcing magnitude.

Owing to the limited scope of this paper, only a few issues related to the SST LRF are examined here. As the GRF experiments simulate all the atmospheric, land, and ice variables, one, in principle, can build LRF for any variable of interest. Research is underway to explore the sensitivities of the ITCZ, jet location and strength, and Hadley cell width and intensity to the q flux, so more results will be forthcoming to address a broader set of questions of climate sensitivities and feedbacks.

Acknowledgments

We thank Professor F.-F. Jin and an anonymous reviewer for their very constructive reviews that helped improve our manuscript substantively. The CESM–SOM warm and cold _q_-flux patch simulation output are archived at the National Energy Research Scientific Computing Center (NERSC) and can be accessed by contacting Jian Lu (jian.lu@pnnl.gov). This study is supported by the U.S. Department of Energy Office of Science Biological and Environmental Research (BER) as part of the Regional and Global Climate Modeling program. Computational resources for the warm and cold _q_-flux patch simulations were provided by NERSC. The Pacific Northwest National Laboratory (PNNL) is operated for DOE by Battelle Memorial Institute under Contract DE-AC05-75RL01830. F. Liu is supported by a Chinese Scholarship Council visiting student fellowship. X. Wan is supported by NSFC (41576004) and the National Basic Research Program of China (2014CB745001).

APPENDIX A

Gaussianity Test

We assess the Gaussianity of the state vectors by checking their skewness and kurtosis (White 1980; Sura and Sardeshmukh 2008). The skewness, which measures the asymmetry of a probability distribution function (PDF), is calculated as , where the overbar denotes the time mean, and T is the time series of the first 62 PCs of the global SST. The kurtosis is defined as , and it describes the “peakiness” of the PDF. To determine if the calculated skewness and kurtosis are significantly different from those of a Gaussian distribution (whose skewness = 0 and kurtosis = 3), we follow the criteria of Brooks and Carruthers (1953) that a distribution is construed significantly different from Gaussian when its skewness or kurtosis are more than twice the corresponding standard errors. The standard errors of skewness and kurtosis are estimated as and , respectively, with the effective sample size _N_eff evaluated using Eq. (5). Fig. A1 shows the scatterplot between kurtosis and skewness for the first 62 leading PCs of the SST data from CESM–SOM. Twelve of them, including the first two leading PCs, are significantly non-Gaussian. Interestingly, the non-Gaussianity is more often caused by kurtosis than skewness: Among all the 12 non-Gaussian PCs, only 2 PCs deviate from Gaussian as a result of large skewness.

Fig. A1.

Fig. A1.

Fig. A1.

Scatterplot of kurtosis vs skewness for the first 62 leading PCs of the global SST. Stars represent the PCs statistically different from a Gaussian distribution at 95% confidence level. The ranks of the PCs are color coded according to the color bar on the right.

Citation: Journal of Climate 31, 9; 10.1175/JCLI-D-17-0462.1

APPENDIX B

Nonnormality of LRF

To illustrate how dimension reduction can introduce error if the included and excluded EOFs are strongly coupled through a nonnormal LRF, a linear system dT/dt =

T + f is projected onto n EOF basis:

eb1

eb1

where a = (_a_1, _a_2, …, a k) is a coefficient vector corresponding to the first _k_th EOFs, Λ is a diagonal matrix containing the eigenvalues λ of LRF,

contains the EOFs of SST,

contains the eigenvectors of LRF, and fe is the projection of f on the EOF basis. For a system governed by a normal LRF,

and

are identical, and the first term on the right-hand side of Eq. (B1) becomes Λa, and the corresponding system becomes uncoupled and can be expressed as

eb2

eb2

If only the first _k_th EOFs are retained after dimension reduction, error can only be attributable to the omission of the projection of f onto the excluded EOFs, but the first _k_th modes would not be affected because of the orthogonality.

On the other hand, if LRF is nonnormal, and hence

is different from

, the matrix in the square brackets in Eq. (B1) will not be diagonal, but have large off-diagonal elements β ij. The corresponding equation for an becomes

eb3

eb3

Consequently, all the EOF coefficients are coupled to one another, and when truncated to the first _k_th EOFs, the equation for the first _k_th coefficients becomes

eb4

eb4

Comparing Eq. (B4) with Eq. (B3), one can see that the accuracy of all the retained EOFs is compromised by not accounting for the coupling with the excluded EOFs:

. Both

and

constructed in our paper show considerable nonnormality as delineated by the off-diagonal elements in Fig. B1.

REFERENCES