Abstract

The parameter estimation problem for the coupled ocean–atmosphere system in the tropical Pacific Ocean is investigated using an advanced sequential estimator [i.e., the extended Kalman filter (EKF)]. The intermediate coupled model (ICM) used in this paper consists of a prognostic upper-ocean model and a diagnostic atmospheric model. Model errors arise from the uncertainty in atmospheric wind stress. First, the state and parameters are estimated in an identical-twin framework, based on incomplete and inaccurate observations of the model state. Two parameters are estimated by including them into an augmented state vector. Model-generated oceanic datasets are assimilated to produce a time-continuous, dynamically consistent description of the model’s El Niño–Southern Oscillation (ENSO). State estimation without correcting erroneous parameter values still permits recovering the true state to a certain extent, depending on the quality and accuracy of the observations and the size of the discrepancy in the parameters. Estimating both state and parameter values simultaneously, though, produces much better results. Next, real sea surface temperatures observations from the tropical Pacific are assimilated for a 30-yr period (1975–2004). Estimating both the state and parameters by the EKF method helps to track the observations better, even when the ICM is not capable of simulating all the details of the observed state. Furthermore, unobserved ocean variables, such as zonal currents, are improved when model parameters are estimated. A key advantage of using this augmented-state approach is that the incremental cost of applying the EKF to joint state and parameter estimation is small relative to the cost of state estimation alone. A similar approach generalizes various reduced-state approximations of the EKF and could improve simulations and forecasts using large, realistic models.

1. Introduction

This is Part II of a two-part study that deals with state and parameter estimation for models of the ocean–atmosphere system in the tropical Pacific Ocean. In the first part of this study, Sun et al. (2002, hereinafter Part I) successfully estimated the model state of an intermediate coupled model (ICM) of this system using an extended Kalman filter (EKF). However, studies on many coupled tropical ocean–atmosphere general circulation models (GCMs; Neelin et al. 1992) have shown that there can be substantial discrepancies between simulated and observed climatology, as well as between the results of a coupled model and the same ocean model forced with observed winds. Incorrect parameter values used in the models are a likely cause for these systematic errors.

The behavior of coupled GCMs beyond the time scales of numerical weather prediction is determined largely by the model’s parameterizations, rather than by the initial state. Predictions of future global climate change, for example, yield a wide range of changes because of large uncertainties in model parameters (Murphy et al. 2004).

Estimating model parameters has always been an important part of the modeling enterprise (Ljung 1987). Successful parameter estimation can be achieved provided that the system is sensitive to the parameters, and state observations are sufficient. The former requirement is related to the concept of identifiability, that is, whether the identification procedure will yield unique values of the parameters and whether the resulting model is identical to the true system that generates the observations. In other words, are the observations informative enough to distinguish between a correct and an incorrect model parameter value (Navon 1997)?

Consequently, the problem of parameter estimation has attracted considerable attention in the meteorological and oceanographic literature (Ghil and Malanotte-Rizzoli 1991; Ghil 1997; Navon 1997). However, to implement systematic parameter estimation, rather than the usual trial-and-error procedures in state-of-the-art coupled GCMs, is still a considerable challenge, and there have been few attempts to do so with observational data.

Smedstad and O’Brien (1991) have used an adjoint method to estimate the wave speed in a linear tropical ocean model. Adjoint methods, however, cannot readily estimate the accuracy of both parameter and state. Annan et al. (2005) applied the ensemble Kalman filter (EnKF; Evenson 2003) to tune the climatology of an ICM by estimating 12 scalar model parameters in identical-twin experiments. The accuracy of the covariance approximated from a limited EnKF ensemble depends, however, on the number of realizations used and several other factors, especially the sophistication of the choice of basis for the ensemble. The high computational burden of state-of-the-art coupled GCMs, though, puts severe constraints on the ensemble size (Keppenne 2000; Keppenne and Rienecker 2002).

Interannual variability in the tropics is dominated by the El Niño–Southern Oscillation (ENSO) phenomenon (Philander 1990; Neelin et al. 1998). Models of various degrees of complexity are capturing different aspects of interannual variability in the tropics with increasing success (Latif et al. 1998). Linear and nonlinear investigations of ENSO dynamics show that distinct parameter values can modify the coupled model’s behavior, both quantitatively and qualitatively (Neelin et al. 1998). One way to improve a coupled model and its forecast skill is to estimate the correct values of its parameters from observations through data assimilation.

In this paper, we apply the EKF to help solve the parameter estimation problem (Ghil and Malanotte-Rizzoli 1991; Ghil 1997), by using the so-called augmented-state approach (e.g., Hao 1994; Hao and Ghil 1995), which incorporates one or more parameters as additional state variables. We will simultaneously estimate the state and two key parameters of an ICM, using first synthetic and then observational datasets. The advantage of using synthetic datasets is that the “true” values of the parameters are known, which helps us to evaluate the methodology. Observational datasets, on the other hand, provide insight on the application of the methodology to real problems.

Ongoing efforts in the Bayesian statistics community aim to develop a theoretical basis for the use of estimation techniques when models have poorly known errors, including those in parameters (Kennedy and O’Hagan 2001; Smith 2002; Goldstein and Rougier 2008). This paper’s Kalman filter approach falls broadly into the same category: from a Bayesian point of view, its objective is to estimate not only the variables of interest, whether the model’s state or its parameters, but also some measure of uncertainty in these estimates (Ihler et al. 2007).

The structure of the present paper is as follows. In section 2, we review briefly the coupled ocean–atmosphere model used in the study. In section 3, the EKF methodology for nonlinear models is briefly reviewed, and the augmented-state approach to parameter estimation using the EKF is presented; the appendix outlines how this approach can be generalized to other sequential estimators, in particular to the partitioned Kalman filter (Fukumori 2002). Our main results on state and parameter estimation studies for the coupled model, using both synthetic and observational data, are reported in section 4. Concluding remarks appear in section 5.

2. The model

We use the ICM of Jin and Neelin (1993), which is essentially a further idealization of the Zebiak and Cane (1987) model. The major simplification is to treat explicitly only the zonal dependence of sea surface temperature (SST) over an equatorial strip, while the meridional structure of the associated atmospheric forcing is specified.

The meridional structure of the mean currents in the tropical ocean’s upper layer is projected onto a basis of Hermite functions, truncated here at M = 14 meridional modes. The coupled model thus boils down to a spatially one-dimensional system of evolution equations in the zonal direction x and time t for the SST, the amplitude q0 of the oceanic Kelvin wave, and the amplitudes qn, n = 2, 4, . . . , 14, of seven Rossby waves. The total length of the ocean basin is 150°, from 130°E to 80°W, with a grid spacing of 6.25°. The grid points in the x direction are numbered from 1 to 25, going from the basin’s western to its eastern boundary; these points are also used to identify the location of observing sections. The coupled system so discretized has (1 + 8) × 24 = 216 degrees of freedom, and the time step is 6 h. Full details on this ICM appear in appendix A of Part I.

