Skill-assessments of statistical and Ensemble Kalman Filter data assimilative analyses using surface and deep observations in the Gulf of Mexico

Abstract

A new data assimilation algorithm (Quasi-EnKF) in ocean modeling, based on the Ensemble Kalman Filter scheme, is proposed in this paper. This algorithm assimilates not only surface measurements (sea surface height), but also deep (~2000 m) temperature observations from the Gulf of Mexico into regional ocean models. With the use of the Princeton Ocean Model (POM), integrated for approximately two years by assimilating both surface and deep observations, this new algorithm was compared to an existing assimilation algorithm (Mellor-Ezer Scheme) at different resolutions. The results show that, by comparing the observations, the new algorithm out-performs the existing one.

Keywords

data assimilation / deep observation / Gulf of Mexico

Introduction

Loop Current and associated rings dominate the upper-layer (z ≈ 0 to-1000 m) circulation in the Gulf of Mexico (Oey et al.,2005). A recent collection of papers in a book edited by Sturges and Lugo-Fernandez (2005) gives useful information on the physical oceanography in the Gulf of Mexico. For modelers and those interested in ocean predictions involving meso-scale features, studying the Loop Current and rings in the Gulf of Mexico will enhance their research due to the area’s following aspects: (1) energetic (speeds 1-2 m/s), hence the large signal-to-noise ratio; (2) relatively large scales (>200 km), hence easily resolvable by satellites and by most modern models with grid sizes of approximately 10 km; (3) entirely contained within the Gulf of Mexico; a basin of large scales (~3000 km) with dynamics that can be affected by planetary beta effects; (4) tendency to be significantly affected by interaction with slopes and shelves; and (5) offer excellent case studies of interaction between energetic ocean currents (eddies) and tropical cyclones which frequently visit the Gulf during the hurricane season. Unfortunately, the same conditions listed above also contribute to complex dynamics, making it challenging to simulate and predict the Gulf’s circulation; particularly for the behaviors of the Loop Current and rings.

Assimilating observations into models is one way to constrain the solution for deriving an initial model state which is closer to the observed, presumably the “truth,” state. However, even with a good initial state, the strong nonlinearity almost certainly guarantees that the simulated Loop Current and rings will, after some time, depart from the (future) observed state. With the use of the breeding method (Kalnay,2003), Yin and Oey (2007) introduced a random field of initial perturbations into their primitive-equation circulation model (the Princeton Ocean Model (POM); Mellor (2004)) of the Gulf of Mexico, and showed that a fully-evolved bred vector, consisting of cyclones and anticyclones, could be developed within approximately six to eight weeks. This was also the time scale in which their ensemble model forecasts tended to significantly depart from the satellite altimeter observations of the positions of the Loop Current and rings. Similar conclusions were also obtained by Oey et al. (2005). The authors focused on predicting the detailed frontal evolutions of the Loop Current and rings, and showed that one large source of error came from the inexact knowledge of the positions of these fronts. It is clear that good analyses from data assimilation are necessary for good forecasts.

Whether or not an analysis is good may be assessed by various means. Kantha et al. (2005) presented a comprehensive comparison of their forecast frontal positions (of the Loop Current and rings) with the observed fronts determined from satellite sea-surface temperature (SST), altimeter data (sea-surface height (SSH)), and surface drifters. Comparisons of Acoustic Doppler Current Profiler (ADCP) currents with the model were also presented. The authors also demonstrated that their relatively simple optimal interpolation (OI) data assimilation scheme, using satellite SSH-anomaly (SSHA), can yield good results. Chassignet et al. (2005) used ocean color data from the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) to assess the skills of five data-assimilative models. Various assimilation schemes were used including OI, statistical, nudging, and the Cooper and Haines’ (1996) technique. Data used included satellite SSHA, and Modular Ocean Data Assimilation System (MODAS) temperature and salinity analyses. As pointed out by the authors, ocean color can characterize the performance of each system, but does not allow for a quantitative assessment of the systems. Nonetheless, the authors showed that the modeled positions of the Loop Current and rings are in good qualitative agreement with the low-chlorophyll regions representative of the Caribbean water. Oey et al. (2005) also used a relatively simple OI assimilation scheme in which the Mellor and Ezer (1991) surface-subsurface correlation functions were used to project satellite SSHA into the density field. The forecast skills were assessed by comparing the frontal positions of the Loop Current and rings with those derived from satellite SSH, SST, and drifters (c.f. Kantha et al.,2005). The authors showed that the inaccurate knowledge of the initial frontal position is the largest source of error, yet the model has the capability to show that after about two weeks, the forecast error is less than the persistence (keeping the front at the initial position for all time) error.

