## Abstract

Transfer function tools commonly used in engineering control analysis can be used to better understand the dynamics of El Niño Southern Oscillation (ENSO), compare data with models and identify systematic model errors. The transfer function describes the frequency-dependent input–output relationship between any pair of causally related variables, and can be estimated from time series. This can be used first to assess whether the underlying relationship is or is not frequency dependent, and if so, to diagnose the underlying differential equations that relate the variables, and hence describe the dynamics of individual subsystem processes relevant to ENSO. Estimating process parameters allows the identification of compensating model errors that may lead to a seemingly realistic simulation in spite of incorrect model physics. This tool is applied here to the TAO array ocean data, the GFDL-CM2.1 and CCSM4 general circulation models, and to the Cane–Zebiak ENSO model. The delayed oscillator description is used to motivate a few relevant processes involved in the dynamics, although any other ENSO mechanism could be used instead. We identify several differences in the processes between the models and data that may be useful for model improvement. The transfer function methodology is also useful in understanding the dynamics and evaluating models of other climate processes.

## 1. Introduction

The dynamics of the El Niño Southern Oscillation (ENSO) are not consistently well captured in general circulation models (GCMs) [1–6]. Furthermore, even if the output statistics (e.g. NINO3 period, amplitude, etc.) are reasonably well simulated, there is no guarantee that there are no compensating errors in the model that lead to seemingly good results for the wrong reason (as noted, for example, in [6]). As a result, even if the present climate is well simulated, then there is no guarantee that changes in the dynamics with climate change will be correctly captured owing to systematic errors not compensating anymore in the changed climate. It is thus important to evaluate models at an individual process level, in addition to matching general output metrics such as those based on NINO3 time series, for example. The same issue of compensating model errors is relevant in any aspect of climate dynamics, e.g. for understanding multi-decadal variability of ENSO, variability of the Atlantic meridional overturning circulation (AMOC), and more. This issue is not restricted to climate variability modes, and models that fit the observed global warming record of the twentieth century diverge in both prediction mode [7] and in trying to simulate AMOC during the last glacial maximum [8]. This again indicates the possibility of compensating errors and requires model validation by process rather than by output metrics such as the globally averaged surface temperature.

We applied the *transfer function* as an additional tool to estimate the process dynamics of ENSO [9] and to analyse AMOC variability [10]. Here, we (i) provide more detail on the transfer function estimation and error bounds, as well as guidelines on its use, and (ii) extend the previous results to consider additional processes involved in ENSO dynamics motivated by the delayed oscillator. In addition to knowledge regarding the processes considered, these two goals illustrate different features of transfer function analysis. In a related work, we have also considered a number of processes based on the alternative recharge-oscillator ENSO description [11].

The term ‘transfer function’ is taken from the control engineering literature [12], and it refers to the linear frequency-dependent input–output relationship between any pair of causally related variables. This is equivalent to finding the output response to an input that is composed of a time series of a single frequency, as a function of the input frequency. The transfer function can be derived from the underlying differential equation relating the input and output variables, or, as here, estimated from data or model time series of the input and output variables. First, this can assess whether the underlying relationship is or is not frequency-dependent—this should be an essential first step in any analysis of the relationship between any two variables, before, for example, least-squares fitting (regression) of the relationship. Further, even if there is no significant frequency-dependence, the transfer function can be used to understand the frequency range over which the data can meaningfully describe an underlying relationship. Finally, if the relationship is frequency-dependent, then the transfer function can be used both to estimate or verify the underlying differential equation and to estimate values of unknown parameters.

In MacMynowski & Tziperman [9], we considered two examples. First, the western boundary reflection process relating westward-propagating Rossby waves and the reflected eastward-propagating Kelvin waves. This relationship was considered as there are ample existing estimates of the reflection coefficient to compare our results with. Second, the relationship between arriving Kelvin waves in the eastern Pacific, and the resulting eastern Pacific sea surface temperature (SST) anomaly (NINO3 index). These relationships were compared for TAO data from ocean buoys [13], model output from the GFDL CM2.1 coupled GCM [14,2], and the Cane–Zebiak model of ENSO [15]. The analysis suggested compensating errors in the GFDL model, and clearly indicated large errors in the more idealized Cane–Zebiak model (see also [16]). Linz *et al.* [11] consider additional ENSO processes motivated by the recharge oscillator description [17] in five CMIP5 models (the Coupled Model Intercomparison Project Phase 5, [18]).