This coupled model presents two different kinds of ENSO oscillations. The nature and properties of the oscillation depend on the values of two key parameters, namely, the relative coupling coefficient μ and the surface layer coefficient δs.

An increase of the coupling coefficient μ enhances the nonlinear impact of local coupling processes between the atmosphere and ocean [see Jin and Neelin (1993), Hao et al. (1993), and Eq. (A20) in Part I]. When μ is very small, the system behaves as an uncoupled model, with SST decaying toward a steady climatological solution. As μ increases, the system undergoes a Hopf bifurcation, and starts to oscillate with many features that resemble the observed ENSO cycle (Jin et al. 1994, 1996).

For small values of the surface layer coefficient δs, the model behaves like the delayed oscillator of Suarez and Schopf (1988) or Battisti and Hirst (1989), while it tends to have westward-propagating features in the SST anomalies when δs approaches unity. Moreover, δs helps determine the period of the resulting oscillation (Hao et al. 1993; Neelin et al. 1998).

In this paper we consider mainly two combinations of parameter values; the pair (δs = 0, μ = 0.76) that corresponds to the model’s delayed-oscillator regime, while (δs = 0.8, μ = 0.56) yields westward-propagating behavior. The remaining model parameter values appear in Table 1 of Part I.

The evolution of the equatorial Pacific SST anomaly field for these two modes is quite different, as can be seen in Fig. 1, where results from a 30-yr model run are shown. In the delayed-oscillator regime (Fig. 1a) the SST anomalies are confined to the eastern part of the basin, with a slight eastward propagation. Since δs is set to zero here, the thermocline feedback dominates over surface layer processes (Hao et al. 1993). Large cyclical variations of the SST in the eastern part of the basin are mainly due to the shallowness of the mean thermocline there, and have a period of about 3.6 yr for this δs = 0. In the western part, the SST anomalies are smaller, since the mean thermocline depth there is much greater. In contrast, the model run for the westward-propagating regime (Fig. 1b), displays much smaller, westward-propagating SST anomalies that are much less well organized.

Fig. 1.

Time evolution of equatorial Pacific SST anomalies (SSTA; °C) for the (a) delayed-oscillator and (b) westward-propagating modes of the ICM; and (c) observations of monthly SST data from the Climate Data Library at the IRI/LDEO (see online at http://ingrid.ldeo.columbia.edu/).

Fig. 1.

Time evolution of equatorial Pacific SST anomalies (SSTA; °C) for the (a) delayed-oscillator and (b) westward-propagating modes of the ICM; and (c) observations of monthly SST data from the Climate Data Library at the IRI/LDEO (see online at http://ingrid.ldeo.columbia.edu/).

Figure 1c shows observational SST anomaly data (Kaplan et al. 1998) for the equatorial Pacific over 30 yr, from 1975 to 2004, with the seasonal cycle removed; these SST anomalies share features of both Figs. 1a,b. It is quite obvious, though, that our ICM is highly idealized and that not all features present in the observations can be adequately captured by using such a model. Several questions arise at this point, within a broader perspective of a hierarchy of models for the understanding, simulation and seasonal-to-interannual prediction of ENSO-related climate phenomena (Ghil and Robertson 2002). First, would a better set of parameters, which allows for the interaction of the two modes, standing and westward propagating, give a better match to the observations? Second, can this match be further improved by allowing for parameter values that vary in time, possibly reflecting changes in processes that are not explicitly captured by this ICM?

The answer to the first question is still often given by running the model many times over to find a “best match” with observations, by using various parameter combinations. Even if aiming merely at better, but still constant parameter values, the trial-and-error approach may not be a practical way to achieve such estimates, since the possible combinations cannot be exhausted when the number of state variables or the number of parameters is large. The results in Fig. 1 thus indicate the need for more accurate, automated techniques of estimating the model parameters by using an optimized combination of observed data and model results. The Kalman filter described in the next section is capable of providing such a combination. As we shall see in section 4, the EKF can even provide a convenient answer to the second question raised above.

3. The Kalman filter and parameter estimation

A practical way to include estimation of model parameters into the Kalman filter is by the so-called state augmentation method (Gelb 1974; Hao and Ghil 1995; Galmiche et al. 2003; Kao et al. 2006; Kondrashov et al. 2007), in which the parameters of interest are treated as additional state variables. For a linear system, the numerical algorithm for advancing the state vector x from time kΔt to time (k + 1)Δt is

 
formula

Here xk = x(kΔt) represents a state column vector, composed of all model variables, and the matrix 𝗠 is obtained by discretizing the partial differential operator. Superscripts “f” and “a” refer to forecast and analysis, respectively, with xak being the best estimate of the state vector at the time k, based on the model and the observations available so far. The evolution of xt, where superscript “t” refers to true, is then assumed to differ from the model by a random error ε. This model error is commonly taken to be a Gaussian white noise sequence, with mean zero and model-error covariance matrix 𝗤, Eεk = 0 and EεkεTl = 𝗤kδkl, where E is the expectation operator and δkl is the Kronecker delta.

For simplicity, let us assume that there is only one model parameter μ: 𝗠 = 𝗠(μ). We can define equations for evolving the parameter’s “forecast” and true values, by assuming, in the absence of additional information, the following persistence model:

 
formula

with Eεμk = 0 and Eεμ2k = qμ. When additional information is available, Eq. (2) can be generalized to allow for more complex spatial and temporal dependence; such dependence may include, for instance, a seasonal cycle (e.g., Kondrashov et al. 2005).

Next, we form an augmented-state vector x, model M and error :

 
formula

and rewrite our model equations for the augmented system:

 
formula

The situation of interest is one in which μ itself is not observed, so

 
formula

The observation matrix 𝗛k accounts for the fact that usually the dimension of yok is less than the dimension of xtk, that is, at any given time observations are not available for all numerical grid locations. In addition, 𝗛k represents transformations that may be needed if other variables than the state vector are observed, as well as any required interpolation from observation locations to nearby numerical grid points. Matrix 𝗛k is the observation matrix 𝗛k augmented by the one column of zeros to account for the fact that parameter μ is not observable.

The observational error εo is also taken to be Gaussian, white in time, with mean zero and given covariance matrix 𝗥, EεokεoTl = 𝗥kδkl. Moreover, one commonly assumes, unless additional information is available, that model error and observational error are mutually uncorrelated, EεokεoTl = 0.

The optimal gain matrix 𝗞k is computed by minimizing the analysis error variance trace(𝗣ak), that is, the expected mean-square error between analysis and the true state. This Kalman gain matrix represents the optimal weights given to the observations in updating the model state vector. The Kalman filter equations for the augmented system then become

 
formula

The analysis step for the augmented system involves only observations of the state:

 
formula

while the augmented error covariance matrices involve cross terms between the state variables and the parameter. Dropping from now on the time subscript k, we have

 
formula

Using the definition of H in Eq. (5), we obtain

 
formula

The augmented model propagates the forecast error of the parameter into the cross-covariance term 𝗣 fμx. By substituting Eq. (9) into Eq. (7), we can readily see that this error propagation enables the EKF to extract information about the parameter from the state observations and to update the unobserved parameter at the analysis step:

 
formula

This formulation can be easily extended to the case when several unknown parameters have to be estimated and μ then becomes a vector instead of a scalar (Ghil 1997).

Equations (9) and (10) show that, if the state estimator is given, the innovation of parameter estimates requires only knowledge of the cross-covariance matrix 𝗣 fμx, which is readily available from Eq. (6) for the propagation of the augmented error covariance matrix. This feature allows one to combine the augmented-state approach with other EKF-type algorithms, including the EnKF (Evenson 2003) and reduced-rank Kalman filters (Fukumori and Malanotte-Rizzoli 1995; Cohn and Todling 1996; Fukumori et al. 1999; Tippett et al. 2000), which rely on various approximations for computing covariance matrices. In the appendix, we show how parameter estimation can be implemented, at a very low additional cost, for a particular reduced-rank EKF, the partitioned Kalman filter (Fukumori 2002).

The Kalman gain is optimal for a linear system, when both 𝗠(x) = 𝗠x and 𝗛(x) = 𝗛x, and under the Gaussian noise assumptions mentioned above. In this case, the gain is based on the correct estimation of forecast error covariances from initial uncertainties, model errors, and model dynamics.

Even when the original model 𝗠 is linear, the augmented model becomes nonlinear as soon as one or more parameters are included into the augmented model state, and the model contains products or other nonlinear functions of the state variables and parameters (Jazwinski 1970; Gelb 1974; Kondrashov et al. 2007). In any such case, one needs to use the EKF formulation or some other sequential-estimation approach that can deal with a nonlinear model. In our application, the SST equation is nonlinear to start with, and so is the augmented model.

In the EKF, the nonlinear model is linearized around the current state when estimating the propagation of the forecast error, while the full nonlinear model is still used to advance the state. So we will use the linearizations and of the augmented (x) and (x), respectively, about the current augmented state x = xfk to propagate the error covariances and compute the Kalman gain matrix in Eq. (6):

 
formula

where indices i and j refer to a particular matrix and state vector entry. The EKF is first-order accurate in general but it may diverge in the presence of strong nonlinearities (Miller et al. 1994; Chin et al. 2007). The identical-twin experiments in section 4a demonstrate that, in our case, we can obtain reliable and robust estimates of both the state and parameters with the EKF.

We apply sequential estimation to the relative coupling coefficient μ and surface layer coefficient δs, along with the state, and so we need to linearize around the current values of the augmented state vector formed by the state vector itself and the two parameters, μ and δs. Since we utilize an explicit scheme, linearization with respect to the state vector and the two estimated parameters is readily available and it follows from the known coefficients of 𝗠. For an implicit scheme, linearization with respect to the estimated parameters can be more laborious. In such cases, one can use small perturbations in the parameter values on the right-hand side of the governing equations and then apply numerical differentiation (Kondrashov et al. 2007).

a. Error covariances

1) Model error covariance 𝗤

The most severe errors in tropical ocean prediction seem to arise from wind stress errors (Leetmaa and Ji 1989; Graham et al. 1992; Hao and Ghil 1994; Miller et al. 1995). We introduce, therefore, model errors at every time step by adding noise to the wind stress that is obtained from the model atmosphere’s response to SST anomalies. The wind stress error is specified as in Miller and Cane (1989), Miller et al. (1995) and Cane et al. (1996). It is taken as white in time and Gaussian correlated and homogeneous in space. The errors are assumed to have the same meridional structure as the atmospheric response to the SST anomalies.

The particular form used for each component of the wind stress error covariance matrix 𝗤w is given by

 
formula

where ek,i and el,j are the wind stress errors at locations i and j, at time steps k and l, respectively; δkl is the Kronecker delta; στ is the magnitude factor of the wind tress error; and Lx is the prescribed decorrelation distance. We follow Miller and Cane (1989) and use στ = 0.02 Pa and Lx = 10° of longitude. Model error is constructed by projecting the wind stress error onto the prognostic model variables, as described in appendix B of Part I.

The weights obtained from Eq. (10) for the innovations in the parameter μ are proportional to the size of its cross-covariance forecast error 𝗣 fμx; the latter, in turn, is related to the parameter’s contribution qμ to the model error 𝗤; see Eqs. (2), (6), and (8). This value qμ should be chosen according to how much variation we allow the estimated parameter to have and also to how much information is needed from the observations of the state. Since a smooth estimation of the parameters is often required, small error values tend to be a good choice: here we used 2% of the parameter’s squared initial values for qμ.

2) Observation error covariance 𝗥