With the exception of the frontal-position comparisons (Kantha et al.,2005; Oey et al.,2005), the above model-observation comparisons are largely qualitative. Moreover, none of these works quantitatively evaluates the modeled currents against in situ observations. Lin et al. (2007) used Oey et al.’s (2005) method and extended it to also assimilate surface drifters. An extensive skill assessment was conducted. In addition to standard comparisons against satellite data, the authors quantitatively evaluated the modeled surface currents against drifter-derived currents for simulations with and without drifter assimilations, and also with and without satellite SSHA assimilations. Quantitative comparisons were also made for the near-surface currents and temperature profiles against ship-board data (ADCP and temperature measurements). In addition, Lin et al. (2007) quantitatively compared the various modeled profiles in the upper layer (0 to -800 m) against ADCP data at one station over the Sigsbee Escarpment (water depth ≈ 2000 m). The first two modeled empirical orthogonal functions (EOFs) not only compared well with the observed, and complex vector correlations but also indicated good model skills.

In this work, we will follow Lin et al. (2007) to assess the skills of two data assimilation schemes applied to the POM. The assimilation schemes include those of Mellor and Ezer (1991) as well as the Ensemble Kalman OI (Evensen,2003) methods. Tests with various horizontal and vertical model grid resolutions are also conducted. A more recent (2003-2004) comprehensive set of observational data obtained by the Minerals Management Service (MMS) (Donohue et al. 2006) from the east-central Gulf of Mexico (Fig. 1) is used to evaluate the models. This paper will focus primarily on both the application of the new data assimilation schemes into the model and the evaluation of the model states by observations. Deep observations are not only used in assimilation (temperature), but also in the evaluati0n of the model (velocity).

Section 2 describes the model and Section 3 the data-assimilation schemes. Section 4 evaluates the modeled velocity fields against the observations. Results are shown in Section 5, and the paper concludes with a discussion and summary in Section 6.

The model

Our circulation model for the Caribbean Sea and Gulf of Mexico is based on the Princeton Ocean Model (Mellor,2004). This model has been extensively tested for process in addition to realistic simulations, especially in the Gulf of Mexico (Oey et al.,2005; Lin et al.,2007; Yin and Oey,2007). A brief description is given below.

The model domain includes the northwestern Atlantic Ocean west of 55°W, as shown in the inset of Fig. 1. At 55°W, estimates of inflow and outflow transports are specified in combination with radiation conditions. The baroclinic velocities are specified using the radiation conditions. Climatological temperature and salinity are specified during inflow and advected out using one-sided differencing at outflow. Details of open boundary conditions are in Oey and Chen (1992a,1992b). When data assimilation is used in the interior of the model domain, tests show that different open-boundary specifications (e.g., perturbations from climatology) have little impact on the interior solutions.

According to the method given in Oey (1995,1996a,1996b), the model is forced by both wind and by monthly discharges from 34 rivers along the northern Gulf coast. Its horizontal grid-size is variable (orthogonal curvilinear grid); it is approximately 10 km in the Loop Current and northwestern Caribbean Sea, and about 5 km in the northeastern Gulf of Mexico. There are 25 terrain-following (the so-called sigma-coordinate) layers, ten of which are located in the top 250 m for local water depth ≈ 2500 m. The Mellor and Yamada’s (1982) turbulence closure scheme modified by Craig and Banner (1994) to affect wave-enhanced turbulence near the surface is used. To account for mixing in stable stratification (e.g., internal waves; MacKinnon and Gregg,2003), Mellor’s (2001) modification of a Richardson-number-dependent dissipation is used. The horizontal viscosity is modeled using the Smagorinsky’s shear-dependent formula with a constant= 0.1, and the corresponding diffusivity set ten times smaller. With nesting, the resulting (doubled-resolution) grid size near the Sigsbee escarpment is approximately 3 km. A fourth-order scheme is used to evaluate the pressure-gradient terms (Berntsen and Oey,2010) and, in combination with high resolution and subtracting of the mean r-profile, results in a very small truncation error of O (mm/s) (c.f. Oey et al.,2003).

