The ensemble Kalman filter (EnKF) has proven its efficiency in strongly nonlinear dynamical systems but is demanding in its computing power requirements, which are typically about the same as those of the four-dimensional variational data assimilation (4DVAR) systems presently used in several weather forecasting centers. A simplified version of EnKF, the so-called ensemble optimal interpolation (EnOI), requires only a small fraction of the computing cost of the EnKF, but makes the crude assumption of no dynamical evolution of the errors. How do both these two methods compare in realistic settings of a Pacific Ocean forecasting system where the computational cost is a primary concern? In this paper the two methods are used to assimilate real altimetry data via a Hybrid Coordinate Ocean Model of the Pacific. The results are validated against the independent Argo temperature and salinity profiles and show that the EnKF has the advantage in terms of both temperature and salinity and in all parts of the domain, although not with a very striking difference.
A review of sequential data assimilation methods with respect to the stochastic systems they solve for can be found in Bertino et al. (2003). The simplest data assimilation method is optimal interpolation (OI) and can be derived within both variational and sequential frameworks. OI has the advantage of simplicity of implementation and low computational cost, and has been applied to many ocean data assimilation studies (Rienecker and Miller 1991; Carton et al. 2000; Etienne and Dombrowsky 2003). However, the OI methods are strongly dependent on an arbitrary choice of background error statistics (Bouttier and Courtier 1999) and are designed for statistical systems. The Kalman filter (KF) is a widespread sequential method for dynamical systems that relies on assumptions of linearity, of an unbiased model and observations, and of initial and model error covariances, instead of the background error covariance matrix. The extended Kalman filter (EKF) can be used in weakly nonlinear systems but instabilities have been reported in stronger nonlinear cases (Evensen 1992).
For high-dimensional systems, handling the error covariances in the KF and EKF methods is practically intractable. The ensemble Kalman filter (EnKF) is an effective data assimilation approach introduced by Evensen (1994) using an ensemble step for the forward integration of error statistics and the traditional update equation of the KF. The EnKF is attractive since it avoids many of the problems associated with the traditional EKF such as instabilities and reduces the numerical needs to a reasonable level (Burgers et al. 1998; Evensen 2003). The EnKF has been employed to assimilate data within a number of different contexts (Evensen and van Leeuwen 1996; Houtekamer and Mitchell 2001; Keppenne and Rienecker 2003), including operational forecasts at the Canadian Meteorological Office, the University of Washington, and in the Towards an Operational Prediction System for the North Atlantic European Coastal Zones (TOPAZ; Bertino and Evensen 2003) system. Leeuwenburgh (2005) has used the EnKF method to assimilate altimetry data in a twin experiment for the tropical Pacific, showing the potential capacity to correct subsurface fields along the equator.
Evensen (2003) proposed an ensemble optimal interpolation (EnOI) method based on a simplification of the EnKF. It is a simple OI scheme where the analysis is computed in the space spanned by a stationary ensemble of model states sampled (e.g., during a long model integration). Since the error statistics are invariant in time, EnOI provides a suboptimal solution compared to the EnKF. On the other hand, the ensemble forecast computation cost of EnKF is N times of that of EnOI if N is the EnKF ensemble size. Typically, in realistic meteorological and oceanic applications, N is approximately 100. This huge difference between the computational costs of the EnKF and the EnOI needs to be justified by the quality of the assimilation results in order to select the EnKF or EnOI in realistic applications where the computational cost is a primary concern. So far, there are very few if any such studies (Oke et al. 2007). The comparison study is also important because the implementation of model errors is still very simple in current applications of the EnKF. In the present application we consider the atmospheric forcing as the only source of model errors. The comparison can help to evaluate this simple implementation of the EnKF and EnOI, then to suggest future improvements for both methods.
In this paper, we compare 3-yr experiments assimilating altimetry data using the EnKF and the EnOI in the Pacific. The assimilation results are verified using independent in situ Argo profiles of temperature and salinity and independent sea surface temperatures (SSTs). Previous validation of the EnKF with a realistic altimetry data assimilation experiment against independent observations can be found for the Atlantic Ocean (Brusdal et al. 2003) as well as the Indian Ocean (Haugen and Evensen 2002). There have been several altimetry assimilation studies of the Pacific Ocean (Verron et al. 2000; Vossepoel et al. 1999; Durand et al. 2003; Leeuwenburgh 2005), but few of them use real in situ observations of temperature and salinity for validation and few experiments extend over a period longer than a year and outside of the equatorial band.
2. EnKF and EnOI
We define ensemble members ψi ∈ ℜn (i = 1, … , N), where N is the ensemble size and n is the dimension of the model state. The model states of the ensemble are held in an ensemble matrix 𝗔:
The ensemble mean is stored in ( = 𝗔N), and its anomaly matrix is
where N ∈ ℜN×N is the matrix in which each element is equal to 1/N.
The observations vector is d ∈ ℜm, and their associated uncertainty is ɛj. We can simulate N random vectors of perturbed observations:
where dj is sampled from a N(d, 𝗥) Gaussian distribution and 𝗥 is the observation error covariance matrix. We store the observation uncertainty ɛj in a matrix γ.
Then, the ensemble of innovation vectors is defined as the matrix
where 𝗛 is the measurement operator that transforms the model state to the observation space and 𝗗 = (d1, d2, … dN).
We can write the forward propagation step for any of the ensemble members as
where k is the time step and qk is a random model error drawn independently for each member from a Gaussian distribution. The operator g is a nonlinear function representing the model integration.
The analysis update equation is given in matrix form:
where 𝗣e is the background error covariance matrix.
Here, the Kalman gain is
Equation (6) can be written using only anomalies as
To compute the denominator efficiently, a singular value decomposition (SVD) is used to calculate the ensemble uncertainties and the measurement errors. Different analysis schemes are detailed in Evensen (2003); we have chosen the original scheme designated as “analysis 2.”
c. Initial ensemble and model errors
The EnKF initial ensemble is generated with random fields following a spectral method (fast Fourier transform) documented in Evensen (1994, 2003) and with more details provided in Natvik (2001) and Wan et al. (2008). At the initial time, we pick a model snapshot after a spinup integration in free-running mode. We assume that the inaccuracy in this snapshot can be attributed to a displacement of the vertical coordinates. We sample the initial ensemble by perturbing the layer interfaces. An improved sampling scheme based on a singular value decomposition of an oversized sample (Evensen 2004) has been used here in order to reduce the redundancy of information in the initial ensemble.
Then, the initial ensemble is integrated forward for 1 month with model errors introduced through the atmospheric forcing fields in order to let the multivariate correlations build up within the ensemble. The model errors are constituted with random red noise added to the atmospheric forcing fields as outlined in Evensen (2003). The spatial correlation is derived from the given spatial scale (also called the decorrelation scale) and typical grid size, while the temporal correlation is formed by a Markov process with a sequence of time-correlated pseudorandom fields with means equal to zero and variances equal to 1. In this study, we take 10 grid cells for the decorrelation scale, 50 km for the typical grid size, and 3 days for the time correlation scale. Equation (9) gives the evolution of the model errors:
where qk is the model error and wk is a white noise process in time. The factor β represents the correlation between two time steps. In addition, Δt is the time step, σ is the standard deviation of the model error, and ρ is a factor to be determined. (In this study, we use ρ = 1.) Detailed discussion about the parameters used here is documented in Natvik (2001) and Wan et al. (2006).
The EnOI analysis is computed in the space spanned by a stationary ensemble of model states (also noted 𝗔), sampled randomly during a long-time integration that includes climatological variability. Only one model integration is required. The computational cost is thus very small compared to the EnKF.
The analysis equation is similar to Eq. (8), as follows:
where α ∈ (0, 1] is the parameter giving different weights for the forecast and measurement error covariances. The variance of a stationary ensemble over a long period usually overestimates the instantaneous variability; therefore α < 1. The EnOI provides a multivariate analysis that has properties that are similar to the EnKF when it comes to the conservation of the dynamical balance (Evensen 2003). The practical implementation is also similar to that of the EnKF but the error statistics do not evolve with time.
3. Model and observations
a. The ocean general circulation model
The Hybrid Coordinate Ocean Model (HYCOM; Bleck 2002) is an upgraded version of the Miami Coordinate Ocean Model (MICOM). HYCOM is a primitive equation, ocean general circulation model with a vertical coordinate that combines the advantages of three vertical coordinate systems: it is isopycnic in the open, stratified ocean and reverts smoothly to Cartesian coordinates (either z levels or terrain-following levels) in unstratified or shallow waters. A sophisticated vertical mixing scheme, the K-profile parameterization (KPP), was included in this study.
The computational domain spans the Pacific Ocean from 28°S–52°N to 95°E–70°W. The standard Mercator horizontal grid configuration used in earlier HYCOM studies has been modified with the conformal mapping tools of Bentsen et al. (1999), and the model grid sizes range from 42 to 72 km in the Pacific, where the sizes are large along the tropical and small in the North Pacific. The vertical direction is split into 22 hybrid layers with reference densities increasing linearly with ⅓ in the upper 10 layers and exponentially with power 2 in the 12 bottom layers. Reference densities range from 18.00 to 27.84 and the three-dimensional state variables are water temperature, salinity, layer thickness, and velocity. The lightest reference densities are used to enforce the use of z coordinates near the surface.
The high-frequency synoptic forcing fields used were temperature, wind, and relative humidity determined from dewpoint temperatures, which were acquired from the European Centre for Medium-Range Weather Forecasts (ECMWF) four times each day. Clouds and precipitation results are based on the climatology of the Comprehensive Ocean and Atmosphere Data Set (COADS) and on Legates and Willmott (1990). At its surface and along the lateral boundaries, the model temperature and salinity are relaxed toward the Generalized Digital Environmental Model (GDEM) climatology (Teague et al. 1990), with a common relaxation time scale of 100 days and over 15 grid cells near the lateral boundary.
b. Experiment design
The state variables that are integrated into the assimilation experiments are all the model prognostic variables: temperature, salinity, current velocity, and layer thickness in the water column, plus the barotropic pressures and velocities. The EnKF initial ensemble is generated following Wan et al. (2008); 100 ensemble members were used, as in Lisæter et al. (2003). The EnOI statistical ensemble is randomly selected from a model run of 20 yr and the size is the same as for EnKF. Table 1 summarizes the parameterization of the initial error, the model errors, and the observation errors. After a 1-month spinup of the initial ensemble, both EnOI and EnKF are run from 1 January 2002 to 31 December 2004. The schematic flowchart of the two experiments is shown in Fig. 1.
A local analysis is used in both experiments to cut off possible long-range covariances and to avoid the problems associated with a large number of the observations. Each grid cell is updated by the M nearest observations within the influence radius r0. In the EnOI experiment, the weighting parameter α is taken as 0.45, which is obtained from other sensitive experiments.
c. The observations
Observational data assimilated in the experiments include sea surface height anomaly (SSHA) satellite maps merged from different missions: the Ocean Topography Experiment (TOPEX)/Poseidon, Jason-1, and the European Remote Sensing Satellite-2 (ERS-2) and the Envisat (Ducet et al. 2000). Datasets are provided by Collecte Localisation Satellites (CLS). The spatial resolution of the datasets is ⅓° × ⅓°. A mean sea surface height (MSSH) is added to the observations and the innovations are computed by interpolating the model SSH to the observation grid. MSSH is calculated from a long-term model integration. The assimilation window is weekly.
The validation data include Argo floats profiles and optimum interpolation sea surface temperatures (OISSTs; Reynolds et al. 2002). OISST data are produced at the National Oceanic and Atmospheric Administration (NOAA) using both in situ and satellite SSTs, mapped weekly on a 1° grid. Temperature and salinity profiles from Argo floats are also used. To compare with the assimilation results, we interpolate the OISSTs to the model grid. The Argo profiles are interpolated linearly to model density layers. All of the Argo float profiles used over the 3-yr assimilation runs are shown in Fig. 2. Several locations can belong to the same float. The distribution of profiles is clearly inhomogeneous with the densest sampling in the northwest Pacific, but we are confident that it is sufficient to validate and compare the different assimilation methods.
For the comparison between the model results and observations, we interpolate the model results to the Argo float positions and interpolate the OISSTs to the model grid in the horizontal direction by bilinear interpolation. In the vertical direction, we interpolate the Argo levels linearly to model density layers.
a. Ensemble statistics
Comparing Eqs. (8) and (10), the analyzed best estimates from the EnOI and the EnKF (the ensemble mean in the case of the EnKF) methods only differ by the ensemble used to calculate the forecast error covariance matrix. The EnOI uses a static ensemble collected from a long model run, while the EnKF ensemble is time dependent and reflects the instantaneous dynamics.
The ensemble standard deviation is calculated by
where σe is the ensemble standard deviation, ϕi is a member of the ensemble, and N is the ensemble size.
To view the correlation pattern, the correlation coefficient rkl is
where ϕi is a member of the ensemble, k and l are the indexes of two model variables at two grid points, and ϕp and ϕl are the corresponding ensemble averages.
Figure 3 shows the standard deviation of each ensemble for surface temperature in the two experiments. To compensate for the large time-averaged variability, the EnOI standard deviation of temperature is weighted by the parameter α in Fig. 3a. Figures 3b–d show three snapshots of standard deviations of temperature in the EnKF experiment on 2 January 2002, 24 April 2003, and 18 August 2004. In Fig. 3a, the large seasonal variability of the North Pacific dominates the forecast error. Along the eastern equatorial Pacific, the standard deviation is similar to that of the EnKF experiment. The model uncertainty of the EnKF experiment is changing with time. In winter, the uncertainties are larger in the southwest Pacific. In spring, they migrate northward to the western equatorial Pacific and later in the summer to the northwest Pacific.
In Wan et al. (2008), it is shown the correlations follow the current patterns very well. We select a point located at 33°N, 137°E on the main axis and near a large meander of the Kuroshio Current. Through the calculation of the correlation coefficient between this point and other points in the domain, we can analyze the variability of the current pattern. We select four snapshots of EnKF correlations with the main current axis in 2004 (Fig. 4). As we know, the main current is expected to be oriented southwest–northeast. However, during the summer of 2004, a larger meander event occurred and lasted for about 1 yr (Usui et al. 2008). In early 2004, the current pattern is narrow, oriented southwest–northeast (Fig. 4a) and then extends to the southeast and northwest (Figs. 4b and 4c). For the EnOI experiment, the correlation is stationary and reflects mostly the seasonal variability but not the current pattern (result not shown). At a relatively low horizontal resolution, it is difficult for the model to reproduce the large meander. However, we still get some information about the changes of the current from the correlation patterns. The correlation of the sea surface height (SSH) of the model in the EnKF experiment (Fig. 5) agrees with the analysis of the current pattern in the same region and the current stream portrays the large meander clearly.
Figures 3 –5 illustrate the potential advantages of the dynamical evolution of the error statistics in the EnKF over the statistical EnOI. In the following, we assess both methods against independent observations.
During the three years, an El Niño episode occurred during the period of late 2002. In Fig. 6, four snapshots were selected to illustrate the correlation pattern of SST located at a point (0°, 135°E) in the eastern tropical Pacific, with two of images taking place during the episode and the others happening during normal episodes. The oceanic Niño indexes (ONIs) of the 3 months are 1.48, −0.12, and 0.24, according to the NOAA Web site. The pattern shows high horizontal correlation during the El Niño episode.
b. Validation against independent data
In this section, we first compare the standard deviation of the ensemble with the root-mean-square errors (RMSEs) against OISSTs to monitor the evolution of the ensemble. Figure 7 displays the results from the 3-yr experiment. For the EnOI experiment, the ensemble is static and its spread is constant, while the ensemble of EnKF is time dependent and the spread fluctuates. In the first months of the experiment, the EnKF spread underestimates the RMS errors, but the two curves converge during the course of the first year. The magnitude of the spread is very close to the RMSE. The initial errors in the stratification are unlikely to play an important part in the sea surface temperatures. The EnKF error statistics for SST are essentially in balance between the uncertainty added by the model errors in the atmospheric fields and their reduction when assimilating the altimetry data. With the parameters chosen here, the combination of the two shows no visible trend over 3 yr, neither growing nor decreasing (“ensemble collapse”), which indicates that the combination generates a sustainable ensemble spread. The time series show a faint seasonal cycle, with the error being slightly larger during winter than summer. This cycle appears both in the EnKF error estimate and in the RMS errors, indicating that the evolution of the model errors in the EnKF is realistic. In the surface, the model errors due to random atmospheric forcing are therefore effective and the variability of the spread shows realistic time-dependent errors.
We compare the temperatures and salinities in different layers of the 3-yr experiments against observations over the whole model domain. The evolution of the RMSE over the 3-yr experiment is shown in Fig. 8. The model surface temperature is compared to the OISST while all other results are compared to Argo profiles. We select three layers, the surface layer (top 3 m in the model); the model layer with potential density 24.00 kg m−3, representative of the bottom of the mixed layer; and the layer whose potential density is 27.02 kg m−3, located below the thermocline. It is noted that the spatial and temporal sampling of the Argo floats is inhomogeneous, in contrast to OISST. This explains why the time series, other than surface temperatures, are slightly more erratic. Furthermore, we have not attempted any declustering of the Argo data, so the RMSE may therefore be weighted toward the northwest Pacific over other parts of the basin and since the northwest is the most dynamically active region, the RMSE against the Argo data might overestimate the actual basin-average errors. The comparisons are however clear enough for evaluating the assimilation methods.
It is obvious that the RMSEs of temperature and salinity at the surface and in the intermediate layer (σ = 24.00 kg m−3) are smaller in the EnKF experiment than in the EnOI experiment, and they are also smaller in the EnOI than in the free run. The comparison holds at nearly all times and at the three selected levels. After assimilating with the EnOI, the RMS errors of the temperature drop by 0.1°C and those of salinity by 0.05 psu in the surface. The reductions are larger with the EnKF, by 0.2°C and 0.1 psu in the surface, respectively. The percentage reduction at the surface from the EnKF experiment is 17.0% for temperature and 33.3% for salinity, exceeding the values of 6.1% and 12.2% from the EnOI experiment. In the intermediate layer (σ = 24.00 kg m−3), the RMSE of the temperature seems to decrease even more than at the surface, but the Argo sampling favors the northwest Pacific where errors are larger. The reduction percentages of EnKF in the intermediate layer are 10.3% for temperature and 17.3% for salinity, while those of the EnOI experiment are 1.8% and 7.2%. In the deeper layer (σ = 27.02 kg m−3), the reductions are smaller, indicating that neither the EnOI nor the EnKF have a strong impact under the thermocline. This could be due to possible degradations of the assimilation imbalances. The reduction percentages from the EnKF and EnOI experiments in the deeper layer are 7.8% and 3.6% for temperature and 4.6% and −0.5% for salinity. The model errors are created by random atmospheric forcing fields. There is no obvious shift (also see in Fig. 7). Though it cannot correct the model enough, the assimilation results do improve. The error reduction is constant over the whole 3-yr period, showing satisfactory convergence for both the EnKF and the EnOI. The EnKF performs better than the EnOI with respect to salinity errors, either due to more realistic correlations between sea surface heights and model salinity fields, or indirectly due to better positioning of the ocean features.
To study the vertical distribution of the RMS errors, we split the domain into five subregions shown in Fig. 2. We select two subregions where Argo observations are dense enough and the characteristics of these subregions are representative to study the vertical distribution of the RMSE. Figure 9 show the profiles of the temperature RMSEs and the salinity RMSEs from the surface to 500-m depth. Region 1, the northwest Pacific, is well covered by numerous Argo floats, while region 4 is an extension of El Niño-3.4. In both regions, the assimilation consistently reduces the RMSEs both at the surface and in deep waters.
Above the thermocline (∼50–100-m seasonal thermocline for region 1 and 200–250 m for region 4), and especially for the top 100-m salinities in region 4, the innovation is larger than that below the thermocline, which shows that altimetry data can improve the control of salinity, in contradiction with Vossepoel et al. (1999). The reason is likely that the ensemble statistics used in both the EnKF and EnOI are multivariate and project equatorial sea level anomalies onto upper 100-m temperature and salinities. The results of the EnKF experiment are overall better than those of the EnOI experiment, especially in deep waters, showing that the EnKF has not only a positive effect in surface waters, directly in contact with the model errors, but also below.
Finally, the mean RMSEs of the temperatures and salinities of different layers averaged over the whole Pacific and for 3 yr in the three experiments are listed in Table 2. In most layers, a reduction in RMSEs is evident for both temperature and salinity, and the EnKF experiment shows larger reductions than does the EnOI experiment.
5. Discussion and conclusions
We have compared two data assimilation methods, the EnOI and the EnKF. Both of these approaches assimilate real altimetry data into a HYCOM model of the Pacific. The experiments spanned 3 yr from 1 January 2002 to 31 December 2004, and a free-running experiment has been run during the same period. The results are validated against the real observations: remote sensing (OISST) and Argo profiles.
The present study confirms the findings of the twin experiments by Leeuwenburgh (2005) but differs by using a real-data assimilation experiment. This study also agrees with Durand et al. (2003) in the ability of the satellite data assimilation to correct ocean fields below the mixed layer as confirmed by independent profiles from the Argo array.
The exploration of a dynamically evolving statistical ensemble in the EnKF shows realistic features that are missing in the EnOI static ensemble. The ensemble statistics in the EnOI experiment near the ocean surface are substantially influenced by the annual cycle in extratropical regions. This influence is absent in EnKF because the temporally evolving statistics are determined over a time window that is much shorter than 1 yr.
Through the comparison of RMS errors in temperature and salinity in different layers, we find that the EnKF yields marginal but consistent improvement over the EnOI in the terms of the temperature and salinity fields at all times and through the entire water column. The improvements made by the EnKF over the EnOI are, perhaps, smaller than expected, but they are still significant with respect to the high sensitivity of the currents to the water mass properties. However, the implementation of the EnKF follows simplifying assumptions, especially the assumption that the atmospheric forcing is the only source of model error, and could be further refined by including uncertainties in other model parameters such as those of the mixing scheme. Moreover, the parameter α is an important factor in the EnOI experiment. Different parameters in different seasons will help ensre that EnOI is implemented well.
The EnKF ensemble spread seems well sustained through the entire EnKF experiment, and fits well with the observed errors against the OISSTs. Since the ensemble spread is a balance between the source of the uncertainties added at the surface boundary conditions and a “sink” of the uncertainties removed by the assimilation of sea surface heights, a balance for the ocean surface seems to have been achieved with the value of the error statistics used in this study. Even though the error statistics applied to altimetry observations and to atmospheric forcing may both be relatively high, the study demonstrates that a nonlinear physical ocean model can well be handled by the EnKF as an observed–controlled system: observed by satellite altimeters and controlled by atmospheric forcing fields. The relative importance of other types of observations and other control parameters can then be evaluated within a similar framework.
This work is supported by the Knowledge Innovation Program of the Chinese Academy of Sciences (Project KZCX1-YW-12-03), the National Basic Research Program of China (2006CB403600), and the National Science and Technology Infrastructure Program (2006BAC03B04). LB is supported by a private donation from Trond Mohn C/O Frank Mohn AS, Bergen, and the MERSEA project from the European Commission (Contract SIP3-CT-2003–502885.). The authors would like to thank two reviewers for their helpful comments.
Corresponding author address: Liying Wan, National Marine Environmental Forecasting Center, State Oceanic Administration, 8 DaHuiSi Rd., Haidian District, Beijing 100081, China. Email: email@example.com