We prescribe the observation error covariance matrix 𝗥 as a diagonal matrix, 𝗥ij = σ2iδij, where σi is the standard deviation of the noise in the observations of the state variable xi. The observations are represented as measurements T of SST, amplitudes qn of the oceanic waves, or atmospheric wind stress data τ at longitudinal locations; they are all taken twice a month if no other specification is given.

The observations are contaminated by white noise, with the following standard deviations: σoT = 0.5 K, σoq0 = 0.02, σoqn = 0.01, and σoτ = 0.01 Pa. Note that the standard deviations for all wave coefficients are in nondimensional units. The equivalent observational errors for the corresponding physical variables, namely, zonal current u, thermocline depth anomaly h, and vertical velocity w, are σou = 0.02 m s−1, σoh = 1.2 m, and σow = 0.275 cm day−1.

3) Initial forecast error covariance 𝗣 f0

The initial errors are white in space, with the following standard deviations at all grid points: σT = 0.9 K, σq0 = 0.06, and σqn = 0.04, for n = 2, 4, 6; and σqn = 0.03, for n = 8, 10, 12, 14. The standard deviations for all wave coefficients are in nondimensional units, as for the observation errors. The initial forecast error covariance 𝗣 f0 is assumed to have a diagonal structure, with its elements equal to the variances of the state variables at t = 0: 𝗣 f0ij = σ2iδij.

4. Results and discussion

a. Identical-twin experiments

To test the parameter estimation scheme described in section 3b, we first conduct identical-twin experiments in which both the true state and model parameters are known. We call the model with correct parameter values the “perfect model,” and the correct parameter values (μ = 0.76, δs = 0) are referred to as the “reference values.” The “observations” are drawn from the perfect model solutions with observation errors superimposed.

The forecast run differs from the “nature” or “control” run not only because of initial state errors and model errors, but also because of incorrect model parameter values. Using standard data assimilation terminology (Bengtsson et al. 1981; Ide et al. 1997), the forecast errors are defined as the instantaneous root-mean-square (RMS) difference between the model forecast xf and the nature run xt at any given time tk [see Eq. (4)]. The time-dependent forecast error variances are estimated by the Kalman filter algorithm as diagonal elements of the covariance matrix 𝗣f, and represent an ensemble average for different realizations of . For the identical-twin experiments, these errors can be computed directly using the known “truth,” which helps to verify the given data assimilation scheme’s optimality; however, there is no way of knowing the true evolution of the system when using real data, given the partial and inaccurate nature of atmospheric and oceanic datasets (see Ghil et al. 1981; Ghil and Malanotte-Rizzoli 1991).