Here, we return to the delayed-oscillator ENSO description; this and Linz *et al.* [11] thus provide complementary descriptions of ENSO dynamics. The two input–output relationships already considered in MacMynowski & Tziperman [9] illustrate the utility of the transfer function analysis, and the importance of estimating process dynamics, but they do not give a complete picture of ENSO dynamics. Here, we present the transfer functions for a total of six processes which span the entire delayed oscillator feedback loop [19,20] shown in figure 1. Each block in this figure represents a process (i.e. relationship between two variables) which will be analysed from appropriate model and data time series in §3. We present and analyse the six transfer functions for TAO array observations and several models, including GFDL CM2.1 as in MacMynowski & Tziperman [9] but using significantly more data, and also CCSM4 [21], and the Cane–Zebiak model.

The transfer function approach represents a tool that may be useful in diagnosing process dynamics, and which complements other approaches for analysing individual physical processes; see a summary in Guilyardi *et al.* [22], Collins *et al.* [23] and other examples [24–31].

Transfer functions are also expected to be useful in identifying compensating model errors in studies of other climate variability modes and climate change and we have made the required code available on our web pages. Section 2 introduces the transfer function analysis tools, and the application of these tools to ENSO is discussed in §3.

## 2. Transfer function estimates and error analysis

We start with a brief introduction to transfer functions following MacMynowski & Tziperman [9], with an example to better illustrate its use, and additional information provided on the error analysis in particular.

A key advantage of transfer functions is that their calculation provides information on the differential equation describing the dynamics of the input–output relationship, even when it is not known *a priori*. It is not unusual to see the application of analyses such as regression methods, that ignore a possible frequency dependence, without an explicit justification. Consider for illustration the following heuristic equation for the NINO3 temperature *T* as a function of the east Pacific Kelvin wave amplitude *K*_{e} [17],
*μ* represents the effects of east Pacific Kelvin waves (*K*_{e}) on the SST (*T*) via both the thermocline feedback [32] and advective feedback [25] and *ϵ* represents dissipation processes. By taking the Fourier transform with a frequency *f*, we have
*s*=2*πif* and *Φ*( *f*), provide useful information,
*b*).

The frequency-dependent transfer function can also be estimated from time series of the input and output without knowledge of a corresponding differential equation ([33 section 6.2]). Given input and output time series *x*(*t*) and *y*(*t*) and their Fourier transforms *S*_{xx}( *f*) and *S*_{xy}( *f*) can either be estimated using multi-taper methods [34] or through periodograms [35]. Here, we use the latter and (i) divide the time series for *x* and *y* into *n* possibly overlapping segments of smaller length, *x*_{k}(*t*) and *y*_{k}(*t*); (ii) compute the Fourier transforms *x*_{k}(*t*), *y*_{k}(*t*) in the *k*th segment, and calculate their products in each segment and (iii) average these products over the segments to eliminate the uncorrelated parts of the response,

In general, any variable in the climate system may depend on many different inputs. With sufficient averaging, the contribution owing to any part of the output signal that is not correlated with the input signal will tend to zero, yielding an estimate due only to the correlated component. However, increasing the averaging using a finite length time series by increasing the number of segments results in shorter segments and thus smaller resolved frequency range. A compromise we use here is to estimate the transfer function at different frequencies using different segment lengths, thus the lowest frequency in our plots uses fewer averages, whereas higher-frequency transfer function estimates are based on shorter segments and more averaging.

Because ‘windowing’ (gradual tapering of the time series to zero near their ends) is required on each data segment to avoid the implicit discontinuity at the segment boundary, non-zero overlaps between segments can be used to increase the number of averages without creating information that was not present in the original data; see Welch [35] for a more detailed discussion. We apply a Hamming window on each segment, and use a segment overlap of 50% (e.g. using 2 year data segments would correspond to using years 0–2, 1–3, 2–4, etc., with a lowest resolved frequencyof 1/(2 years)).