Data assimilation

We assimilate satellite SSHA (δη_o, from www.aviso.oceanobs.com) and SST (from the United States GODAE, www.usgodae.org) using two schemes (below) to derive analyses. In addition to the Mellor-Ezer scheme, a new data assimilation algorithm (called Quasi-EnKF) is introduced based on the Ensemble Kalman Filter scheme.

The Mellor-Ezer scheme

In this scheme, satellite data are assimilated into the model according to the methodology given in Mellor and Ezer (1991) and Ezer and Mellor (1994). In this method, SSHA is projected into the subsurface temperature field using pre-computed correlation factors derived from a long-time (≈ 10 years) prognostic integration that has yielded a statistical equilibrium eddy field. Thus the resulting temperature anomaly (δT) is (operator<>is time-averaging, and T is the potential temperature):

(1)\delta T(x,y,z,t)=FT(x,y,z)\delta {\eta }_o(x,y,t),

where the correlation factor is (δη = model SSHA)

(2a)F_T=\langle \delta T\delta \eta \rangle /\langle \delta {\eta }^2\rangle ,

and the corresponding correlation coefficient is

(2b)C_T=\langle \delta T\delta \eta \rangle /{(\langle \delta T^{\mathrm{2}}\rangle \langle \delta {\eta }^{\mathrm{2}}\rangle )}^{\mathrm{1}/\mathrm{2}}.

Ezer and Mellor (1994) assimilate along-track δη_o data assuming a linear-saturation error growth model for the first-guess error. Our experience has been that if AVISO δη_o maps are assimilated, the following simplified formula (due originally to Ezer et al. (2003), private communication; see Wang et al. 2003) suffices:

(3)T^a=T+(\mathrm{2}{R}_A{C}_T^{\mathrm{2}}/(\mathrm{1}+\mathrm{2}{R}_A{C}_T^{\mathrm{2}}-C_T^{\mathrm{2}}))(T_O-T),

where T is the model (first-guess) temperature, T^a denotes the analysis temperature, R_A is the ratio of the assimilated time step Δt_A to the de-correlation time scale Δt_E of the model eddy field, and T_O is the “observed” temperature inferred from Eq. (1),

(4)T_O=\langle T\rangle +F_T\delta {\eta }_o.

Instead of using the model mean for<T>in (4), our past experience has been that setting<T>= T_C, the observed temperature climatology, helps to control long-term (~10 years) drift in the model. For the present application, the differences are small. Eq. (3) assumes that the AVISO map errors are small compared to the model errors, and that Δt_A<<Δt_E. We follow Ezer and Mellor (1994) and set Δt_A = 1 day. The Δt_E is estimated from the above-mentioned 10-year prognostic model run and is ≈ 30 days in regions of the Gulf of Mexico dominated by the Loop Current and rings. This may be compared with the value of 20 days used by Ezer and Mellor (1994) for the Gulf Stream, which therefore appears to have shorter meander and eddy evolution time scales. The Δt_E is also proportional to the time scale of the model error growth, and the 30-day value is consistent with the Oey et al (2005; see also Yin and Oey,2007) findings of predictability time scales of approximately one month for the Loop Current and its associated rings. As pointed out by Ezer and Mellor (1994), the assimilation of Eq. (3) is such that T^a ≈ T_O in regions where the correlation is high (C_T² ≈ 1), but T^a ≈ T where the correlation is low. A similar assimilation of SST is also carried out after Eq. (3) with C_T and F_T replaced by the corresponding functions that use δ(SST) in place of δη in Eq. (2). The SSHA and SST assimilations complement each other: SSHA assimilation is most effective over deep waters (for isobath>500 m), while SST assimilation influences waters on shallow shelves. For more details see Wang et al. (2003), Oey et al. (2005), Lin et al. (2007), and Yin and Oey (2007).

The Ensemble Kalman Filter scheme

A sequential data assimilation method is used:

(5a)x^a=x^f+K(x^o-H{x}^f),

(5b)K=P^f{H}^T{(H{P}^f{H}^T+R)}^{-\mathrm{1}},