Two types of assimilation runs with the EKF are performed: one updates the model state only; the other updates both model state and parameters. We use in either case SST and subsurface currents at a single location, situated at 111°W, as well as a more plentiful dataset from 15 equatorial Pacific locations that resembles the Tropical Atmosphere–Ocean (TAO) array; these two datasets are referred to as “1 data” and “TAO data,” respectively. We start the forecast model with incorrect parameter values and prescribed model error (i.e., the wind stress errors added at each time step).

1) Relative coupling coefficient μ only

Here, the relative coupling coefficient μ is initially set to 0.56, while its reference value is 0.76; the value of δs is kept fixed, δs = 0. Without including parameter estimation, μ stays the same at the initial value. Since this smaller value of μ places the model in the climatological steady regime, the forecast SST anomalies simply decay to zero.

Five-year-averaged RMS errors for the SST and oceanic waves are given in Table 1. The forecast errors in the first row of the table, (FCST), are for a model run with incorrect μ and without data assimilation.

Table 1.

Five-year-averaged RMS errors for FCST, the forecast run with the incorrect coupling coefficient μ = 0.56; ASML (1 data), state estimation only using the 1 data and the incorrect μ value; ASML (TAO data), same as previous row, but using the full TAO data; ASML(1 data, μ), assimilation run with the μ corrected using the single-section data; and ASML(TAO data, μ), same as previous row, but using the full TAO data.

Five-year-averaged RMS errors for FCST, the forecast run with the incorrect coupling coefficient μ = 0.56; ASML (1 data), state estimation only using the 1 data and the incorrect μ value; ASML (TAO data), same as previous row, but using the full TAO data; ASML(1 data, μ), assimilation run with the μ corrected using the single-section data; and ASML(TAO data, μ), same as previous row, but using the full TAO data.
Five-year-averaged RMS errors for FCST, the forecast run with the incorrect coupling coefficient μ = 0.56; ASML (1 data), state estimation only using the 1 data and the incorrect μ value; ASML (TAO data), same as previous row, but using the full TAO data; ASML(1 data, μ), assimilation run with the μ corrected using the single-section data; and ASML(TAO data, μ), same as previous row, but using the full TAO data.

With the incorrect μ, we performed state estimation using the TAO-data and 1-data datasets. The results, shown as ASML (TAO data) and ASML (1 data) in Table 1, indicate that the EKF estimation of the model state alone does reduce the forecast errors significantly. In the case of ASML (TAO data), these errors, in both SST and waves, are just slightly larger than those using the same data with the perfect model. The improvements suggest that assimilation of observations can compensate for the deficiency in the model parameter μ, even without correcting it. This is good news for constructing consistent fields from observations when using a model that might have some errors in its parameters, provided one applies an advanced data assimilation method like the EKF.

For the purposes of seasonal-to-interannual forecasting, however, the error in μ will prevent the model from showing the correct oscillatory features even when the initial data have been corrected through data assimilation; therefore, parameter estimation is still of the essence. In updating μ, the dataset 1 data can be very useful, as shown under ASML (1 data, μ) in Table 1. Even this minimal amount of data, when using the EKF to estimate μ, helped restore the oscillation with the correct phase and period, as well as the amplitude of the SST anomalies (not shown). The improvements in all model fields here are also comparable to those using the same dataset with the perfect model.

The results above suggest that the observations of both SST and ocean waves at a single location, preferably in the eastern basin (see Hao and Ghil 1994 and Part I), are sufficient for both the state and parameter estimation, when only the μ value is poorly known in this ICM.

The ASML (TAO data, μ) results in the table, though, show larger improvements of both model state and parameter, when a larger dataset, the TAO data, is available. In fact, the improvements in model states are identical to those with the perfect model using the same data (not shown).

2) Surface layer coefficient δs only

In this experiment, the surface layer coefficient is initially set to a large value of δs = 0.8, while its reference value is zero; the value of μ is kept at its reference value: μ = 0.76. Five-year-averaged forecast errors are given in Table 2 (FCST); these are larger than those with the incorrect μ value in Table 1. Updating only model states with the EKF does improve substantially the state fields; see ASML (1 data) and ASML (TAO data) in Table 2.

Table 2.

Five-year-averaged RMS errors as in Table 1, but for the surface layer coefficient δs, which is kept at the incorrect value of 0.8 in the first three rows, and estimated along with the state in the last two rows. The results of state estimation only using the 1 data and TAO data appear in the second and third row, respectively.

Five-year-averaged RMS errors as in Table 1, but for the surface layer coefficient δs, which is kept at the incorrect value of 0.8 in the first three rows, and estimated along with the state in the last two rows. The results of state estimation only using the 1 data and TAO data appear in the second and third row, respectively.
Five-year-averaged RMS errors as in Table 1, but for the surface layer coefficient δs, which is kept at the incorrect value of 0.8 in the first three rows, and estimated along with the state in the last two rows. The results of state estimation only using the 1 data and TAO data appear in the second and third row, respectively.

When using the TAO data, the assimilation errors for the oceanic waves are similar to those in the corresponding cases when the perfect model was employed, but SST errors are much larger (cf. Table 2 in Part I and also Table 1 here). This difference in the degree of improvement for the SST and oceanic wave fields is due to the fact that δs values affect mainly the SST equation, but do not impact the oceanic waves directly. These results suggest that, with oceanic wave observations alone, it would be difficult to obtain an accurate SST analysis with the value of δs left uncorrected.

Estimating δs along with the state by using the 1-data dataset improves the SST field relative to ASML (TAO data), while the results for the ocean waves are only slightly poorer, see ASML (1 data, δs) in Table 2. The 1-data dataset, however, is barely sufficient for the estimation of δs itself. On the other hand, the last row ASML (TAO data, δs) in Table 2 shows that δs and the model state can be greatly improved when a larger number of observations, as in TAO data, are used.

3) Both parameters μ and δs

Figure 2 shows results from an ICM experiment with synthetic SST observations, but not subsurface currents, from the “TAO” array to estimate both the μ and δs parameters. The assimilation run starts with incorrect model parameters that correspond to the westward-propagating regime (μ = 0.56, δs = 0.8); these parameter values are corrected by assimilating observations taken from the control run, that uses delayed-oscillator parameter values (μ = 0.76, δs = 0).

Fig. 2.

Time evolution of estimated parameters (a) μ and (b) δs for the identical-twin assimilation run with the observations for SST data taken from 15 equatorial locations in the east-central Pacific. The correct (true) values of μ and δs are equal to 0.76 and 0, respectively.

Fig. 2.