The magnitude of the error in *T*_{xy}( *f*) that results from contributions to the output signal not correlated with the input signal can be estimated from the coherence,
*d*). In practice, there are rarely sufficient data to make these errors negligible.

Nonlinear effects in the relationship between input and output variables will also reduce coherence; information at one frequency in an input *x*(*t*) would in this case also result in signal content at a different frequency in the output *y*(*t*), in general uncorrelated with the signal content in *x* at that frequency. In many cases, it is difficult to distinguish between poor coherence caused by nonlinearity, or by some other variable affecting the output. If the true relationship is nonlinear, then the presumed linear relationship will also depend on the signal amplitude. The estimated value will not in general yield the linearization of the relationship about *x*=0, but will be closer to the average slope of the input–output relationship at the typical amplitudes of the input signal in the data. Much greater detail on the effects of nonlinearity can be found, for example, in Schoukens *et al.* [36] and Nordsjo & Wigren [37].

There are additional sources of uncertainty in the interpretation of the estimated transfer function. A correlation between *x* and *y* may be due to the input *x* affecting the output *y* (the presumed interpretation), owing to a closed feedback loop where the output *y* affects the input *x*, or owing to some independent factor affecting both *x* and *y*. As an example, consider western Pacific Rossby waves influencing the reflected Kelvin waves. A closed feedback loop effect exists, because the reflected Kelvin waves, in turn, affect east Pacific SST, the SST influences winds, which in turn excite Rossby waves again. An independent factor affecting both input and output could be western Pacific wind anomalies that simultaneously excite both Kelvin and Rossby waves, biasing the estimate of the western boundary reflection coefficient. The phase of the transfer function can be useful in diagnosing these factors, by noting which signal leads the other, or whether the two signals are in phase.

As an example of critically evaluating a transfer function estimate, consider the following system, modified from (2.1) by the addition of a ‘leakage’ path where the output *y*(*t*) is directly influenced by the input *x*(*t*), and also influenced by an uncorrelated noise *ν*,
*μ*=1 and *ϵ*=4 (years^{−1}). Assume that the input *x* is a known random Markov process with an auto-correlation time four times larger than the time constant *ϵ*^{−1}, and that *ν* is white noise, which therefore dominates the response at high frequencies. The term *γx* in the output corrupts the estimate of the desired process. The transfer function is estimated from 20 years of simulated data, comparable to the TAO record. Figure 2 shows the magnitude and phase estimates (with error bars indicating two standard deviations), as well as the frequency-dependent coherence. Two different cases are plotted, corresponding to *γ*=0.1 and *γ*=0. In each case, the lowest frequency is estimated based on averaging over 4 year segments, and higher frequencies from 2 year, 1 year or six month segments.

As an example of how the transfer function phase may be used to critically evaluate the results, note that with *γ*=0 the phase corresponding to (2.10) should be −90° at high frequencies. This is because our equations reduce in the absence of noise to *f*≫*ϵ*, this may be approximated as *i*=*e*^{π/2} implying a 90° phase difference between the input *x* and output *y*. However, introducing the additional pathway from input-to-output results in an average phase closer to 0°, indicating that the output at these frequencies is not dominated by the dynamics governed by *μ* and *ϵ* in (2.10). Thus, the transfer function phase indicates that because of the leakage path, the high-frequency information is not useful in understanding the dynamics associated with (2.10) at lower frequencies. This is clear despite the high uncertainty at high frequencies (evident in the coherence) owing to noise *ν*.

This example demonstrates how estimating the transfer function at each frequency is immediately useful in understanding the underlying differential equation that relates the variables. The transfer function estimates can also be used to estimate the parameters of the differential equation, e.g. *μ* and *ϵ* in (2.10). In principle, these parameters could be identified directly from the time domain data *x*(*t*) and *y*(*t*) either using least-squares in the case of fitting a constant relationship, or by fitting an AR model. However, if some fraction of the frequency range is dominated by noise or by other processes (as in the above example at high frequencies), then this would bias the time-domain estimate. Such a problem should be immediately apparent in the transfer function frequency response, and can then be taken into account either by choosing which frequency range to use for estimating parameters, or by filtering the time-domain data before fitting parameters of an AR model.

