Application of a hybrid EnKF-OI to ocean forecasting

Application of a hybrid EnKF-OI to ocean forecasting F. Counillon, P. Sakov, and L. Bertino Mohn Sverdrup Center/Nansen Environmental and Remote Sensing Center Thormøhlensgate 47, Bergen, Norway Received: 20 March 2009 – Accepted: 3 April 2009 – Published: 15 April 2009 Correspondence to: F. Counillon (francois.counillon@nersc.no) Published by Copernicus Publications on behalf of the European Geosciences Union.

For the GOM application, we believe that the info was already given: "The EnOI uses a historical ensemble that is composed of 122 weekly model outputs over a 2.5 years period". The missing information for the QG experiment has been added (p 3, 3 rd paragraph): "The hybrid covariance blends 200 static members gathered randomly over a period of 500000 model time steps with …. ". Thank you.
The authors need to be clear about the meaning of the term "diverging". I think they use it to refer to the errors -that is, the errors diverge, or grow (see line 26 of pg 661 for eg). But for an ensemble, it could mean the ensemble diverges, or for the application, it could mean that the circulation diverges. I suggest the authors adopt a different phrase. They tend to say "the method diverges". Perhaps they simply mean the system fails.
Thank you. We have replaced the sentence by the following (Page 4, 3 rd paragraph): Note that the hybrid covariance method with beta=1 is failing when inflation (or deflation) is used. re: equation 4 -did you test what happens if you just augment the stationery ensemble with a few dynamic members? This would expand the sub-space of the ensemble (following the arguments on pg 656, following line 18), but would eliminate the need for the "blending factor" beta. Also, assuming the number of dynamic members is always small, c_d will always suffer from sampling error -but augmenting the static ensemble should only act to improve sampling error due to finite ensemble size. I suspect the ratio of the number of static members (n_s) to the number of static + dynamic members (n_s+n_d), (n_s)/(n_s+n_d) _ optimal beta. Note the the "optimal" beta for the HYCOM case is 0.95. n_d = 10, n_s = 122, so (n_s-n_d)/n_s = 112/112 _ 0.92 -pretty close to just augmenting the static ensemble with the dynamic members.
We have tried different methodology on the QG model, (i.e. Dressing EnKF, augmenting the state vector, scaling the alpha locally with the variance of the dynamic ensemble) but they all lead to larger RMSE that the method proposed here.
In the following figure, the dashed black line shows your suggested formula for the case of the QG model. It is different from the optimal values estimated. We suspect that the optimal beta, will always have an exponential decrease, but that the horizontal range of the decrease will be dependent on the model subspace.
Here, it appears that with localization of 25 grid cells, 40 members seems almost sufficient to reach maximum accuracy, which might suggest a model subspace of about 40. I also add a plain black line where instead of ns=200, I am assuming that ns is the model subspace in the formula suggested and that ns=40. This curves is then matching much better the optimal.
The results in Figure 2 are quite unsatisfying. I suggest they be excluded. Figure  3 gives the same story without the uncertainty. The configurations that succeed/fail seem quite random. There's not a very clear trend for configurations that succeed or fail. For example, for EnKF-OI (beta = 0.8) the system works ok when the inflation is <1.2, equals 1.45 or 1.55, but fails if inflation is 1.25-1.4, and fails for 1.5 and 1.6. Perhaps the experiments weren't run for long enough, so the statistics still suffer from sampling error. Figure 3, when more dynamic members are used, doesn't have the same problem.
For 5 members, although a value of beta=0.7 shows the minimum error, a value of beta=0.9 gives almost as accurate result and is more stable, as all the experiments are converging (independently of the scheme or of the innovation used). In addition this values would better match the exponential curve. But we agree with your comment and the Figure has been removed.
On pg 662, line 18 the authors state that "The flow dependent covariance is expected to be more consistent than static covariance." I don't agree with this. Perhaps the flow dependent covariance could be more optimal, but I'd expect it to be less dynamically consistent. Recall that the static covariance is generated from a long run that is not frequently updated, or initialised. The static ensemble is therefore dynamically consistent. By contrast, because the dynamic ensemble is regularly updated, or initialised, it includes some noise associated with the shock of the initialisation. The dynamic ensemble is therefore probably less dynamically consistent than the static ensemble. I therefore don't accept the explanation in the paragraph following line 18 of pg 662. Perhaps this is why these results contradict those of Wang et al. (2007) -see line 28 of pg 662.
We disagree with some of your statements. It is true that small ensembles suffer from sampling error, however historical ensemble assumes that temporal variability is representative of instantaneous forecast error. When the dynamic ensemble is large enough, the sampling error becomes negligible. On the contrary, the assumption made in the EnOI can remain wrong regardless of the ensemble size. For example, if there are two distinct processes, occurring at different time of year, then the EnOI will provide a combination of the two processes as an optimal estimate. However, we agree that the sentence was not clear and we have chosen to reformulate it (Page 4, last paragraph): "When the dynamic ensemble is large enough, the sampling error becomes negligible. The EnOI makes the additional assumption that historical ensemble is representative of instantaneous forecasting error. Therefore, when dynamical ensemble is self-sufficient, one merely expects the adjunction of static covariance to deteriorate the results..." Assimilation shocks in our case are small and does not last more than 2 days even when alpha is purposely strengthen (see Counillon et al. 2009b). Our perturbation system does not produce much noise neither (only small noise has been observed on the first days when perturbing the lateral boundary condition, Counillon et al 2009a). We are still confident that the reason Wang et al. didn't reach as accurate result with ETKF-OI than ESRF, is because they do not use localization nor the hybrid covariance for updating the ensemble anomalies.
On pg 672, line 15, the authors say that as the size of the dynamic ensemble is increased, the optimal beta (blending coef) decreases -that is the dynamic ensemble is given more weight. But the authors also find that the hybrid system always outperforms the dynamic system. How can these conclusions be reconciled?
To us these two conclusions are not contradictory. The optimal beta reflects the diminishing returns from additional dynamical members as the size of the dynamical member increase. We believe that if ensemble size tends to infinity, the optimal parameter beta will converge to 0 remaining positive.
Minor comments: -the statement on pg 656, line 10, that the cost of the assimilation will remain negligible is not accurate. This depends on many factors that vary with each applications -the number of obs assimilated, the localisation radii, the model resolution etc. This statement should be removed. We softened our statement including the term often (page 2, 2 nd paragraph): "For ocean applications, observations are typically less frequent than for the atmospheric applications, so that the model integration step often dominates the computational cost relative to the assimilation step. In this case, an update of both the ensemble mean and covariance using the hybrid covariance will lengthen the assimilation step, but will remain negligible with respect to the total computational cost." Thank you.
-The captions for Table 1 and 2 need to include the application's they are valid for. ie Table 1 for QG model, Table 2 for HYCOM.
Thank you -it's meaningless to include labels 5.2, 5.4 etc on fig 7b. The values for the horizontal axis are run numbers. For clarity, you could extend the horizontal range to 4.5-7.5. Thank you. This is correct, we have removed the meaningless labels, but we did not feel necessary to extend the horizontal range. Thank you.