Time evolution of estimated parameters (a) μ and (b) δs for the identical-twin assimilation run with the observations for SST data taken from 15 equatorial locations in the east-central Pacific. The correct (true) values of μ and δs are equal to 0.76 and 0, respectively.

The parameters adjust within 1 yr (heavy blue line) to fluctuate about their true values (dashed black line) and remain within the EKF’s estimated standard deviation (solid red line in Fig. 2), thus recovering the delayed-oscillator regime in the evolution of the state variables (not shown). The convergence time depends largely on the assumed model error of the parameters: increasing this error will usually accelerate the convergence. Numerical sensitivity experiments (not shown) confirm that other combinations of true and “incorrect” parameter values did not produce any adverse effects on the convergence of the parameter estimation process.

In the EKF framework, corrections of the full state vector, consisting of both the SST and ocean wave fields, are made in the analysis step in accordance with the new SST measurements alone. Figure 3 shows RMS errors in the zonal surface current from EKF experiments with and without parameter estimation. When no data assimilation is performed, the errors are very large (black curve). Even though the wave fields are not directly observed, our assimilation scheme still provides considerable improvement in their estimates, with the best results obtained when we estimate both the state and parameters simultaneously (blue curve).

Fig. 3.

RMS errors in the zonal surface current for the identical-twin experiment of Fig. 2. Blue and red lines are for actual errors with and without parameter estimation, respectively; and the black line indicates the pure forecast error.

Fig. 3.

RMS errors in the zonal surface current for the identical-twin experiment of Fig. 2. Blue and red lines are for actual errors with and without parameter estimation, respectively; and the black line indicates the pure forecast error.

The surface current errors in Fig. 3 contain a very large periodic component; such a component is present also, albeit much weaker, in the errors made in estimating the two parameters in Fig. 2. This periodic component in the errors is due to our identical-twin setup, where the truth is dominated by the delayed-oscillator mode, while the incorrect model contains no such highly periodic signal (cf. Figs. 1a,b).

The rapid convergence of the parameter estimates—in spite of the obviously chaotic, unstable character of the ENSO model dynamics—suggests extending the theory for state estimation in nonlinear forecast-assimilation systems (Carrassi et al. 2008) to the combined, state and parameter estimation problem. Such an extension could further optimize the convergence of parameter estimation in the augmented-state approach.

b. Experiments using actual SST observations

To explore the applicability of this methodology using real-world observations, we assimilated monthly SST observations over the 30-yr interval from January 1975 to December 2004, obtained from the Climate Data Library at the International Research Institute for Climate Prediction/Lamont-Doherty Earth Observatory (IRI/LDEO). This dataset is a statistically homogenous concatenation of Kaplan et al.’s (1998) SST fields and the National Centers for Environmental Prediction (NCEP) optimal-interpolation SST fields of Reynolds and Smith (1994); the latter dataset is on a finer grid (i.e., 1° × 1°) and was regridded to the resolution of the former (i.e., 5° × 5°).

We first interpolated the observations to the equatorial grid points of our ICM and removed the seasonal cycle. To match the climatology of the model and observations, we removed the climatology of the observations at each grid point and then added the model climatology to the observational anomalies. The climatological mean of the modified observational dataset thus obtained is the same as the basic state of the model.

To achieve robust parameter estimates over the whole assimilation interval, we perform the parameter estimation iteratively, using the previous parameter estimates in the next forward estimation. Recall from Eq. (10) that parameter changes are related to their cross covariances with the model state; these cross covariances are not generally known in advance, so they are set to zero initially and modified on the next pass through the observations. Also, unlike in the identical-twin case, the true model parameters are neither known nor necessarily constant in time, and so their initial values have to be educated guesses at best. The first pass of the EKF helps to establish the cross covariances in the assimilation cycle, while providing preliminary parameter estimates. As we shall see, the second pass achieves a robust, albeit time-varying parameter estimate over most of the record.

Convergence of parameter estimates within two such iterations, for both δs and μ, is illustrated in Fig. 4. The estimates of μ (Fig. 4b) do not change much between the first and second iteration (blue and red curves, respectively, in Figs. 4b,c), while δs (Fig. 4c) varies significantly over an initial, 8-yr time interval (1975–83). When using quite different initial parameter values δs0, the second iterates of the δs estimation agree very well (Fig. 4d); similar results (not shown) are obtained when estimating μ. The remaining discrepancy in the initial transition interval (1977–83), in Figs. 4c,d, could be attributed to the very small SST anomalies in that interval (see Figs. 1a and 4a); these small anomalies do not constrain the surface layer parameter δs in an effective manner.

Fig. 4.

Time evolution of estimated parameters, during an assimilation run: (a) observations of monthly SST data from IRI/LDEO Climate Data Library, taken from 15 equatorial locations in the east-central Pacific; (b) evolution of the nonlinear coupling parameter μ; and (c), (d) evolution of the surface layer parameter δs. In (c) two successive iterations of our δs estimate are illustrated for an initial parameter value of δs0 = 0.3, while (d) shows the second iterate of the δs estimation for three distinct initial values, δs0 = 0.3, 0.5, and 0.8.

Fig. 4.

Time evolution of estimated parameters, during an assimilation run: (a) observations of monthly SST data from IRI/LDEO Climate Data Library, taken from 15 equatorial locations in the east-central Pacific; (b) evolution of the nonlinear coupling parameter μ; and (c), (d) evolution of the surface layer parameter δs. In (c) two successive iterations of our δs estimate are illustrated for an initial parameter value of δs0 = 0.3, while (d) shows the second iterate of the δs estimation for three distinct initial values, δs0 = 0.3, 0.5, and 0.8.

The parameters can vary on fairly short time scales and switch between values that are associated with one or the other of the two distinct modes of ENSO behavior, the delayed-oscillator and westward-propagating mode. Rapid adjustments of these parameters occur, in particular, during the two strong ENSO events, 1982/83 and 1997/98, during which the uncertainty in both parameters becomes very small (thin black curves in Figs. 4b,c).

Figure 5 shows the evolution of the equatorial Pacific SST anomaly field obtained from state estimation only (Fig. 5b) and simultaneous state and parameter estimations (Fig. 5a). The experiment with only state estimation held the parameters δs and μ constant at values in the delayed-oscillator regime. The main differences between the two panels are the westward spatial extent, intensity, and duration of the strong ENSO events of 1982/83 and 1997/98, which are better captured by the EKF with parameter estimation.

Fig. 5.

Time–longitude plots (Hovmöller diagrams) of SST anomaly results (°C) for three EKF experiments: (a) combined state and parameter estimation, (b) state estimation only with constant parameters from the delayed-oscillator mode, and (c) model simulation with a time history of estimated μ(t) and δs(t), but with no further data assimilation.

Fig. 5.

Time–longitude plots (Hovmöller diagrams) of SST anomaly results (°C) for three EKF experiments: (a) combined state and parameter estimation, (b) state estimation only with constant parameters from the delayed-oscillator mode, and (c) model simulation with a time history of estimated μ(t) and δs(t), but with no further data assimilation.