The fit of model parameters such as *μ* and *ϵ* in the above example also requires careful treatment. Consider first the case where the assumed transfer function is a constant *α* over frequency, and we are trying to fit its value. Then, the best estimate of the constant *α* can be obtained from a weighted least-squares, using the variance (2.8) to weight the quality of the transfer function information *T*_{xy}( *f*) at each frequency,

The more general case of a frequency-dependent fit is illustrated by the transfer function (2.3). Fitting the squared magnitude of the transfer function, the data at all *N* frequencies are related to the unknown parameters *μ* and *ϵ* by,
*μ* and *ϵ*, where each measurement (information at each frequency *f*_{i}) is weighted by the inverse of its variance. This algorithm is used to obtain the fits in the following analysis of ENSO models and observations.

## 3. Application to El Niño Southern Oscillation

### (a) Choice of transfer functions to evaluate

The first step in using transfer functions to diagnose the dynamics is to define appropriate input–output variable pairs. Returning now to the schematic representation in figure 1, define the following variables.

—

*T*as the SST averaged over the NINO3 region,—

*W*_{r}and*W*_{k}as the central Pacific atmospheric wind stress projected onto the Kelvin and Rossby meridional structure (as in [38]),—

*K*_{c},*R*_{c}as the central Pacific Kelvin and Rossby waves excited by wind anomalies,—

*K*_{w},*R*_{w}as the western Pacific Kelvin and Rossby wave amplitudes, and—

*K*_{e}as the eastern Pacific Kelvin wave amplitude.

Figure 1 can be used to motivate the following delay differential equation describing the dynamics of ENSO [19,20],

The projection of data onto Kelvin and Rossby waves is done following Boulanger & Menkes [38], Boulanger *et al.* [39], using a single baroclinic mode based only on the layer thickness anomaly (and not currents) [39]. The atmospheric wind forcing on the Kelvin and Rossby waves is computed as in Boulanger & Menkes ([38, equation A 7]). For TAO data, the wind stress is computed using the expression for friction velocity derived by Vera (unpublished manuscript, 1983) and published as Large *et al.* ([40, equation 8]).

The appropriate longitude ranges of ‘east’, ‘central’ and ‘west’ (figure 1) can be chosen through correlation analysis, shown in figure 3 for the TAO array data, weekly averaged GFDL CM2.1 model output, and CZ model output. Similar correlation plots have been used by Zelle *et al.* [41] to estimate time delays.

The left columns of figure 3 show the maximum influence of Kelvin waves on NINO3 to occur east of 240°, whereas the right columns of figure 3 show most of the influence of NINO3 on creating Kelvin waves to be west of 240°. Hence this longitude can be used to understand the influence of Kelvin waves on SST without a confounding feedback influence of SST on Kelvin waves. The resulting transfer function is found to be not strongly dependent on the chosen longitude.

Similar correlation plots indicate that most of the influence of SST on the forcing of ocean wave modes occurs east of 156°. Thus, this longitude can be used to evaluate the Rossby and Kelvin modes near the western boundary in order to calculate the reflection coefficient *T*_{RK} without interference from the forcing of the waves by the atmospheric response to eastern Pacific SST. This longitude also has the advantage of being sufficiently removed from land masses close to the equator.

Finally, the peak excitation of winds by NINO3 is near 180° (in the GFDL model) and 190° (in the TAO data and CZ model), and the latter is therefore used in evaluating the corresponding transfer functions *T*_{TWk,r} and *T*_{WrR}, *T*_{WkK}. In order to estimate the effect of wind on the oceanic Kelvin wave directly, without a confounding influence from wind causing Rossby waves that reflect off the western boundary as Kelvin waves, it is necessary to subtract a delayed version of the western Pacific Kelvin wave from the central Pacific Kelvin before computing the transfer function *T*_{WkK}.

### (b) Results for El Niño Southern Oscillation transfer function estimates