Introduction
Data assimilation methods can use ensembles to obtain and propagate the system state and the background error covariance. Two different approaches are often used. The first one, referred to as ensemble optimal interpolation (EnOI, , uses a static ensemble of model states. The second and theoretically more consistent approach, uses a dynamic ensemble, as for example the ensemble Kalman filter (EnKF, Evensen, 2006). Dynamic ensembles provide a flow-dependent background error covariance, but they can require of the order of 100 model realizations for realistic oceanic applications (Natvik and . Therefore, in practice one has to either favor the high model resolution combined with an inferior data assimilation method, or a more optimal data assimilation method at the expense of the model resolution. However, in many applications it is important to have a sufficient model resolution for obtaining a realistic representation of the dynamics. For example, Chassignet et al. (2005) show the importance of the model resolution for placing accurately the fronts in the Gulf of Mexico. Yin and Oey (2007) and later Counillon and Bertino (2009a) investigate ensemble forecasting with a small ensemble (10 members) and a high resolution model of the Gulf of Mexico. Yin and Oey (2007) show that a probabilistic forecast provides a better accuracy than a single forecast, and Counillon and Bertino (2009a) show using an advanced perturbation system that the ensemble spread is correlated in space and time with the model error. This indicates that even small dynamic ensembles can be useful for data assimilation purposes. Hamill and Snyder (2000) suggest a hybrid scheme called EnKF-3DVAR that combines the covariance from a dynamic ensemble with the static background covariance from 3DVAR. Each ensemble member is updated variationally with perturbed observations. The method is tested with a quasi-geostrophic model, and it shows improvements rela-tive to 3DVAR. The improvements are largest in case of a sparse observation network. A similar approach was also suggested by Lorenc (2003). Etherton and Bishop (2004) suggest a hybrid scheme called the ETKF-OI. It uses the ensemble transform Kalman filter (ETKF, Bishop et al., 2001) analysis scheme for updating the dynamic ensemble. In order to reduce the computational cost, only the ensemble mean is updated with the hybrid covariance, whereas the ensemble is updated using the ETKF. The ETKF-OI is shown to outperform the 3DVAR in a two-dimensional turbulent model. Wang et al. (2007) use the ETKF-OI, with localization for the ensemble mean, and compare it to the ESRF with localization (Whitaker and Hamill, 2002) for a two-layer primitive equation model. The ETKF-OI outperforms the ESRF for small dynamic ensemble size (5 members), produces similar results for intermediate dynamic ensemble size (20 members), but is outperformed by the ESRF for larger dynamic ensembles. Finally, Wang et al. (2008a,b) demonstrates that the benefit of hybrid covariance remains when assimilating real observation on a coarse model.
For ocean applications, observations are typically less frequent than for the atmospheric applications, so that the model integration step often dominates the computational cost relative to the assimilation step. In this case, an update of both the ensemble mean and covariance using the hybrid covariance will lengthen the assimilation step, but will remain negligible with respect to the total computational cost. If using the hybrid covariance is beneficial for an update of the ensemble mean, it should also be beneficial for updating the ensemble covariance. Furthermore, an update of the ensemble mean and covariance is more in line with the Kalman Filter. To analyze the benefit of the hybrid covariance, we compare the hybrid covariance EnKF (called hereafter EnKF-OI) to the EnKF, and the hybrid covariance ESRF (called hereafter ESRF-OI) to the ESRF. The ESRF is a deterministic formulation of the EnKF, and yield to better performance than EnKF for small ensemble size (Whitaker and Hamill, 2002).
We see the main reason for using the hybrid covariance is expanding the subspace of ensemble anomalies produced by a small dynamic ensemble. While systems based on the ensemble Kalman filter can yield a theoretically optimal update, in practice this can only happen with an ensemble of a sufficient rank to span the system error subspace. With an ensemble of insufficient rank, the analysis becomes not only suboptimal, but degenerative, resulting in a collapse of the ensemble. In contrast, a theoretically suboptimal EnOIbased system is often able in similar circumstances to yield meaningful updates and to provide an overall robust data assimilation system. The hybrid covariance can be therefore viewed as a compromise between a theoretically superior but computationally expensive data assimilating system based on a dynamic ensemble, and a computationally cheap and robust but theoretically inferior system based on a static ensemble.
A primary objective of our study is to test the hybrid covariance when localization and hybrid covariance are applied to both the ensemble mean and ensemble covariance, by doing a clear cut comparison for different analysis schemes (EnKF, ESRF). A secondary objective resides in applying the hybrid covariance to a realistic application and to an ocean application.
The outline of this paper is as follows. The hybrid covariance method is presented in Section 2. The method is then validated on a simple 1.5-layer reduced gravity quasigeostrophic model in Section 3. Finally, we demonstrate the benefit of the hybrid covariance for a realistic application in the Gulf of Mexico, in Section 4, and present our conclusion in Section 5.

Hybrid covariance methodology
In sequential data assimilation, the system error covariance is often calculated from an ensemble of model states. With the EnOI, the static ensemble A s and the centered static ensemble A ′ s are defined as: where ψ is a model state vector, N s is the size of the static ensemble, n is the size of the model state vector, and the overbar denotes ensemble average. The square brackets denote a horizontal concatenation of the matrices separated by a comma. A static ensemble may contain the model states sampled from a long model integration. The integration should be long enough to contain a wide variety of possible model states.
Similarly A ′ d ∈ R n×N d represents the centered dynamic ensemble matrix, where N d is the size of the dynamic ensemble.
The ensemble covariance matrix calculated from the static ensemble is denoted C s , and the one calculated from the dynamic ensemble C d . They are both assumed to represent the forecast error ǫ: and The superscript T denotes a matrix transpose. The variance of the static ensemble is usually different from the forecast error variance, so that a scaling factor α is introduced in the traditional EnOI framework . The parameter α is part of the tuning of an EnOI system, and is kept constant in the following. As suggested by Hamill and Snyder (2000), we compute a linear combination of the two covariance matrices with an adjustable blend parameter β 1 : 3 When β=0 (resp. β=1) the hybrid covariance is expressed entirely by the dynamic ensemble (resp. static ensemble). Wang et al. (2007) manipulate the matrix CH T ∈ R n×m , with H being the measurement operator relating the prognostic model state variables to the measurements, and m the number of measurements. However, for a large number observations CH T remains a large matrix. Here the matrix A is the matrix of combined ensemble anomalies: The Kalman filter equation is then solved as: where is the Kalman gain, R is the observation error covariance matrix, I is the identity matrix, d is a vector of measurements. An ensemble data assimilation system gets suboptimal due to the limited ensemble size and partially inadequate priors assumptions. Such suboptimalities can lead to an excessive reduction of the ensemble spread, which can be maintained pragmatically by multiplying the assimilation anomalies with a term δ called the ensemble inflation. The superscript "a" refers to the analysis state and "f" to the forecast. Different schemes can be used for solving Equation 8. In the following we apply the two most widely used: the EnKF and the ESRF. The EnKF is based on a Monte Carlo sampling of Equation 7 and applies perturbations to the observations that can impair the stability of the results for small ensembles. To circumvent this problem, the deterministic ESRF solves the analysis Equations 7 and 8 without perturbing the observations, and can provide more accurate and stable results (Whitaker and Hamill, 2002).
The cost of the assimilation time step is usually in O(N × m × n). When the hybrid covariance is used, N = N s + N d , and the time needed for assimilation becomes longer than both the dynamic method and the static method. However, usually N d ≪ N s , so that the cost remains similar to the EnOI.

A quasi-geostrophic model
In order to analyze the capability of the hybrid covariance, we first apply it to a simple 1.5-layer reduced-gravity quasigeostrophic (QG) model with double-gyre wind forcing and biharmonic friction. It is a non-linear model with dimension (127×127) (see Figure 1), and model subspace dimension of order of 10 2 -10 3 . The model is eddy resolving as it generates eddies of size O(10) of the model grid (see Figure 1). More details about the model are given in Sakov and Oke (2008). The model is run over 1000 model time steps, and is assimilating 300 sea surface height (SSH) observations every 10 time steps. The observations are extracted from a model run with lower viscosity, to which white noise is added with a variance 2 of 4. The observations are distributed uniformly over the domain, with a different random offset for each assimilation in order to mimic the typical distribution of satellite tracks (represented by dots on Figure  1). Both the model code and the framework of the data assimilating system used in the experiments are available from http://enkf.nersc.no/Code/EnKF-Matlab.
The parameter α is fixed at a value for which the EnOI is robust and performs with lowest error 3 (α=0.04). This value is also optimal with the hybrid covariance method (not shown). The localization is applied with a Gaussian localization function to the state error covariance matrix by means of a Schur product (Houtekamer and Mitchell, 2001). The localization radius is set to 25 grid cells. The localization and the assimilation time step were both chosen large enough to challenge the data assimilation method. For comparison Sakov and Oke (2008) assimilate every four time steps, and show that a smaller localization radius (of approximately 15 grid cells) provides more stable and accurate results.
The hybrid covariance blends 200 static members gathered randomly over a period of 500000 model time steps with an increasing number of dynamic members until the performance in terms of RMSE saturates. Four data assimilation schemes are compared here: the EnKF, the ESRF and the hybrid covariance EnKF-OI and ESRF-OI.
The accuracy of the system is assessed by the time average root mean square error: ε is calculated from the ensemble mean before assimilation and start at iteration p 0 (here p 0 =10 in order to remove the data assimilation spin-up time). p f is the total number of assimilation steps (here p f =100). ψ t represents a known true field, from which synthetic observations are extracted. It is similar for all runs, and the random seed is fixed, so that all runs of similar ensemble size use the same random perturbation.
We vary the three adjustable parameters of the method, the linear blending coefficient β, the inflation parameter δ and the size of the dynamic ensemble N d , and evaluate the resulting error ε. The parameter β is chosen from 0 to 1 with an increment of 0.1, the parameter δ from 1 to 1.6 with an increment of 0.05, as shown in Figure 2 where every circle represents the run average forecast error ε calculated for a given β, δ and N d . The values of N d are chosen as reported in Table 1. The optimal setting minimizes the error ε. The values at optimum are denoted ε * , β * , δ * , and are reported in Table 1.
When β=0, the EnKF-OI (resp. ESRF-OI) coincides with the EnKF (resp. ESRF) method. When β=1, every member of the dynamic ensemble is updated by the static covariance. The best guess is still provided by the ensemble average in this case, so that the method does not strictly coincide with the EnOI. In these experiments we observe that the EnKF-OI and ESRF-OI with β=1 provide a slightly better estimate than the traditional EnOI, which is in agreement with the results of Wang et al. (2007). This is a consequence of the nonlinearity of the QG model, because in a linear model, the integration of the ensemble mean coincides with the mean of the integrated ensemble. Note that the hybrid covariance method with β=1 is failing when inflation (or deflation) is used.
We start our experiment with 5 dynamical members. It is characterized by a big proportion of failed runs. However, there is a clear core of runs that complete for hybrid covariance methods (EnKF-OI and ESRF-OI) with high values of β, whereas neither the EnKF nor the ESRF (i.e. with β=0) have completed. This indicates that the hybrid covariance methods avoid divergence in some configurations with a small dynamic ensemble, in agreement with Wang et al. (2007). Using the hybrid covariance also reduces the error ε by approximately 26% compared to the EnOI. The hybrid covariance methods become more stable with 10 dynamic members, as fewer runs fail. At the same time, only a single run have complete with the ESRF, but it is of similar accuracy to the EnOI. For 15 members (see Figures 2), all of the four schemes converge. When more than 25 dynamic members are used, the benefit from using the hybrid covariance over the EnKF/ESRF is only slight, and the best ε * is loosely defined on a broad range of values of β and δ.
When the dynamic ensemble is large enough, the sampling error becomes negligible. The EnOI makes the additional assumption that historical ensemble is representative of instantaneous forecasting error. Therefore, when dynamical ensemble is self-sufficient, one merely expects the adjunction of static covariance to deteriorate the results, therefore in this case β * =0. In Table 1, this occurs two times: for 30 dynamic members with the EnKF-OI and for 40 members with the ESRF-OI. The first occurrence with the EnKF-OI is likely to be the result from random variations, as β * is again positive with a larger dynamic ensemble. The second occurrence with the ESRF-OI seems more reliable because the ESRF with an ensemble of 40 members does converge to its maximum accuracy. After refining the discretization in β and δ, we obtain β * =0.05. It indicates that even when the EnKF has nearly converged to its maximum accuracy, the use of hybrid covariances is still beneficial (see Figure 3). This seems to contradict Wang et al. (2007) where the ESRF gives better result than the ETKF-OI for a large dynamic ensemble. However, their result may be caused by the difference between the ESRF with full localization and the ETKF with localization only for updating the ensemble mean.
In Table 1 and Figure 4, the relationship between the blend parameter β and the size of the dynamic ensemble N d is analyzed. The value of β * decreases with increasing N d . This result seems natural, as with an increase of the dynamic ensemble size, the need for a static ensemble is reduced. In Figure 4, the curve of β * for the ESRF-OI is relatively regular. The relationship is noisier with the EnKF-OI due to its Monte Carlo nature, but overall the curves match, and can be fitted by an exponential (green line on Figure 4). This indicates that the parameter β * may be independent of the analysis scheme chosen. We suggest that the value of β * reflects the balance between the rank of the dynamic ensemble and the model error subspace dimension, which implies that one can predict the optimal value of β * for a given size of dynamic ensemble, and the dimension of the model error subspace.

A realistic application: Gulf of Mexico nested model
The dynamics in the Gulf of Mexico (GOM) are dominated by the northward Yucatan Current flowing into a semienclosed basin. This current forms a loop (called the Loop Current, LC) and exits through the Florida Straits. At irregular intervals (Vukovich, 1988;Sturges and Leben, 2000) the LC sheds large eddies that propagate westward across the GOM. The eddy shedding involves a rapid growth of non-linear instabilities (Cherubin et al., 2005b), and these are difficult to forecast (Chassignet et al., 2005;Counillon and Bertino, 2009b). In particular, Eddy Yankee (2006) was problematic for the offshore industry operating in the northern shelf of the GOM. Thus, it is used for investigating the benefits of hybrid covariance for a realistic system conducted below. Chassignet et al. (2005) demonstrate the skill of HYCOM for the GOM. They emphasize the importance of both horizontal resolution and lateral boundary conditions for accurate prediction of the fronts there. A nested configuration can satisfy Counillon François: Application of a hybrid EnKF-OI to ocean forecasting 5 these two requirements with reasonable computing cost. The TOPAZ3 system provides lateral boundary conditions to a high-resolution model of the GOM (Chapter 15 in Evensen, 2006) using nesting techniques described in Browning and Kreiss (1982).

Experimental setup
TOPAZ3 is a real-time forecasting system for the Atlantic and Arctic basins with a configuration of HYCOM, capable of monitoring the circulation patterns in the Atlantic. The grid is created using a conformal mapping of the poles to two new locations by the algorithm outlined in Bentsen et al. (1999). The horizontal resolution varies from 11 km in the Artic to 18 km near the Equator (approximately 1/8 • ). The model is initialized from the GDEM3 climatology (Teague et al., 1990) and spin-up for 16 years. The TOPAZ3 prototype used in this work only provides the boundary conditions and does not include data assimilation, because of computer costs.
The high-resolution GOM model is set-up with 5 km horizontal resolution, which is sufficient to resolve the mesoscale features considering the Rossby radius in the area (R o ≃ 30 km). The model uses a fourth order numerical scheme for treating the advection of momentum in the primitive equations (Winther et al., 2007). To minimize the spin-up time, the initial state is interpolated from an equilibrium state of TOPAZ3, and spin-up for 3 years.
For both models, there are 22 layers in the vertical. The bathymetry is specified using the General Bathymetric Chart of the Oceans database (GEBCO) with 1' resolution, interpolated to the model grids. The models are forced by the 6-hourly and 0.5 • analyzed fields from the European Center for Medium range Weather Forecasting (ECMWF, see http://www.ecmwf.int). The system is described in more detail in Counillon and Bertino (2009b).
At the time of writing the article, we can afford a dynamic ensemble of 10 members in real-time. Here we compare a 10-member EnKF-OI, a 10-member EnKF, and an EnOI by assimilating altimetry data. The EnOI uses a historical ensemble that is composed of 122 weekly model outputs over a 2.5 years period, without data assimilation. All methods use a fixed localization radius 4 of 150 km. A smooth transition is ensured at the edges of the localization area by applying a weight function to the innovations. This function depends on the distance between the observation and the target point (Counillon and Bertino, 2009b). A value of α=0.23 is found to be optimal for both the EnOI and the EnKF-OI set-up at the nowcast forecast horizon, and has is used in the experiments below. In order to identify the optimal value of β, we run the EnKF-OI with the values shown in Table 2.
The dynamic ensemble for the EnKF and the EnKF-OI is spin-up in order to provide realistic correlation. The initial set of restart files is randomly extracted from an a three weeks EnOI run that precedes the starting date, in order to re-duce the spin-up time of the ensemble. The ensemble spread is then maintained by using perturbations and inflation over the successive assimilation cycle. In Counillon and Bertino (2009a), perturbations are applied to the atmospheric forcing, the lateral boundary conditions, and the assimilated state for the GOM model. The perturbations of the assimilated state appear to control the main position of the large features of the GOM (e.g. LC, eddies), whereas perturbations of the boundary conditions (atmospheric and lateral) control the growth of cyclonic eddies at their boundary. The same technique is used here for perturbing the boundary conditions, but the assimilated state is perturbed similarly to the EnKF; i.e. by perturbing the measurements during assimilation. The perturbations of the atmospheric fields (wind stress, air temperature) are simulated with a spectral method , using a 50 km decorrelation radius, which is too small to be represented in the ECMWF data, but large enough to stimulate the ocean model. The perturbations of the lateral boundary conditions are simulated by introducing a time lag to the boundary conditions, which is set randomly at the beginning of each model run, between -37 and +37 days from the actual date. This perturbation system introduces an initial barotropic imbalance, which is damped within a day. In addition, we use an ensemble inflation of 15 %, which increases the ensemble spread and slightly improves the accuracy over time.
The assimilated SLA maps are provided by Ssalto/Duacs on a 1/3 • Mercator grid (Traon et al., 2003), with a one week delay. The standard deviation of the measurements is assumed to be constant, and it is using the average value specified by the provider in the GOM area (3 cm). The measurements are less accurate in the coastal area, therefore measurements are selected only in regions deeper than 300 m, which correspond in the GOM to distances of at least 50 km from the coast. Accordingly, a Gaussian covariance with a decorrelation radius of 50 km is used for the observation error. As the SLA needs to be referred to a mean SSH, a two-years average of TOPAZ3 SSH is interpolated to the high-resolution grid. Qualitatively, it compares well with the mean dynamic topography based on satellites and in situ measurements (Rio and Hernandez, 2004).
The experiment starts six weeks prior to the shedding of Eddy Yankee, in order to spin-up the perturbation system (as in Counillon and Bertino (2009a)). Seven runs are processed with weekly assimilation of sea level anomaly (SLA) and are hereafter referred to as Run 1 to Run 7 (see Figure 5). The first assimilation is applied on the 7 th of June and the last one on the 19 th of July 2006. Each run is integrated for 7 days, which corresponds to a nowcast with respect to the availability of the near real-time altimeter data.

Forecast Errors
The benefit of the hybrid EnKF-OI over the EnOI and the EnKF is analyzed by calculating the RMSE between the model daily average SLA with the daily SLA data maps. The calculation starts from the 5 th of July (i.e. Run 5) to allow sufficient ensemble spin-up prior to the pre-shedding (with the non-linear growth of cyclonic eddies), the shedding and the near-reattachment of Eddy Yankee. Furthermore, this starting date corresponds to the start of available daily altimetry maps. During the shedding event, the fast dynamics of the eastern GOM contrast with the slower activity of the remaining domain. In order to clearly identify the benefit in the shedding area, the statistics are calculated over a restricted area such that the longitude is within [91 • W ; 83.5 • W ] and the latitude north of 22.25 • N. The areas shallower than 300 m are also excluded from the statistics due to inaccurate measurements.
In order to analyze different aspects of the performance of the data assimilation methods that use the hybrid covariance, we form the following characteristics of the system errors: -The overall average error ε 1 , for obtaining the optimal value of β * (see Table 2).
-The average of the run error ε 2 , for analyzing the stability of the methods over the successive runs (see Figure  6 b).
-The error average at a given forecast horizon ε 3 , for characterizing the persistence during model integration (see Figure 6 a). In particular, the error average at the nowcast is reported in Table 2.
The EnKF with 10 members yields a larger ε 1 than all the other schemes (see Table 2). This result is expected because the dynamic ensemble has too few members. In particular, we note that the error ε 2 (Figure 6 b) should decrease in the last run because the eddy has almost reattached to the LC, and the complexity of the prediction is reduced. This does not occur for the EnKF system, for which the run error average ε 2 increases with successive runs. The error for the EnKF-OI with β=0.2 also does not decrease for the last run. This indicates that with 10 members, the EnKF might diverge, and that a minimum value of β >0.2 is necessary for reaching stability with the hybrid covariance method.
In Table 2, there is a clear trend indicating that high values of β are preferable, but the minimum is once again loosely defined, possibly due to the random perturbations of observations and boundary conditions. However, the three statistics ε 1 , ε 2 , and ε 3 indicate that a value of about β = 0.95 is optimal for the hybrid EnKF-OI. On average over the three runs, the EnKF-OI with β=0.95 reduces the total error average ε 1 by approximately 10 %, and the nowcast error average ε 3 by approximately 14 % compared to the EnOI. Although the blending parameter is largely on the EnOI side, the gain in accuracy seems important.
In Figure 6 a, the error ε 3 is doubling during one week of model integration with all schemes. However, one can observe that the EnOI is diverging faster than the run with the hybrid covariance method. Counillon and Bertino (2009a) show that a high value of α allows for higher accuracy initially but diverges faster with model integration. A lower value of α tested for both the EnKF-OI and the EnOI leads to more stable but less accurate predictions at the nowcast stage.
The EnKF-OI with β = 1 is performing better than the EnOI, as shown with the QG model. It is more accurate to 5% on average, and to 8% at the nowcast stage to the EnOI (see Table 2). This indicates that imbalances caused by the perturbation system are too low to deteriorate the performance, and that the sources of variability are efficient.

Frontal analysis
Our main objective is to represent accurately the position of eddies (orientation, center, shape) during the shedding events. Large velocities are located at the fronts, at the outer edge of eddies; with maximal values where the cyclonic and anticyclonic eddies interact. In this Section, we conduct a frontal analysis for the EnOI and for the EnKF-OI (with β = 0.95) in order to qualitatively interpret the statistical gain observed above. Chassignet et al. (2005) show that ocean color (OC) data is useful for identifying the position of the fronts of the LC and the eddies. Furthermore, it provides an independent source of validation, and has higher resolution than the altimetry data. In Figure 7 a the deep blue contour area represents the low chlorophyll water (< 0.3 mg/m 3 ) that originates from the Caribbean Sea. The light green areas characterize the water with higher chlorophyll concentration (> 0.5 mg/m 3 ), which is usually found in areas of high biological production along the coast or within cyclonic eddies. The high chlorophyll water often propagates along the outer edge of the LC eddies, and clearly defines the front.
The position of the front is analyzed for the last three runs at the nowcast stage (represented with the letters a-c in Figure  5) for the EnKF-OI with β = 0.95 and for the EnOI. This period covers the shedding and the near reattachment of Eddy Yankee (2006). On the 12 th of July (Figure 7 a) there is a pre-shedding scenario with a deep cyclonic penetration from the east. On the 19 th of July (Figure 7 b), Eddy Yankee is on the verge of being shed from the LC. On the 26 th of July (Figure 7 c) Eddy Yankee has shed, rotated around a cyclonic eddy located on its eastern side, and started to reattach from the east.
In Figure 7, the ensemble fronts from the EnKF-OI are represented by pink lines over the OC data. The EnKF-OI ensemble mean provides the best guess and is represented by the thick red line. The front calculated from the EnOI is indicated by the thick yellow line. The front position derived from altimeter SSH maps (hereafter referred to as SSH data) is represented by a thick black line. Although OC data has a higher resolution, the SSH data may also be a useful indicator of the error in the assimilated SSH data, or be used to 7 locate the front when clouds mask the OC data. The fronts of the model and SSH data are the 10 cm SSH isolines, which appears to fit best with the front from OC data.
On the 12 th of July, the shape of the front obtained with the EnKF-OI is better positioned than that obtained with the EnOI. In particular, the eddy sheds too early with the EnOI. Furthermore, the northern front is slightly too far to the north. One can also notice that both nowcasts are relatively accurate with respect to the altimetry data.
On the 19 th of July, both methods misrepresent the northern front. Indeed in the OC data, the shape of Eddy Yankee is deformed by a cyclonic eddy that is squeezed between its northern front and the northern shelf. The reason is that this cyclonic eddy was poorly represented in the gridded SLA data assimilated on the 12 th of July, because of its proximity to the coast. In the rest of the domain, the EnKF-OI describes the remaining connection of Eddy Yankee with the LC, whereas the eddy appears as already shed in the EnOI.
On the 26 th of July, the EnOI has placed Eddy Yankee slightly too far to the west, and the EnKF-OI better represent its northern front. The "comet" shape of the eddy seems to be better represented in the EnKF-OI than in the altimetry data.
One can notice that the eastern LC base is systematically too broad in the model, with both assimilation schemes. This may be caused by a model bias. Overall, the EnKF-OI improves the representation of the fronts over the EnOI during the shedding of Eddy Yankee.

Hybrid correlation
In order to obtain a better understanding of the benefits from the hybrid covariance, the correlation of SSH from the static, dynamic and hybrid ensembles are analyzed below. In Counillon and Bertino (2009b), a statistical analysis of the static ensemble throughout the GOM is conducted at two typical locations, one in the center of the basin and one on the upper-shelf. The correlation is found to be realistic in the center of the basin, but has some limitations near the uppershelf area. Here, we analyze the spatial correlation of SSH at the same locations, marked with letters A and B in Figure 8. This analysis is repeated for the static, dynamic, and hybrid ensembles (with β=0.95). Figure 8 shows the variance of the static ensemble and that of the dynamic ensemble. As expected, the variance of the static ensemble shows a maximum value induced by the resident LC base, and a positive track induced by the passage of eddies that drift westward. However, the variance at locations closer to the shelf is very small, in particular at point B. The dynamic ensemble shows a large variance near the front of Eddy Yankee, because a small displacement of a front induces a large difference in the SSH. The dynamic ensemble is expected to be most useful in the dynamic areas and in the upper-shelf coastal area, but its contribution in the middle of the basin should be small because the variance from the static ensemble dominates there.
In the center of the basin (point A in Figure 8), the correlation of SSH and velocity from the static ensemble clearly reflects the geostrophy, as a higher SSH increases the anticyclonic currents around the target point (see Figure 9), and there is a slight anisotropy due to the interaction with the resident LC. The SSH correlation from the dynamic ensemble highlights the position of the southern front of the eddy, but it is noisy and uneasy to interpret. In total, the hybrid correlation is nearly identical to the static ensemble.
In the coastal area, (point B in Figure 8), the correlation of SSH and velocity from the static ensemble is dominated by the seasonal variability. An increase of SSH at the target point increases the SSH almost uniformly in the whole area, and induces a weak shelf current (see Figure 10). In the dynamic ensemble, there is a positive correlation with a small anticyclonic ring, and with a vortex of opposite sign at its boundary. This correlation seems realistic because small anticyclones and cyclones 5 interact with the shelf there (Hamilton and Lee, 2005), but it is noisy due to a small number of dynamic members. In the hybrid ensemble, the correlation with the anticyclonic ring is clearer than in the static ensemble, both for the SSH and currents. Furthermore, the hybrid correlation is smoother than that of the dynamic ensemble.
The benefits of using hybrid covariance compared to the static covariance may not appear as large, but they are still important as they occur at locations where the static covariance is known to be inadequate.

Conclusions
This study follows the work from Hamill and Snyder (2000); Lorenc (2003); Etherton and Bishop (2004); Wang et al. (2007Wang et al. ( , 2008a and analyzes the benefits of using the hybrid covariance. The hybrid covariance combines the covariance from a dynamical ensemble and from a static ensemble. Unlike previous works, the hybrid covariance and localization are applied to updates both the ensemble mean and the covariance. We also evaluate the hybrid covariance with the EnKF and the ESRF analysis schemes instead of the ETKF (Etherton and Bishop, 2004;Wang et al., 2007), or of a variational approach (Hamill and Snyder, 2000). As a consequence, we found the hybrid covariance to be beneficial compared to the dynamical approach; even when a large dynamical ensemble is used. In addition, the hybrid covariance methods are applied for the first time to a realistic application (i.e. assimilating real observations, using a state of the art model at a resolution capable of resolving the dynamics, and thus applicable in an operational setting), and to the ocean.
The results with the QG model indicate that: 5 with a radius smaller than 75 km 8 Counillon François: Application of a hybrid EnKF-OI to ocean forecasting -Using the square-root scheme improves the stability and the accuracy compared to the EnKF. However, the conclusions regarding the hybrid covariance are independent of the analysis scheme used.
-The hybrid covariance method provides better results than both the static ensemble and the dynamic ensemble methods at equal dynamic ensemble size. The benefits are large for a small dynamic ensemble, and small for a large dynamic ensemble.
-The hybrid covariance allows a smooth transition between the EnOI and the EnKF (or ESRF).
-The optimal linear blending coefficient β * decreases exponentially with the size of the dynamic ensemble. In other words, more weight is being put towards the dynamical ensemble when its size increases.
In a realistic HYCOM configuration of the GOM, 10 dynamic members are insufficient for the EnKF to converge. The hybrid covariance EnKF-OI improves both the SSH RMSE (by 14 %) and the position/shape of the fronts compared to the EnOI at the nowcast stage. Using a comparison with ocean color, we show that the improvements are related to a better representation of the fronts. However, the overall results are not as encouraging as with the QG model. It is clear, considering the very high value of β, that a further increase of the number of members should be largely beneficial. In two characteristic locations studied, the correlation calculated from a hybrid ensemble seems more realistic than that of the static and of the dynamic ensembles, and it is possible that an increase of the localization radius can further improve the accuracy of the results. Furthermore, Hamill and Snyder (2000) indicated a larger benefit when a sparse observation network data were assimilated. It seems probable that the benefit from the hybrid covariance could increase if SLA track data were assimilated instead of the interpolated maps.
The results obtained from both models indicate that the hybrid covariance methods can improve both the accuracy and the stability of the system. For 10 members, the optimal value β * is larger with the GOM model than with the QG model, which possibly indicates that the GOM model has more degrees of freedom. It could be interesting to extend the study for an increasing number of dynamic members for the GOM application when more computing power is available. It would allow for a better comparison between the QG and the GOM experiments. In particular, it would be interesting to analyze whether an exponential decrease of β * with the number of dynamic members also occurs in the GOM. This could provide a useful indication of the expected error when setting up an operational forecasting system.    Fig. 4. Quasi-geostrophic summary plot of optimal blending coefficient β * vs. the dynamic ensemble size, for the EnKF-OI and the ESRF-OI. The green solid line is an exponential fit to the ESRF-OI experiment.    7. Overlay of model ensemble fronts with the non-assimilated ocean color data (contour in mg/m 3 ) for the 12 th of July (a), the 19 th of July (b), and the 26 th of July (c). Blue color (resp. green) indicates low (resp. high) concentration of chlorophyll, and cloud covered areas are in white. The thick black line represents the front derived from SSH altimeter data. The pink lines represent the nowcast of each ensemble member of the EnKF-OI, and the red thick line is the ensemble mean. The thick yellow line is the EnOI front. The points A and B represent the locations where correlations are analyzed in Section 4.4