Figure 5c shows the SST anomalies obtained when forcing the model with the time history of estimated δs(t) and μ(t), but without assimilating any other data. It is interesting to note that this model simulation eliminated the regular ENSO oscillations in the delayed-oscillator regime. The model response is rather modest, except for the two strong ENSO events during 1982/83 and 1997/98, when the “correct” history of the parameters seems to reproduce the SST anomalies obtained from EKF with state estimation alone.

Assimilating SST thus improves the model’s SST estimates, as expected. This result does not guarantee, however, that unobserved model fields, such as currents, would also be improved. Figure 6 shows the anomaly field for the equatorial zonal surface current that is obtained from simultaneous state and parameter estimation (Fig. 6a), as well as the result when performing state estimation only, with the fixed parameters in the delayed-oscillator regime (Fig. 6b). The main difference between these two assimilation experiments is the intensification of the currents during the strong ENSO event of 1997/98, which is better captured by the EKF with state and parameter estimation.

Fig. 6.

Hovmöller diagrams of zonal surface current anomalies (m s−1). Two EKF experiments: (a) combined state and parameter estimation and (b) state estimation for constant parameters from the delayed-oscillator mode. (c) Observed surface current anomalies from OSCAR data for the 12-yr interval 1993–2004.

Fig. 6.

Hovmöller diagrams of zonal surface current anomalies (m s−1). Two EKF experiments: (a) combined state and parameter estimation and (b) state estimation for constant parameters from the delayed-oscillator mode. (c) Observed surface current anomalies from OSCAR data for the 12-yr interval 1993–2004.

This intensification is in qualitative agreement with the observational data displayed in Fig. 6c, which plots 1993–2004 surface current anomalies from the Ocean Surface Current Analyses–Real Time (OSCAR) dataset, derived from satellite altimeter and scatterometer data (Bonjean and Lagerloef 2002). The OSCAR data along the equator (0.5°N–0.5°S) were interpolated to the ICM’s grid and then smoothed by using a 3-month running mean and spatial averaging over the 3 nearest grid points.

The results in Figs. 4 –6 suggest that our ICM is too idealized to represent in detail the complex evolution patterns of the observed SST and surface currents, but it has skill in the eastern tropical Pacific, especially when the ENSO signals are strong. Moreover, missing mechanisms in the physical model may be at least partially compensated for by skillful parameter estimation in seasonal-to-interannual prediction.

These results encourage us to extend parameter estimation to large, state-of-the-art GCMs, where the improvement might be more substantial, since these models are much more realistic and hence closer to the observations. The main obstacle in doing so is the size and nonlinearity of the problem in this case (Ghil 1997). We propose in the appendix an approach to solving the size problem.

5. Conclusions

In this paper, we demonstrated parameter estimation in an intermediate coupled model (ICM) using the EKF and the “augmented state” approach. The EKF method was first used to estimate the state and parameters in an identical-twin framework, based on incomplete and inaccurate observations of the model state. Model-generated oceanic datasets (Fig. 1) were assimilated to produce a time continuous, dynamically consistent description of the model’s El Niño–Southern Oscillation (ENSO). Two unknown model parameters, the coupling coefficient μ and surface layer coefficient δs, were estimated by treating them as additional state variables. The results of section 4a show that the augmented-state approach in the EKF methodology can estimate these two parameters successfully (Figs. 2 and 3). In addition, state estimation with an imperfect model can be improved by including the estimation of model parameters in the assimilation process.

Errors in the model parameters can impact the assimilation results in different ways. An error in the relative coupling coefficient μ did little harm in estimating the state (Table 1). A large error in the surface layer coefficient δs, however, made it difficult to assimilate SST anomalies well without correcting the parameter (Table 2). Assimilation of observations can thus be improved if we identify first the erroneous model parameter(s). On the other hand, all parameter errors will hamper successful model forecasts. Identifying a priori an erroneous model parameter and assigning a proper error to it are left for further work.

In section 4b, we performed additional experiments to evaluate the applicability of this method to real-world SST observations. In this case, the true parameters are not known and are not necessarily constant over time. The parameters adjust on fairly short time scales, in particular during strong ENSO events (Fig. 4); they switch between values that are associated with one or the other of the two distinct modes of ENSO behavior, the delayed-oscillator and westward-propagating mode. The results show that the estimation of the unknown parameters still helps improve the state estimation (Fig. 5), and that missing mechanisms in the physical model may be partially compensated for by skillful parameter estimation.

The impact of parameter estimation on SST evolution when using SST observational data is less notable. An encouraging result is that our assimilation scheme does provide considerable improvement in current estimation (Fig. 6) even though the flow fields are not directly observed. Performing both state and parameter estimation yields again the best results.

Future work will include applying EKF parameter estimation efficiently to much larger, state-of-the-art, coupled ocean–atmosphere general circulation models (GCMs). The key idea in our approach to this problem (see the appendix) is that the number of parameters one wishes to estimate, even for a global coupled GCM, is small relative to the number of state variables. Thus, given any sequential-type filter (Ghil and Malanotte-Rizzoli 1991) for state estimation, the increase in rank for the augmented-state approach to parameter estimation is quite small, which offers a good opportunity for efficiently estimating parameters in a GCM. To do so, we have shown that the evolution of the parameter estimates in the EKF estimation process can be separated from that of the state estimates. It is this separation that is crucial and extending it to any efficient EKF implementation, such as the partitioned Kalman filter of Fukumori (2002), can help to estimate the parameters of the coupled GCM, along with its state.

Acknowledgments

It is a pleasure to thank I. Fukumori and J. D. Neelin for stimulating discussions and concrete suggestions, and two anonymous reviewers for constructive comments. The sea surface temperature (SST) data for 1975–2004 were obtained online at http://ingrid.ldeo.columbia.edu/, while the ocean current data for 1993–2004 were obtained online at http://www.oscar.noaa.gov/. This work was supported by NASA’s Modeling, Analysis and Prediction (MAP) Program through Grant 1281080.

REFERENCES