where x^f is forecast state, x^o is observation, x^a is analysis state, H is observational operator, K is gain matrix, P^f is forecast error covariance, R is observation error covariance, ( )^T denotes matrix transpose, and ( )^-1 is matrix inverse. We assume that R is known. The ensemble optimal interpolation (EnOI) method is used for assimilation; however, the error covariance P^f is not fixed (below). Instead of Eq. (5b), we use the following formula to calculate K:

(6)K=\alpha P^f{H}^T{(\alpha H{P}^f{H}^T+R)}^{-\mathrm{1}},

where \alpha \in (\mathrm{0},\mathrm{1}) is a parameter that controls the weight between P^f and R.

Estimating the forecast error covariance

The forecast error covariance (P^f) is estimated by an ensemble method. In the full Ensemble Kalman Filter (EnKF) method, the P^f is time-dependent and calculated at every assimilation time-step Dt_a (= 1 day, say). One way is to compute, at every Dt_a, the P^f from the solutions of a set of perturbed experiments. In EnOI, the P^f is fixed; one way is to take the temporal mean of error covariances from a multi-year integration, and then use that mean as the time-independent P^f in Eqs. (5) and (6) for the analyses covering the same multi-year period. The EnKF is very cpu-intensive; the EnOI is simple, but it does not take into account the “error of the day.” In this work, we have designed a slightly different algorithm that partially takes into account the evolution of the background field. Instead of computing the P^f at every Dt_a as done in EnKF, the ensemble of perturbed experiments are conducted at every t_E, where t_E corresponds to the time scale in which the eddy field has evolved “significantly.” Yin and Oey’s (2007) bred-vector analysis shows that random perturbations would evolve into a well-organized eddy field in approximately eight weeks. This time scale agrees with past experiences of modeling the Loop Current and rings (e.g., Oey et al.,2003); we take, then, t_E = 60 days.

The ensemble is then comprised of a set of model solutions corresponding to the analysis time interval Dt_a. The initial states of different ensemble members are given small random perturbations (see Yin and Oey,2007 for details). These ensemble calculations are done every t_E = 60 days. There are N_ens= 30 ensemble solutions denoted as x_i(i=\mathrm{1},\mathrm{2},\cdots ,N_{ens}), where x_i is a column vector including all model variables (η, T, u, v, S, etc.). The forecast error covariance P^f is then estimated by:

(7)P^f={\displaystyle \frac{(X-\mu ){(X-\mu )}^T}{N_{ens}-\mathrm{1}}},

where matrix X=[x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_{N_{ens}}], vector \mu ={\displaystyle \frac{1}{N_{ens}}}{\displaystyle \sum _{i=1}^{N_{ens}}{x}_i}.

P^f contains the covariance (or correlation) between any two model variables, as well as that between two model grids of one model variable.

Theoretically, the estimated P^f will approach the true P^f as N_ens approaches infinity. However, this is not practical for a very large N_ens due to computational time, so there is likely some spurious covariance in P^f at large spatial distances. The following piecewise continuous function (Gaspari and Cohn, 1999) is adopted to filter out that kind of spurious covariance.

(8)C_{\mathrm{0}}\left(z,{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},c\right)=\{\begin{array}{ll}\mathrm{0},\hfill & \mathrm{2}c\le \left|z\right|\hfill \\ -{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{5}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{4}}+{\displaystyle \frac{\mathrm{5}}{\mathrm{8}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{3}}-{\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{2}}+\mathrm{1},\hfill & \mathrm{0}\le \left|z\right|\le c\hfill \\ {\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{5}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{4}}+{\displaystyle \frac{\mathrm{5}}{\mathrm{8}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{3}}+{\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{\mathrm{2}}-\mathrm{5}\left({\displaystyle \frac{\left|z\right|}{c}}\right)+\mathrm{4}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\left({\displaystyle \frac{\left|z\right|}{c}}\right)}^{-\mathrm{1}},\hfill & c\le \left|z\right|\le \mathrm{2}c\hfill \end{array}

where c is a given positive constant to control the correlation length, and z is the input distance.

In summary, our assimilation method is a little more sophisticated than the standard EnOI method. It adopts the EnKF idea in that the P^f changes with time, although discretely at 60-day intervals. The P^f is calculated by ensemble experiments with 30 ensemble members each of which is initially given random perturbations. An example is shown in Fig. 2 of the time-dependent correlation coefficient between η (SSHA) and T (temperature), i.e.,CC(\eta ,T)={\displaystyle \frac{P^f(\eta ,T)}{{\sigma }_{\eta }{\sigma }_T}}, where σ denotes the standard deviation. We will refer to our scheme as the QEnKF – Quasi-EnKF.

Assimilating altimetry SSHA and moored temperature data

We first describe how a two-dimensional (x,y) dataset, such as the altimetry SSHA, is assimilated. To assimilate at the model grid {(x_{i,j},y_{i,j})}^m, the following algorithm (called QEnKF SSHA) is used:

(1) Given a small (x,y) region (e.g., a circle of radius<1.0°; Fig. 3(a)) containing the model grid, find all the SSHA observational grids {(x,y)}_l^{\eta } inside this region;

(2) For each {(x,y)}_l^{\eta }, find the smallest quadrilateral that encloses it, formed by four model grids {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}})}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}})}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}})}^m, and {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}})}^m; one of which is {(x_{i,j},y_{i,j})}^m, the point to be assimilated (assumed to be {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}})}^m in the figure);

