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.
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.
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
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:
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 :
and rewrite our model equations for the augmented system:
The situation of interest is one in which μ itself is not observed, so
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
The analysis step for the augmented system involves only observations of the state:
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
Using the definition of H in Eq. (5), we obtain
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:
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):
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
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 q′n 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.
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.
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).
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).
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.
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.
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.
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.
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.
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.
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, x′1, x′2, . . . , x′l, each with a much smaller dimension than the full model state itself:
where 𝗕i denotes a transformation approximating particular elements of the model state x by x′i. Now we will consider x to be the augmented state, including the parameters to be estimated, and each partitioned substate x′i will be a separate, reduced-order approximation of the augmented state x. Here x′i 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
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:
where 𝗕*i is the generalized inverse of the transformation 𝗕i. Following Eq. (10), the analysis step for updating the single parameter μ is approximated by
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.
* 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: email@example.com