REFERENCES
Annan
,
J. D.
,
J. C.
Hargreaves
,
N. R.
Edwards
, and
R.
Marsh
,
2005
:
Parameter estimation in an intermediate complexity earth system model using an ensemble Kalman filter.
Ocean Modell.
,
8
,
135
154
.
Battisti
,
D.
, and
A.
Hirst
,
1989
:
Interannual variability in a tropical atmopshere–ocean model: Influence of the basic state, ocean geometry, and nonlinearity.
J. Atmos. Sci.
,
46
,
1687
1712
.
Bengtsson
,
L.
,
M.
Ghil
, and
E.
Källén
,
1981
:
Dynamic Meteorology: Data Assimilation Method.
Springer-Verlag, 330 pp
.
Bonjean
,
F.
, and
G. S. E.
Lagerloef
,
2002
:
Diagnostic model and analysis of the surface currents in the tropical Pacific Ocean.
J. Phys. Oceanogr.
,
32
,
2938
2954
.
Cane
,
M. A.
,
A.
Kaplan
,
R. N.
Miller
,
B.
Tang
,
E. C.
Hackert
, and
A. J.
Busalacchi
,
1996
:
Mapping Tropical Pacific sea level: Data assimilation via a reduced state space Kalman filter.
J. Geophys. Res.
,
101
,
22599
22617
.
Carrassi
,
A.
,
M.
Ghil
,
A.
Trevisan
, and
F.
Uboldi
,
2008
:
Data assimilation as a nonlinear dynamical systems problem: Stability and convergence of the prediction–assimilation system.
Chaos
,
18
.
023112, doi:10.1063/1.2909862
.
Chin
,
T. M.
,
M. J.
Turmon
,
J. B.
Jewell
, and
M.
Ghil
,
2007
:
An ensemble-based smoother with retrospectively updated weights for highly nonlinear systems.
Mon. Wea. Rev.
,
135
,
186
202
.
Cohn
,
S.
, and
R.
Todling
,
1996
:
Approximate data assimilation schemes for stable and unstable dynamics.
J. Meteor. Soc. Japan
,
74
,
63
75
.
Evenson
,
G.
,
2003
:
The ensemble Kalman filters: Theoretical formulation and practical implementation.
Ocean Dyn.
,
53
,
343
367
.
Fukumori
,
I.
,
2002
:
A partitioned Kalman filter and smoother.
Mon. Wea. Rev.
,
130
,
1370
1383
.
Fukumori
,
I.
, and
P.
Malanotte-Rizzoli
,
1995
:
An approximate Kalman filter for ocean data assimilation; an example with an idealized Gulf Stream model.
J. Geophys. Res.
,
100
,
6777
6793
.
Fukumori
,
I.
,
R.
Raghunath
,
L.
Fu
, and
Y.
Chao
,
1999
:
Assimilation of TOPEX/POSEIDON data into a global ocean circulation model: How good are the results?
J. Geophys. Res.
,
104
,
25647
25665
.
Galmiche
,
M.
,
J.
Sommeria
,
E.
Thivolle-Cazat
, and
J.
Verron
,
2003
:
Using data assimilation in numerical simulation of experimental geophysical flows.
C. R. Acad. Sci. (Mecanique)
,
331
,
843
848
.
Gelb
,
A.
,
1974
:
Applied Optimal Estimation.
The MIT Press, 374 pp
.
Ghil
,
M.
,
1997
:
Advances in sequential estimation for atmospheric and oceanic flows.
J. Meteor. Soc. Japan
,
75
,
289
304
.
Ghil
,
M.
, and
P.
Malanotte-Rizzoli
,
1991
:
Data assimilation in meteorology and oceanography. Advances in Geophysics, Vol. 33, Academic Press, 141–266
.
Ghil
,
M.
, and
A. W.
Robertson
,
2002
:
“Waves” vs. “particles” in the atmosphere’s phase space: A pathway to long-range forecasting?
Proc. Natl. Acad. Sci. USA
,
99
,
2493
2500
.
Ghil
,
M.
,
S.
Cohn
,
J.
Tavantzis
,
K.
Bube
, and
E.
Isaacson
,
1981
:
Applications of estimation theory to numerical weather prediction.
Dynamic Meteorology: Data Assimilation Methods, L. Bengtsson, M. Ghil, and E. Källén, Eds., Springer-Verlag, 139–224
.
Goldstein
,
M.
, and
J. C.
Rougier
,
2008
:
Reified Bayesian modelling and inference for physical systems.
J. Stat. Plann. Infer.
,
doi:10.1016/j.jspi.2008.07.019 (in press)
.
Graham
,
N. E.
,
T. P.
Barnett
, and
M.
Latif
,
1992
:
Considerations of the predictability of ENSO with a lower-order coupled model.
TAO Notes
,
7
,
11
15
.
Hao
,
Z.
,
1994
:
Data assimilation for interannual climate-change prediction. Ph.D. dissertation, Dept. of Atmospheric Sciences, University of California, Los Angeles, 223 pp
.
Hao
,
Z.
, and
M.
Ghil
,
1994
:
Data assimilation in a simple tropical ocean model with wind stress errors.
J. Phys. Oceanogr.
,
24
,
2111
2128
.
Hao
,
Z.
, and
M.
Ghil
,
1995
:
Sequential parameter estimation for a coupled ocean–atmosphere model. Proc. WMO Second Int. Symp. on Assimilation of Observations in Meteorology and Oceanography, WMO/TD-651, Tokyo, Japan, World Meteorological Organization, 181–186
.
Hao
,
Z.
,
J. D.
Neelin
, and
F. F.
Jin
,
1993
:
Nonlinear tropical air–sea interaction in the fastwave limit.
J. Climate
,
6
,
1523
1544
.
Ide
,
K.
,
P.
Courtier
,
M.
Ghil
, and
A. C.
Lorenc
,
1997
:
Unified notation for data assimilation: Operational, sequential and variational.
J. Meteor. Soc. Japan
,
75
,
181
189
.
Ihler
,
A.
,
S.
Kirshner
,
M.
Ghil
,
A. W.
Robertson
, and
P.
Smyth
,
2007
:
Graphical models for statistical inference and data assimilation.
Physica D
,
230
,
72
87
.
Jazwinski
,
A. H.
,
1970
:
Stochastic Processes and Filtering Theory.
Academic Press, 376 pp
.
Jin
,
F. F.
, and
J. D.
Neelin
,
1993
:
Modes of interannual tropical ocean–atmosphere interaction—A unified view. Part I: Numerical results.
J. Atmos. Sci.
,
50
,
3478
3503
.
Jin
,
F. F.
,
J. D.
Neelin
, and
M.
Ghil
,
1994
:
El Niño on the Devil’s staircase: Annual subharmonic steps to chaos.
Science
,
264
,
70
72
.
Jin
,
F. F.
,
J. D.
Neelin
, and
M.
Ghil
,
1996
:
El Nino–Southern Oscillation and the annual cycle: Subharmonic frequency locking and aperiodicity.
Physica D
,
98
,
442
465
.
Kao
,
J.
,
D.
Flicker
,
K.
Ide
, and
M.
Ghil
,
2006
:
Estimating model parameters for an impact-produced shockwave simulation: Optimal use of partial data with the extended Kalman filter.
J. Comput. Phys.
,
214
,
725
737
.
doi:10.1016/j.jcp.2005.10.022
.
Kaplan
,
A.
,
M.
Cane
,
Y.
Kushnir
,
A.
Clement
,
M.
Blumenthal
, and
B.
Rajagopalan
,
1998
:
Analyses of global sea surface temperature 1856–1991.
J. Geophys. Res.
,
103
,
18567
18589
.
Kennedy
,
M. C.
, and
A.
O’Hagan
,
2001
:
Bayesian calibration of computer models.
J. Roy. Stat. Soc.
,
B63
,
425
464
.
Keppenne
,
C.
,
2000
:
Data assimilation into a primitive-equation model with a parallel ensemble Kalman filter.
Mon. Wea. Rev.
,
128
,
1971
1981
.
Keppenne
,
C.
, and
M.
Rienecker
,
2002
:
Initial testing of a massively parallel ensemble Kalman filter with the Poseidon isopycnal ocean general circulation model.
Mon. Wea. Rev.
,
130
,
2951
2965
.
Kondrashov
,
D.
,
S.
Kravtsov
,
A. W.
Robertson
, and
M.
Ghil
,
2005
:
A hierarchy of data-based ENSO models.
J. Climate
,
18
,
4425
4444
.
Kondrashov
,
D.
,
Y.
Shprits
,
M.
Ghil
, and
R.
Thorne
,
2007
:
Kalman filter technique to estimate relativistic electron lifetimes in the outer radiation belt.
J. Geophys. Res.
,
112
.
A10227, doi:10.1029/2007JA012583
.
Latif
,
M.
, and
Coauthors
,
1998
:
A review of the predictability and prediction of ENSO.
J. Geophys. Res.
,
103
,
14375
14393
.
Leetmaa
,
A.
, and
M.
Ji
,
1989
:
Operational hindcasting of the tropical Pacific.
Dyn. Atmos. Oceans
,
13
,
465
490
.
Ljung
,
L.
,
1987
:
System Identification: The Theory for the User.
Prentice Hall, 519 pp
.
Miller
,
R. N.
, and
M. A.
Cane
,
1989
:
A Kalman filter analysis of sea level height in the tropical Pacific.
J. Phys. Oceanogr.
,
19
,
773
790
.
Miller
,
R. N.
,
M.
Ghil
, and
F.
Gauthiez
,
1994
:
Advanced data assimilation in strongly nonlinear dynamical systems.
J. Atmos. Sci.
,
51
,
1037
1056
.
Miller
,
R. N.
,
A. J.
Busalacchi
, and
E. C.
Hackert
,
1995
:
Sea surface topography fields of the tropical Pacific from data assimilation.
J. Geophys. Res.
,
100
,
13389
13425
.
Murphy
,
J. M.
,
D. M. H.
Sexton
,
D. N.
Barnett
,
G. S.
Jones
,
M. J.
Webb
,
M.
Collins
, and
D. A.
Stainforth
,
2004
:
Quantification of modelling uncertainties in a large ensemble of climate change simulations.
Nature
,
430
,
768
772
.
Navon
,
I. M.
,
1997
:
Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography.
Dyn. Atmos. Oceans
,
27
,
55
79
.
Neelin
,
J. D.
, and
Coauthors
,
1992
:
Tropical air–sea interaction in general circulation models.
Climate Dyn.
,
7
,
73
104
.
Neelin
,
J. D.
,
D. S.
Battisti
,
A. C.
Hirst
,
F. F.
Jin
,
Y.
Wakata
,
T.
Yamagata
, and
S. E.
Zebiak
,
1998
:
ENSO theory.
J. Geophys. Res.
,
103
,
14261
14290
.
Philander
,
S. G.
,
1990
:
El Niño, La Niña, and the Southern Oscillation.
Academic Press, 293 pp
.
Reynolds
,
R. W.
, and
T. M.
Smith
,
1994
:
Improved global sea surface temperature analyses.
J. Climate
,
7
,
929
948
.
Smedstad
,
O. M.
, and
J. J.
O’Brien
,
1991
:
Variational data assimilation and parameter estimation in an equatorial Pacific Ocean model.
Prog. Oceanogr.
,
26
,
179
241
.
Smith
,
L.
,
2002
:
What might we learn from climate forecasts?
Proc. Natl. Acad. Sci. USA
,
99
,
2487
2492
.
Suarez
,
M.
, and
P.
Schopf
,
1988
:
A delayed action oscillator for ENSO.
J. Atmos. Sci.
,
45
,
3283
3287
.
Sun
,
C.
,
Z.
Hao
,
M.
Ghil
, and
J. D.
Neelin
,
2002
:
Data assimilation for a coupled ocean–atmosphere model. Part I: Sequential state estimation.
Mon. Wea. Rev.
,
130
,
1073
1099
.
Tippett
,
M.
,
S. E.
Cohn
,
R.
Todling
, and
D.
Marchesin
,
2000
:
Low-dimensional representation of error covariance.
Tellus
,
52
,
533
553
.
Zebiak
,
S. E.
, and
M. A.
Cane
,
1987
:
A model El Niño–Southern Oscillation.
Mon. Wea. Rev.
,
115
,
2262
2278
.