(3) Use {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}})}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}})}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}})}^m, {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}})}^m and {(x,y)}_l^{\eta } to define the observational operator H;

(4) Compute HP^fH^T, the covariance of SSHA at observational grids {(x,y)}_l^{\eta }. Compute P^fH^T = the covariance between temperature at model grids {(x_{i,j},y_{i,j},z_k)}^m and the SSHA at observational grids {(x,y)}_l^{\eta }.

(5) Compute gain matrix K using HP^fH^T, P^fH^T and R (assumed known) by Eq. (6);

(6) Compute analysis temperature at {(x_{i,j},y_{i,j},z_k)}^m by Eq. (5a).

To assimilate the moored temperature observations (three-dimensional data), the following algorithm (called QEnKF mooring) is used (see Fig. 3(b)):

(1) Given a small (3-D) region (e.g., a ball) containing the model grid, find all the moored temperature observations {(x,y,z)}_l^T inside this region;

(2) For each {(x,y,z)}_l^T, find the smallest quadrilateral with four model grids in horizontal direction: {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}})}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}})}^m,

{(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}})}^m

and {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}})}^m. In addition, find two nearest model vertical levels z_d and z_u, such that . Then come up with a 3-D grid box with eight model grids: {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}},z_d)}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}},z_d)}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}},z_d)}^m, {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}},z_d)}^m, {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}},z_u)}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}},z_u)}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}},z_u)}^m and {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}},z_u)}^m.

(3) Use {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}},z_d)}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}},z_d)}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}},z_d)}^m, {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}},z_d)}^m, {(x_{i\mathrm{1},j\mathrm{1}},y_{i\mathrm{1},j\mathrm{1}},z_u)}^m, {(x_{i\mathrm{1},j\mathrm{2}},y_{i\mathrm{1},j\mathrm{2}},z_u)}^m, {(x_{i\mathrm{2},j\mathrm{1}},y_{i\mathrm{2},j\mathrm{1}},z_u)}^m, {(x_{i\mathrm{2},j\mathrm{2}},y_{i\mathrm{2},j\mathrm{2}},z_u)}^m and {(x,y,z)}_l^T to model the observational operator H.

(4) Compute HP^fH^T, the covariance of temperature at observational grids {(x,y,z)}_l^T , and then compute P^fH^T, the covariance of temperature at between model grids {(x_{i,j},y_{i,j},z_k)}^m and observational grids {(x,y,z)}_l^T.

(5) Compute gain matrix K using HP^fH^T, P^fH^T and R (assumed known) by Eq. (6).

(6) Compute analysis temperature at {(x_{i,j},y_{i,j},z_k)}^m by Eq. (5a).

Experiments

We carried out three assimilation experiments by using an identical initial state on Nov. 04, 2002. All experiments ended on Nov. 24, 2004. The first experiment was a QEnKF SSHA assimilation, the second a QEnKF SSHA+ QEnKF mooring assimilation, and the third a QEnKF SSHA assimilation with double resolution in the vertical direction. SSHA observations are available every day during the whole model integrations, while mooring observations (including deep observations) are only available daily from March 2003 to April 2004. Thus the assimilation time interval is one day for all assimilation experiments.