Transfer function estimates for each of the six processes represented by the six boxes in figure 1 are computed from TAO data, GFDL CM2.1 and CCSM4 GCM output, and CZ model output, and plotted in figures 4–6. As implied by figure 1, the analysis assumes that each process relates scalar input and output signals described by a (forced) linear differential equation. This is an approximation both because the relationship might not be linear, and perhaps more significantly because the relevant state of the ocean at each point in the conceptual diagram of figure 1 cannot be captured by a single scalar variable. Nonetheless, much insight can be gained through these simplifications; indeed, the latter assumption is also implicit in simplified delayed-oscillator or recharge oscillator descriptions of ENSO.

The western boundary reflection coefficient and *T*_{KT} (panels *a* and *b* of figures 4–6) were discussed in MacMynowski & Tziperman [9]; calculations here with GFDL CM2.1 are updated to use 500 years of model output rather than only 22 years. Transfer function estimates, and the correlation estimates in figure 3, are based on weekly averaged TAO data from 1994 to 2013, 500 years of monthly averaged GFDL CM2.1 and CCSM4 model output simulated under pre-industrial conditions, and 75 years of CZ model output, all with the seasonal cycle removed. The GFDL model captures the general features of ENSO reasonably well but with too high an amplitude and possibly too short a period [14]. Because the time scale of the processes estimated here is short (weeks to a few months), the roughly 20 years of data available are sufficient to reasonably constrain these transfer functions even if calculating ENSO event statistics generally requires a longer record [42].

Consider first the TAO array estimates. Because of the relatively short data record, we use 2 years as the longest data segment for estimating transfer functions that depend on frequency only weakly, and a 4 year maximum segment length for frequency-dependent transfer functions in order to better resolve low-frequency characteristics. The western boundary reflection coefficient is reasonably constant at low frequencies; the higher frequency content in figure 4*a* (and figure 5*a*) is due to correlated excitation of Kelvin and Rossby waves from wind disturbances [9]. This can be verified from the roughly constant phase at higher frequencies (figure 7), rather than from the linearly increasing phase with frequency expected from a time delay. The magnitude of the reflection coefficient is consistent with other estimates, e.g. 0.33–0.37 in Boulanger *et al.* [39]. Also described in MacMynowski & Tziperman [9], the excitation of SST owing to eastern Pacific Kelvin wave anomalies is consistent with (2.1), validating the form of the underlying dynamics for this process. Fitting the transfer function yields estimates for both the feedback strength *μ* in (2.1) and the dissipation *ϵ*; these are indicated in the figure.

The transfer function describing the forcing of wind anomalies projected on the Rossby and Kelvin modes by East Pacific SST, *T*_{TWk,r}, are given in panels (*c*) and (*d*) of figures 4–6 and show that NINO3 anomalies produce a roughly proportional wind–stress anomaly (i.e. independent of frequency). This is consistent with the short atmospheric response time of a few days.