APPENDIX

Parameter Estimation for the Partitioned Kalman Filter

The partitioned Kalman filter (Fukumori 2002) consists of using a series of reduced-state approximations with physically or statistically independent errors; this partitioning allows one to compute a number of much smaller estimation problems, independent of each other. The key step consists of approximating model state x by a sum of independent substates, x1, x2, . . . , xl, each with a much smaller dimension than the full model state itself:

 
formula

where 𝗕i denotes a transformation approximating particular elements of the model state x by xi. Now we will consider x to be the augmented state, including the parameters to be estimated, and each partitioned substate xi will be a separate, reduced-order approximation of the augmented state x. Here xi may be the amplitudes of various physical modes or approximations corresponding to different parts of a spatial partition of the full domain.

With these slight changes in interpretation, we can now follow the formalism of Fukumori (2002) without any further change. Assuming that the observation error 𝗥 is much larger then forecast errors, the Kalman filter 𝗞 and the error covariance matrix 𝗣 for the full augmented state can be approximated as

 
formula

The forecast error covariance matrix of the partitioned augmented state 𝗣′i is propagated separately from the others by the standard Kalman filter algorithm [cf. Eq. (6)], based on the partitioned state transition matrix 𝗠′i, observation matrix 𝗛′i, and model error 𝗤′i:

 
formula

where 𝗕*i is the generalized inverse of the transformation 𝗕i. Following Eq. (10), the analysis step for updating the single parameter μ is approximated by

 
formula

where 𝗣′μi is the cross-covariance term of the partitioned augmented state available to us from error covariance propagation. This equation effectively combines parameter estimates in distinct partitions and thus naturally accounts for their expectedly different influence on the partitioned state.

Footnotes

* Current affiliation: CSIRO Marine and Atmospheric Research, Floreat Park, Western Australia, Australia

+ Additional affiliation: Geosciences Department, and Laboratoire de Météorologie Dynamique (CNRS and IPSL), Ecole Normale Supérieure, Paris, France

Corresponding author address: Dmitri Kondrashov, Department of Atmospheric and Oceanic Sciences, and Institute of Geophysics and Planetary Physics, University of California, Los Angeles, 405 Hilgard Ave., Los Angeles, CA 90095-1565. Email: dkondras@atmos.ucla.edu