To compare these three experiments, we added four experiments using other algorithms: (ⅰ) the Mellor-Ezer scheme; (ⅱ) the Mellor-Ezer scheme with double resolution in the vertical direction; (ⅲ) the Mellor-Ezer scheme with half resolution in the horizontal direction; and (ⅳ) a model run without assimilation. We then have the following seven cases to compare (Table 1). (Please note that the Mellor-Ezer scheme only uses SSHA observations.)

Results

Compare (u,v) time series between model and observations at four mooring locations

A complex correlation coefficient (CC) is computed between model solutions and observations for (u,v) time series. There are three values corresponding to a complex CC: magnitude, angle, and real part in Tables 2-5. Results show that none of the cases produce good estimations at deep levels. However, case C is the best near the surface.

Comparison of mean spatial correlation

In comparison with the no assimilation experiment (case G), all assimilation analysis significantly improved SSHA solutions (mean correlation coefficient>0.8 in Table 6)

Comparisons of temperature time series

To compare the temperature between model and mooring observations, the mean of the absolute value of their difference was calculated over all observational time and at all observational grids of one mooring. The values are shown in Table 7 as “Mean(abs)”. In addition, the standard deviation of their difference was shown as “Std(ori)” in the table. These two statistics show the model accuracy in terms of observations.

Comparison of temperature at depths

To validate the model solutions at different locations (depths of moorings), the correlation coefficients of two time series (model and observation) were calculated and shown in Table 8.

Discussion

In this paper, we discuss the data-assimilative analyses carried out on POM in the Gulf of Mexico using two schemes. The goal is to quantify the accuracy of currents obtained from a model that assimilates altimetry SSHA data and mooring temperature data.

As seen in Table 8, all assimilations are quite good above the 1000 m level, with the best being the QEnKF SSHA+ QEnKF mooring (case C). This is reasonable given it is the only case that assimilates mooring temperature observations into the model. Below the 1000 m level, the correlation coefficients for all cases are less than 0.3, some of them even at negative numbers. There are two possible reasons for the small correlation coefficients: (π) the analysis temperature is far away from the mooring observations, and (θ) the observations are not accurate in deep water. As seen in Fig. 5, the L4 mooring observations at 2000 m, 2500 m, and 3250 m are almost constant during the whole observational period; thus, we could suspect those observations are not accurate.

From the mean absolute value of T (model-observation) from Table 7, the best is still the QEnKF SSHA+QEnKF mooring case (C), which proves assimilation with temperature observation has the greatest impact on model temperature. The second best is one of three cases: A, B, and E. Surprisingly, the analysis temperature of case F is as good as case A.

Inside the Gulf of Mexico, all assimilations are quite good with the exception of four locations: (86°W, 27.5°N), (93.5°W, 22.5°N), (95.5°W, 24°N), (94°W, 20.5°N). As stated earlier, the best analysis belongs to QEnKF SSHA (Fig. 4).

In addition, from the correlation coefficient between SSHA and temperature (Fig. 2), most of the positive values occur above 500 m. Below that, most are negative, which means increasing SSHA will decrease temperature at the deep water level in QEnKF scheme.

After comparing seven experiments, including the new algorithm QEnKF, and taking into consideration that no experiment shows good estimations of (u,v) in deep water, the new algorithm performs best in terms of observation comparisons of SSHA and temperature.

References

[1] Berntsen J, Oey L Y (2010). Estimation of the internal pressure gradient in σ-coordinate ocean models: Comparison of second-, fourth-, and sixth-order schemes. Ocean Dyn, 60(2): 317-330

[2] Chassignet E P, Hurlburt H E, Smedstad O M, Barron C N, Ko D S, Rhodes R C, Shriver J F., Wallcraft A J, Arnone R A (2005). Assessment of data assimilative ocean models in the Gulf of Mexico using ocean color. Geophysical Monograph-American Geophysical Union, 161: 87

[3] Cooper M, Haines K (1996). Altimetric assimilation with water property conservation. J Geophys Res, 101(C1): 1059-1077

[4] Craig P D, Banner M L (1994). Modeling wave-enhanced turbulence in the ocean surface layer. J Phys Oceanogr, 24(12): 2546-2559