However, an interesting result we obtain here is that the forcing of ocean waves from applied wind perturbations (*T*_{WkK} and *T*_{WrR}, panels (*e*) and (*f*)) does involve some frequency dependence, indicating a prognostic relationship between the two (i.e. they are related through an expression involving a time derivative) which is normally ignored in the ENSO literature. The system responds more strongly to low-frequency perturbations than to high frequencies. This frequency dependence suggests (see examples in §2) a first-order differential equation of the form *τ* d*R*/d*t*=*gW*_{r}−*R* and therefore a corresponding transfer function *T*_{WrR}=*g*(1+*τs*)^{−1}. The time constant *τ* weighting the time derivative term for generating Rossby waves is estimated to be approximately five months and is fast enough that it may be a reasonable approximation to ignore it in simple heuristic models of ENSO such as the delayed oscillator. However, based on the transfer function shape, we can predict that the low-frequency relationship between the wind and the forced Rossby waves will be underestimated if the time series (involving the full frequency content of the signal) is used to naively estimate the value of an assumed frequency-independent relationship between the wind and the forced Rossby waves. The time constant involved in generating Kelvin waves is longer and difficult to reliably estimate from this data (and also from GFDL or CCSM model output); the fit indicated in the figure only indicates the integral part of the relationship, d*K*/d*t*=*αW*_{k} and transfer function *T*_{WkK}=*α*/*s*.

Next, consider the two GCMs analysed here, which have roughly similar characteristics for all six transfer functions considered here. With longer model time series, the transfer functions are more readily extended to lower frequencies (with monthly averaged data, 1/(two months) is the highest frequency resolved). Both the GFDL and CCSM model appear to have slightly too high a reflection coefficient; this may indicate other processes involved in creating low-frequency correlations between the Kelvin and Rossby wave (whether the latter is causally forced by the Kelvin wave or by other processes that influence both variables); in any case, this is not expected to be dominant in ENSO dynamics [39]. Note that the value of the reflection coefficient differs from that estimated in MacMynowski & Tziperman [9] owing to the inclusion of many more years of model output here. Small errors in the excitation of NINO3 perturbations from Kelvin wave anomalies in both the strength of the response and the dissipation are not likely significant relative to the errors in estimating transfer functions from the short TAO data record in particular.

Additional compensating errors are found in the transfer functions considered here for the first time: we find that both the GFDL and CCSM4 models simulate too weak a wind stress response to NINO3 anomalies (transfer functions *T*_{TWk} and *T*_{TWr}; figure 5*c*,*d*). This is again compensated by another error of the model wind stress having too large an effect on ocean waves (transfer functions *T*_{WkK} and *T*_{WrR}; figure 5*e*,*f*), particularly for the excitation of Kelvin waves in both models. The inertia effect in wave excitation (i.e. the time constant *τ* in *T*_{WrR} in particular) might be too strong, with high-frequency perturbations resulting in a smaller ocean response than the corresponding process estimated from data; however, the uncertainty in estimating this parameter is substantial.

None of these process discrepancies corresponds directly to individual parameters that can be modified in the GCM. However, understanding errors in individual processes is still invaluable in tuning the model to improve the simulation of ENSO.

The CZ model western boundary reflection coefficient is specified in the model, and MacMynowski & Tziperman [9] showed that the transfer function analysis successfully yields close to the correct value of a theoretical perfect reflection of 0.41 (figure 6*a*). The NINO3 response to Kelvin wave (thermocline) anomalies in the CZ model is more than a factor of two too large. The ability to tune parameters in this model to improve the match between model and data is illustrated in MacMynowski & Tziperman [9], where both the advective and thermocline feedbacks are re-tuned in the model to improve the match of this process transfer function with data.

We now proceed to other systematic CZ model errors that are revealed by the additional transfer functions calculated here. The excitation of wind stress by NINO3 anomalies (transfer functions *T*_{TWr} and *T*_{TWk}) in the CZ model is 50% larger than in observations. The excitation of Rossby waves by the winds (*T*_{WrR}, compare figures 6*f* and 4*f*) has both an incorrect frequency dependence and again too strong an amplitude at the relevant low frequencies. There is also a significant error in the frequency dependence of the Kelvin wave response to wind stress anomalies, though the low-frequency amplitude is similar to the observed (*T*_{WkK}, figures 6*e* and 4*e*). Overall, the product of the three transfer functions involved in the East Pacific Kelvin wave feedback (*T*_{KT}×*T*_{TWk}×*T*_{WkK}) together corresponding to the feedback parameter *α* in (3.1) is a factor of 3 larger than in TAO data. This is a significant difference, indicating the need for observationally motivated model tuning.

Finally, the phase of the transfer functions also provides useful information, often allowing us to evaluate the reliability of the amplitude estimate. An example is shown in figure 7 for the phase of the transfer function between western Pacific Rossby waves and Kelvin waves, corresponding to the western boundary reflection coefficient. The waves are estimated some distance away from the boundary and therefore one expects some delay between the arriving Rossby waves and the reflected Kelvin waves. In general, the input–output relationship may therefore be written as *K*(*t*)=*a*×*R*(*t*−*τ*). Taking the Fourier transform, the transfer function is *T*_{xy}(*s*)=*a* *e*^{−i2πfτ} which has an amplitude equal to the reflection coefficient *a*, and a phase increasing linearly with frequency, *ϕ*_{xy}( *f*)=−2*πτf*. This clearly fits the transfer function from Rossby to Kelvin wave evaluated from the CZ model in figure 7*c*, with a time-delay of two months corresponding to the wave propagation delay from 156°W where the waves are estimated, to the western boundary (124° in this model) and back. The frequency axis is plotted on a log scale, hence the phase does not show as a straight line. As noted earlier, the corresponding transfer functions evaluated from either TAO data or the GFDL model output show a constant phase at higher frequencies (figure 7*a*,*b*), indicating that there is no time delay between the signals at these frequencies, and suggesting that something other than wave reflection is involved. One possibility is a simultaneous excitation of western Pacific Rossby and Kelvin waves by local wind perturbations, leading to the waves having the same phase. The observational record is not sufficiently long to confidently determine the time delay at lower frequencies owing to wave travel from the longitude where the waves are estimated to the western boundary and back. The apparently larger delay in observations than in the CZ model may be due to the dominant reflection occurring further west than the CZ model boundary, differences in wave propagation speed, other physical processes that result in Kelvin waves or simply insufficient data for a reliable estimate.

## 4. Conclusion

Transfer functions are a useful tool for understanding the dynamics of processes involved in ENSO, or potentially other aspects of climate variability. Using time series of input and output signals, we can identify frequency–domain relationships between them, and use this both to understand the underlying differential equations describing the dynamics, and to estimate parameters. Applying this to estimates of processes involved in ENSO using TAO data, the GFDL CM2.1 and CCSM4 GCM models, and the CZ model, we are able to quantitatively compare different individual processes between the observations and models, and therefore identify several significant compensating errors in both models.

For both the GFDL-CM2.1 and the CCSM4 models, we find differences between the model and observations indicating compensating model errors in processes involving the wind stress. The wind stress induced by SST changes is too weak, and this is compensated by too large an ocean wave response to a given wind stress perturbation. In addition, the response time involved in the ocean response to applied wind stress may be too large in the models, although this is more difficult to reliably estimate. These identified model errors are in addition to compensating errors in the SST response to Kelvin wave (thermocline) anomalies that were found for the GFDL model by MacMynowski & Tziperman [9].

Using the original values for all CZ model parameters, the process errors in this more idealized model are found to be larger than in the GCMs examined here. The SST response to an arriving Kelvin wave is too strong, the wind stress anomalies that result from SST anomalies are also too strong, and these are partially compensated by a too weak ocean response to wind stress anomalies. These process errors suggest that the seemingly realistic simulation of ENSO using the original model parameters is likely influenced by compensating errors.

The first and perhaps most crucial step in a transfer function analysis is the choice of input and output variables. Our choices were motivated by the delayed-oscillator mechanism. The advantage of the methodology is that if the choice is wrong (say because this mechanism has nothing to do with the dynamics of ENSO), then one obtains an immediate indication of this in the results, as the transfer function will not be significantly different from zero. Thus, the method can not only be used to quantify individual processes, but also to find out if a given mechanism is consistent with the model output or observations. We find the transfer functions spanning the entire delayed-oscillator ENSO mechanism are significant, and thus that this mechanism is a viable description of the ENSO cycle. This does not exclude additional effects or alternative descriptions such as the recharge oscillator [17] or others [43,44]; in companion work, we also examine transfer functions corresponding to the recharge description of ENSO [11].

This analysis of ENSO models and observations here demonstrate that transfer functions can be used to quantify individual physical processes in climate models in a way that allows the identification of model errors and compensating errors. Given that such compensating errors are a major concern for the reliability of model future predictions, we expect that this tool can be used to study and improve climate models in general.

## Funding statement

Funding for this work is from DOE grant no. DE-SC0004984 and NOAA awards NA13OAR4310129 and NA13OAR4310130. E.T. is also supported by the NSF climate dynamics programme, grant no. ATM-0754332 and by the NASA ECCO-II project.

## Acknowledgements

The ocean data were made available by the TAO Project Office of NOAA/PMEL. E.T. thanks the Weizmann Institute for its hospitality during parts of this work.

## Footnotes

One contribution to a Special feature ‘Climate dynamics: multiple scales, memory effects and control perspectives’.

- Received April 3, 2014.
- Accepted May 28, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.