[5] Donohue K, Hamilton P, Leaman K, Leben R, Prater M, Watts D R, Waddell E (2006). Exploratory Study of Deepwater Currents in the Gulf of Mexico, Vol. I and II. US Dept. of the Interior, Minerals Management Service, Gulf of Mexico OCS Region, New Orleans, LA. OCS Study MMS, 2006-073

[6] Evensen G (2003). The Ensemble Kalman Filter: theoretical formulation and practical implementation. Ocean Dyn, 53(4): 343-367

[7] Ezer T, Arango H, Hermann A (2003). Terrain-Following Ocean Models Users Workshop, Seattle, WA, <month>Aug</month>. <day>4-6</day> AOS Program, Princeton University, 19

[8] Ezer T, Mellor G L (1994). Continuous assimilation of Geosat altimeter data into a three-dimensional primitive equation Gulf Stream model. J Phys Oceanogr, 24(4): 832-847

[9] Gaspari G, Cohn S E (1999). Construction of correlation function in two and three dimensions. Q J R Meteorol Soc, 125(554): 723-757

[10] Kalnay E (2003). Atmospheric Modeling, Data Assimilation and Predictability. New York: Cambridge University Press, 341

[11] Kantha L, Choi J K, Schaudt K J, Cooper C K (2005). A Regional Data-Assimilative Model for Operational Use in the Gulf of Mexico. Geophysical Monograph-American Geophysical Union, 161, 165

[12] Lin X H, Oey L Y, Wang D P (2007). Altimetry and drifter data assimilations of Loop Current and eddies. J Geophys Res, 112(C5): C05046

[13] MacKinnon J A, Gregg M C (2003). Shear and baroclinic energy flux on the summer New England Shelf. J Phys Oceanogr, 33(7): 1462-1475

[14] Mellor G L (2001). One-dimensional ocean surface layer modeling, a problem and a solution. J Phys Oceanogr, 31(3): 790-809

[15] Mellor G L (2004). User’s Guide for a Three-dimensional, Primitive Equation, Numerical Ocean Model. Program in Atmospheric and Oceanic Sciences, Princeton University, 42

[16] Mellor G L, Ezer T (1991). A Gulf Stream model and an altimetry assimilation scheme. J Geophys Res, 96(C5): 8779-8795

[17] Mellor G L, Yamada T (1982). Development of a turbulence closure model for geophysical fluid problems. Rev Geophys, 20(4): 851-875

[18] Oey L Y (1995). Eddy- and wind-forced shelf circulation. J Geophys Res, 100(C5): 8621-8637

[19] Oey L Y (1996 a). Flow around a coastal bend: a model of the Santa Barbara Channel eddy. J Geophys Res, 101(C7): 16,667-16,682

[20] Oey L Y (1996 b). Simulation of mesoscale variability in the Gulf of Mexico: sensitivity studies, comparison with observations, and trapped wave propagation. J Phys Oceanogr, 26(2): 145-175

[21] Oey L Y, Chen P (1992a). A model simulation of circulation in the Northeast Atlantic shelves and seas. J Geophys Res, 97(C12): 20,087-20,115

[22] Oey L Y, Chen P (1992b). A nested-grid ocean model: with application to the simulation of meanders and eddies in the Norwegian Coastal Current. J Geophys Res, 97(C12): 20,063-20,086

[23] Oey L Y, Lee H C, Schmitz W J (2003). Effects of winds and Caribbean eddies on the frequency of Loop Current eddy shedding: a numerical model study. Journal of Geophysical Research: Oceans (1978-2012), 108(C10)

[24] Oey L, Ezer T, Lee H (2005). Loop Current, rings and related circulation in the Gulf of Mexico: a review of numerical models and future challenges. Geophysical Monograph-American Geophysical Union, 161, 31

[25] Sturges W, Lugo-Fernandez A (2005). Circulation in the Gulf of Mexico: Observations and Models, Geophys. Monogr. Ser., Vol. 161. Washington, D C: AGU, 347

[26] Wang D P, Oey L Y, Ezer T, Hamilton P (2003). Near-surface currents in DeSoto Canyon (1997-99): comparison of current meters, satellite observation and model simulation. J Phys Oceanogr, 33(1): 313-326

[27] Yin X Q, Oey L Y (2007). Bred-ensemble ocean forecast of Loop Current and rings. Ocean Model, 17(4): 300-326

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg