1 Introduction

A three-dimensional variational (3DVAR) assimilation scheme for Lagrangian trajectories was developed and tested by Taillandier et al. (2006b) and Taillandier et al. (2010) using Argo float surfacing positions during a 3-month period in the North-Western Mediterranean Sea. In these studies, encouraging results were obtained although only four Argo floats were used in the numerical experiments. Here, a numerical study, taking into account all Argo float surfacing positions available in the Mediterranean Sea, has been undertaken for a 3-year period (2005–2007). The model results, with and without Argo float trajectory assimilation, were assessed with special focus on the intermediate velocity fields at 350 m depth.

The numerical experiments have been performed using the MFS, which has been in operational use for a decade and has continuously been providing forecasts and analyses for the region (Pinardi et al. 2003; Tonani et al. 2008). These forecasts yield local- and basin-scale information of the state of the sea, e.g., sea level, velocity, temperature, and salinity fields, and can be useful for search-and-rescue missions and for tracking oil spills (Coppini et al. 2010). For these purposes, it is crucial to provide state-of-the-art model output, and it has been established that OceanVar (Dobricic and Pinardi 2008) is capable of improving significantly the overall quality of the MFS model fields (Dobricic et al. 2007). At present, the assimilated observational datasets are obtained from both remote sensing and in situ measurements, such as satellite-observed sea-level anomalies (SLA, Le Traon et al. 2003), temperature profiles from expendable probes (XBT, Manzella et al. 2007), as well as temperature and salinity profiles from Argo floats (Poulain et al. 2007). In the present study, the Argo float positions obtained at the sea surface (Poulain et al. 2007; Menna and Poulain 2010) will be added to this ensemble, with the primary aim of constraining and improving the modeled intermediate velocity fields in the Mediterranean Sea.

The Mediterranean Sea is a deep semi-enclosed basin connected to both the Atlantic Ocean and the Black Sea through narrow straits. Its bathymetry is characterized by two major sub-basins, the Western and Eastern Mediterranean, which are separated by a shallow sill (∼400 m) in the Sicily Channel, cf. Fig. 1. The thermohaline circulation is mainly driven by buoyancy loss due to the evaporation as well as the low precipitation and fresh-water inflow in the area, thus, making the Mediterranean Sea one of the largest concentration basin in the world. This inverse estuarine circulation can be schematically described as a balance between the inflow from the Atlantic and the Levantine Intermediate Water (LIW) outflow (Benzohra and Millot 1970; Millot 1999). Most of the dense water is formed in the Gulf of Lions and the Levantine Basin subduction zones, in particular, the LIW is formed in the North-Eastern part of the Levantine Basin, where, due to buoyancy loss it sinks to its characteristic depth of ∼200–500 m depth. Thereafter, it spreads across the Mediterranean through various pathways and finally exits into the Atlantic Ocean through the Straits of Gibraltar, hereby closing the loop of what is known as the Mediterranean conveyor belt, cf. Pinardi and Masetti (2000).

Fig. 1
figure 1

Observed Argo float positions in the Mediterranean Sea in 2005–2007

Full size image

The Mediterranean circulation is, furthermore, characterized by well-defined coastal currents, such as the Algerian Current (AC), as well as small-, meso-, and sub-basin-scale eddy structures in the interior of the basin, cf. Millot (1991), Robinson et al. (1991). Modeling the meandering of these currents and calculating the evolution of the eddies is highly complex, thus the velocity-field forecasts are often associated with large uncertainties due to the inherent nonlinear dynamics of the circulation (cf. Molcard et al. 2002). In this context, assimilation of Lagrangian trajectories in operational ocean models could offer a possibility to reduce some of these errors and to ultimately yield more accurate ocean forecasts.

The contents of the manuscript are disposed as follows: Section 2 provides a general description of the Mediterranean Forecasting System and the OceanVar assimilation scheme, Section 3 describes the numerical experiments, Section 4 subsequently presents and discusses the results, and finally Section 5 offers some conclusions.

2 The Mediterranean Forecasting System

The MFS consists of three fundamental constituents: the data collection network, the ocean general circulation model (OGCM), and the OceanVar assimilation scheme. The MFS daily cycle and its coupling is schematically illustrated in Fig. 2. Next, the Argo trajectory dataset will be presented, whereafter the complete forecasting system is described.

Fig. 2
figure 2

Illustration of the MFS daily cycle and the coupling of the forcing (F), initial conditions (IC), the OGCM, OceanVar, and the observations (Obs). The OGCM produces 1-day forecasts (background fields) and OceanVar calculates model field corrections (analyses) daily

Full size image

2.1 Argo float trajectories

Data from both Apex (manufactured by Webb Research Corporation, USA) and Provor (produced by NKE Electronics, France) profiling floats were provided by the Coriolis Operational Oceanography data center and extensively quality checked by The National Institute of Oceanography and Experimental Geophysics (OGS) in Trieste. These floats, commonly referred to as MedArgo floats, started to be deployed in the Mediterranean Sea in 2003 within the framework of the MFSTEP project. The MedArgo floats are programmed to perform continuous cycles, in which the float descend from the sea surface to 350-m parking depth, where it drifts for a 4.5-day period. The cycle is completed by a 700-m dive (2,000 m every ten cycles), whereafter the float re-emerges to the sea surface and makes contact with the Argos satellite system and transmits the data (e.g., surfacing coordinates, temperature, and salinity profiles). When the data-transfer procedure has been completed, after approximately 6 h, the float begins the next cycle by descending back to the parking depth (cf. Menna and Poulain 2010). Each Argo float file holds information of time, float surfacing positions as well as the corresponding water depth at these locations. The depths were retrieved from the Smith and Sandwell Bathymetry (Smith and Sandwell 1997) which is based on Satellite Altimeter data.

When Argo float coordinates are available, the OceanVar trajectory model computes 5-day trajectory forecasts. Moreover, an algorithm in the daily assimilation cycle searches the observational time series for pertinent data on the “present day” (t = t f f=final) and on the “preceding Argo cycle day” (t i  = t f -Δt; i=initial, Δt = 5 days), in order to perform trajectory data assimilation. In conjunction with this, the depths on these two occasions were examined to preclude erroneous trajectory estimates when the float might have been stuck to the bottom. Problems of this type were avoided by constraining the data selection so that only coordinates obtained at locations with depths greater than or equal to 400 m were accepted (approximately 4% of the data was rejected). The coverage of the Argo float positions during 2005–2007 is displayed schematically in Fig. 1, and the OceanVar trajectory model is presented further in Section 2.4.

2.2 The ocean general circulation model

The Océan Parallélise code (Madec et al. 1998) was adapted to the Mediterranean Sea by Tonani et al. (2008) and served as OGCM in the numerical experiments. This model is based on the primitive equations subjected to the Boussinesq and incompressibility approximations, and thereafter, discretized on a spherical grid with a horizontal resolution of 1/16° × 1/16° (∼5–7 km depending on latitude). The model depths are described by 72 unevenly spaced levels, with a 3-m layer thickness near the surface and 300 m near the bottom, hereby resolving the Mediterranean bathymetry in a realistic manner. The OGCM was spun up from a state of rest and forced by 1/2° horizontal resolution, 6-h atmospheric fields from the European Centre for Medium-Range Weather Forecasts.

The water exchanges through the Gibraltar and Dardanelles straits are dealt with in different manners, the latter being implemented as a river (surface boundary condition) with monthly inflow and salinity climatology parameterizations (Adani et al. 2011; Kourafalou and Barbopoulos 2003). The open boundary in the Atlantic, on the other hand, is described by box (outer limit at 18° W) in which the flow across Gibraltar Strait is relaxed to climatology and vanishing currents are prescribed at the Atlantic model boundaries. The temperature and salinity along these borders are relaxed to Atlantic climatology at all depths. Moreover, the heat flux in the Mediterranean Sea is corrected by relaxation of the modeled surface-layer temperatures (Dobricic et al. 2005) towards sea-surface temperature observations (SST, Marullo et al. 2007).

Due to the fact that Mediterranean basin is a concentration basin, i.e., it has a net water loss since the evaporation exceeds the precipitation, the water balance is defined to conserve mass under the condition that the net water flux at the sea surface is negative. The changes in the sea-surface salinity due to evaporation and precipitation is, thus, modeled by relaxation to the monthly MEDATLAS climatology values (MEDAR/MEDATLAS Group 2002), cf. Tonani et al. (2008).

2.3 The OceanVar assimilation scheme

The development of variational assimilation schemes for atmospheric forecasting models started in the late 1980s (Lorenc 1986), and research efforts since have resulted in the highly advanced four-dimensional variational schemes. The variational assimilation techniques for operational ocean forecasting purposes are not, at present, as sophisticated as the methods applied in operational weather forecasting; however, continuous progress is being made in this field (cf. Dobricic 2009). Implementation of 3DVAR assimilation schemes in ocean models is not straight forward, due to the variable model-domain lateral boundaries (coast lines) and bathymetry and to overcome these challenges new practical solutions have had to be found (cf. Dobricic and Pinardi 2008).

Although the amount of oceanic observations is quite modest compared with those available for the atmosphere, the assimilation of these somewhat limited datasets was found to yield more accurate ocean forecasts (Dobricic et al. 2007). In this context, it should be mentioned that the computational cost related to the 3DVAR assimilation procedures is less dependent on the number of observations, but mostly on the size of the model state vector, (i.e., the number of model output variables). This implies that the CPU time will not be significantly affected if, at a later stage, further observational datasets were to be added to the forecasting system.

In this numerical study, OceanVar has been extended to also assimilate Argo float positions using a trajectory model as the observational operator (Taillandier and Griffa 2006; Taillandier et al. 2010). Hereby, Argo float trajectories are being added to the wide range of variables that are already being assimilated operationally, e.g., SLA and in situ observations of temperature and salinity profiles from Argo floats and temperature profiles from XBTs.

The OceanVar minimizes by iterations a cost function formulated as:

$$ \begin{array}{rll} J &=& \frac{1}{2}\left({\bf x}-{\bf x}_{b}\right)^T{\bf B}^{-1}\left({\bf x}-{\bf x}_{b}\right)\\ &&+ \frac{1}{2}(H({\bf x}) - {\bf y})^T{\bf R}^{-1}(H({\bf x}) - {\bf y}) . \label{eq:cost} \end{array} $$
(1)

Here, x is the model state vector, x b the background state vector, y the observational vector, B the background error covariance matrix, R the observational-error covariance matrix, H the nonlinear observational operator, and T denotes the vector transpose (e.g. Lorenc 1997). The model state vector contains the temperature, salinity, velocity, and sea-level model output in matrix form as x = [T S U V η]T.

In order to set-up the minimization routine, a quadratic cost function is created by linearizing Eq. 1 around the background state vector x b , yielding:

$$ J = \frac{1}{2}\delta{\bf x}^T{\bf B}^{-1}\delta{\bf x} + \frac{1}{2}({\bf H}(\delta{\bf x})-{\bf d})^T{\bf R}^{-1}({\bf H}(\delta{\bf x})-{\bf d}), \label{eq:cost_lin} $$
(2)

where δ x = x − x b are the increments and H is the linearized observational operator. The so-called misfits (the differences between the observations and the background fields) are contained in \(\bf d = [{\bf y}-{\it H}({\bf x}_b)]\), where the nonlinear operational operator, \({\it H}\), transfers the model variables onto the observational grid, hereby allowing direct comparisons between the two datasets. J is rapidly minimized in iterations, in which it is necessary to calculate J and its gradient \(\nabla J\). The calculation of \(\nabla J\) requires the application of adjoint operators for each linear operator appearing in J.

In OceanVar, the increments are written in terms of the control vector v and the matrix V:

$$ \delta{\bf x} = {\bf Vv}, \label{eq:an} $$
(3)

where V is constructed as B = VV T. Thus, the linearized cost function J can be reformulated for the transformed space as:

$$ J = \frac{1}{2}{\bf v}^T{\bf v} + \frac{1}{2}({\bf H}{\bf Vv}-{\bf d})^T{\bf R}^{-1}({\bf H}{\bf Vv}-{\bf d}). \label{eq:cost_lin2} $$
(4)

The background error covariance matrix, B, is modeled as a sequence of linear operators (cf. Dobricic and Pinardi 2008) as:

$$ {\bf V} = {\bf V}_D{\bf V}_{\rm uv}{\bf V}_{\eta}{\bf V}_H{\bf V}_V, \label{eq:cov} $$
(5)

where V V contains multi-variate temperature (T) and salinity (S) Empirically Orthogonal Functions (EOFs) computed from model output. The EOFs include information on both the temporal and the spatial variability in the Mediterranean; the year is divided in four seasons, and the basin into 13 regions due to the large differences in seasonal and local dynamics (Demirov et al. 2003; Dobricic et al. 2005). V H applies Gaussian-distributed horizontal covariances with constant correlation radii on the vertical T and S fields, hereby yielding three-dimensional spatial covariances as V S  = V H V V . Sea-surface height (SSH) error covariances are provided by V η based on the T and S corrections and a barotropic linear model explained by Dobricic and Pinardi (2008). V uv calculates the baroclinic components of the velocity correction in geostrophic balance with the surface pressure gradient, and V D applies a divergence damping filter on the velocity field.

The corrections of the velocity fields due to trajectory assimilation enter in the V uv operator, and can thereafter influence the sea level through V η , and the full 3D mass fields through V S . The observational-error covariance matrix, R, contains information of the observational errors, and hence sets the weight (related to the data reliability and representativity) in the corrections of the model state estimate.

In conclusion, in each iteration OceanVar perturbs the coefficients that are multiplied with the vertical EOFs. By the application of linear operators, those perturbations are transformed into velocity perturbations, which thereafter are used as input for the linearized trajectory model represented by the observational operator H. Furthermore, the integration of the linear trajectory model yields perturbations of the last Argo float position HVv that can be subtracted from the misfit d. The linearized trajectory model represented by H in Eqs. 2 and 4, as well as the computations of predicted and analyzed float positions will be described next.

2.4 The OceanVar trajectory model

A trajectory model was implemented in the nonlinear observational operator, thus providing a possibility to correct the modeled velocity fields by the observed Argo surfacing coordinates. The forecasted trajectories were calculated from the Eulerian model velocity fields by 5-day integrations of the particle advection equation (Eq. 6), starting on the latest observed float positions. The advection equation states that the fluid velocity in a fix point (in space and time) is equal to the velocity of a fluid parcel that is located in that position at that time; this relation is described by the nonlinear first-order differential equation:

$$ \frac{{\rm d}{\bf r}}{{\rm d}t} = {\bf u_L}({\bf r}(t),t). \label{eq:adv} $$
(6)

where r is the float position, and u L represents the Lagrangian velocities at the float parking depth during the drift period. The Eulerian velocities u can be described in the Lagrangian framework as u L (r(t),t) =\( \mathcal{L} \)(r(t))u(t), where \( \mathcal{L} \) is the bilinear Lagrange interpolator. The time-integrated advection equation yields the fully nonlinear trajectory model H(u), to be discretized and implemented in the nonlinear observational operator H for OceanVar, and is here presented for one step of integration:

$${\bf r}(t_{\!f}) = {\bf r}(t_i) + \int_{t_i}^{t_{\!f}} \mathcal{L} ({\bf r}(t)){\bf u}(t) {\rm d}t, \label{eq:dr1} $$
(7)

where t i and t f  = t i  + Δt indicates the limits of time interval (here, Δt = 5 days). However, the nonlinearity of this equation imposes a severe analytical problem when the background velocity fields at the observational positions are to be retrieved. Hence, a tangent-linear approximation was applied to Eq. 7 as proposed by Taillandier et al. (2006a), thus yielding the linearized perturbation equation, which will provide the linearized observational operator H with the Eulerian velocity increments δ u:

$$ \delta{\bf r}(t_{\!f}) = \delta{\bf r}(t_{i}) + \int_{t_{i}}^{t_{\!f}} \left( \frac{\partial{\bf u_L}}{\partial {\bf u}} \right|_{r=r_b} \left.\delta {\bf u} + \left. \frac{\partial{\bf u_L}}{\partial {\bf r}} \right|_{u=u_b} \delta {\bf r} \right) {\rm d}t. \label{eq:dr2} $$
(8)

where the position and the Eulerian velocity increments (δ r = r − r b and δ u = u − u b ) are evaluated around the background velocity u b and background position r b . Moreover, transforming the Lagrangian velocities (u L ) to Eulerian (u), the partial derivatives in Eq. 8 can be rewritten as \(\frac{\partial{\bf u_L}}{\partial {\bf u}}=\) \( \mathcal{L} \)(r b ), and \(\frac{\partial{\bf u_L}}{\partial {\bf r}}={\bf L\cdot u}(r_b)\), where L is the derivative of \( \mathcal{L} \) around the background position r b at the time t b . The equality of Eq. 8 is assumed to be fulfilled when the higher-order (nonlinear) terms in Eq. 7 are negligible (O(δ r 2,δ u 2) ∼0). This equation was given in discretized form by Taillandier et al. (2006a) in Eq. 2.

In summary, 5-day trajectory predictions are computed for each Argo float from the model velocity fields by the nonlinear trajectory model Eq. 7, starting at their respective surfacing positions. When a float has fulfilled an “Argo cycle,” the OceanVar computes the analyzed float position and the analyzed trajectory, a procedure which requires the present and prior float coordinates. The analyzed position is obtained by minimizing the distance between the present observed float position and the “background position,” i.e., the last position of the “float” trajectory produced by a 5 days long integration of the trajectory model, in the cost function (Eq. 2) through the linear operator H(δ u). From this analyzed position, the adjoint operator thereafter recalculates the trajectories between the analyzed positions and the prior observed positions.

The OceanVar assimilates data in a daily cycle while trajectories are 5 days long. This inconsistency is neglected by assuming that the innovation is constant throughout the trajectory integration time. That is, in Eq. 8, the background velocity fields are stored during the 5 days long period with the temporal frequency of 6 h, but it is assumed that δ u (Eulerian) does not change with time. After OceanVar has finished its daily routine, the initial float position for the next Argo cycle is re-set with the observed Argo float position, \({\bf r}(t_i) ={\bf r}^{\rm obs}(t_i)\), i.e., the initial float position are held fixed in the OceanVar (δ r(t i ) = 0), and only the final positions are perturbed by Eq. 8.

Comparisons between the trajectory predictions from numerical experiments (with and without trajectory assimilation) and the observed Argo float positions allows an evaluation of the consistency of the velocity field corrections. The impact of the corrections of the velocity fields can furthermore be assessed quantitatively by calculating the distance between the end points of the trajectories produced at an arbitrary t = t i and the corresponding observed float positions at t f  =t i  + Δt. A schematic figure of the predicted and analyzed trajectories at parking depth during one Argo cycle is provided in Fig. 3.

Fig. 3
figure 3

Schematic figure of observed and modeled trajectories starting from an observed float position at t = t i . The precision of the MFS velocity fields is evaluated from the differences between the observed float (obs) and the predicted (CTRL and TRAJ fcst) as well as the analyzed (TRAJ an) float positions at t = t i  + Δt, where Δt = 5 days

Full size image

3 The numerical experiments

3.1 Experimental model setup

Two numerical experiments (cf. Table 1) were undertaken in order to evaluate the impact of the assimilation of Argo float trajectories on the model state variables. Both experiments run the OGCM daily and produce 24-h mean three-dimensional temperature, salinity, and velocity fields, as well as two-dimensional SSH fields. These fields are typically denoted model background fields (Daley 1991), and in this case they are actually 1-day forecasts. After the computation of the background fields, OceanVar assimilates the “present-day” observations and the model analysis is calculated, cf. Fig. 2.

Table 1 Design of the numerical experiments
Full size table

Both the control (CTRL) and trajectory assimilation (TRAJ) experiments assimilate on a daily basis satellite-measured SLA and SST as well as in-situ observations of temperature and salinity, moreover, in the TRAJ experiment are also Argo float positions assimilated. In the case when no observations were available, no corrections were calculated and the next daily MFS cycle was initialized.

3.2 Model-result evaluation

Comparisons of the CTRL and TRAJ daily forecasts, as well as their corresponding analyses can give an indication of where and how the model fields have been corrected by OceanVar. The differences in the CTRL and TRAJ background fields are due to the different daily initial conditions, since the assimilated observational datasets are not identical, cf. Table 1 and Fig. 2. Hence, by studying the discrepancies between their respective background fields, the impact of the previously assimilated Argo trajectories on the model fields as well as the propagation of the velocity field corrections can be evaluated.

The quality of the model fields is assessed by comparisons with observed datasets, and the differences between the background fields and the observations are known as “misfits.” In this study, the misfits have been calculated daily from the OGCM output using the “present-day” SLA, temperature, salinity and Argo float positions, and since the daily assimilation procedure has not yet been performed this comparison can be regarded as independent. In the next step, the observations are assimilated by OceanVar, whereafter the analysis is compared with the same observations. In this case, the comparison is not independent, but the differences between the analyses and the observations are still interesting as they provide a measure of how the model fields have converged towards the observed state of the ocean.

3.3 Observational error sensitivity study

Before starting the numerical study, a year-long (2005) sensitivity test focusing on the observational position error for the Argo floats was carried out. The observational float position error was first set a horizontal distance of 500 m, based on the inherent Argos doppler-based positioning errors (∼250–1,500 m, cf. Menna and Poulain 2010) of the measurements. Root mean square (RMS) float position misfits were calculated and averaged over the test period and the Mediterranean, and the preliminary results indicated that adding trajectory assimilation makes the intermediate-current forecasts more accurate (CTRL, 30 km and TRAJ500 m: 25.6 km). However, the OceanVar scheme experienced convergence problems due to the “small” observational error, and it was found that a 2,000-m observational error made the iteration procedure in OceanVar more stable and yielded slightly better estimates of the modeled trajectories (TRAJ2,000 m: 24.4 km).

4 Results

4.1 Statistical analysis

4.1.1 The modeled float positions

The results from the sensitivity study suggested general improvements of the MFS float position forecasting skill. This was corroborated by the results from a more extensive statistical study on the differences between the modeled and observed float trajectories for the entire 3-year period in the Mediterranean Sea.

RMS float position misfits were calculated between the observations and the CTRL and TRAJ output, and these diagnostics were thereafter averaged over 2-week intervals (∼3 float cycles) in order to assure reliable statistics for the Mediterranean Sea. The results presented in Fig. 4 indicate that the use of OceanVar trajectory assimilation yield more accurate velocity field forecasts at intermediate depth in the Mediterranean. This improvement was quantified in terms of relative differences (%) between the 3-year averages of the CTRL and TRAJ RMS float position misfits (28 and 23.6 km, respectively) and found to be around 15%. Moreover, RMS differences were calculated between the TRAJ analyses and the observed float positions, and it was found that the quality of the trajectory analyses is on the order of the model horizontal grid resolution (∼5–7 km). The statistics were based on a fairly homogenous supply of trajectory observations (approximately two communicating floats per day), with a total number of 795, 822, and 682 Argo cycles in the years 2005, 2006, and 2007, respectively.

Fig. 4
figure 4

Two-week Mediterranean mean float position RMS misfits from the CTRL (black) and the TRAJ (red) experiments, as well as RMS differences in float positions between the TRAJ analyses and the observations (green)

Full size image

These results ultimately confirm that the trajectory model was implemented in a satisfactory manner into the OceanVar routines.

4.1.2 The sea level and mass fields

The velocity field corrections can influence the sea level, temperature, and salinity fields through the linear operators in Eq. 5, thus RMS misfits were calculated between the observations and the CTRL and TRAJ SSH, T and S fields. These results are presented in Fig. 5, and it can be deduced that the trajectory assimilation does not significantly influence the SLA, T (at 350 m) and S (at 350 m) RMS misfits in the Mediterranean during the 3-year period. The same conclusion could be drawn from the differences between the CTRL and TRAJ analyses and the observations, respectively. Hence, the quality of the modeled SSH, T, and S remain at the former levels. A possible explanation why the SSH fields are not significantly changed by the introduction of trajectory assimilation in OceanVar, might be the vastly larger amount of satellite-observed sea-level data in the Mediterranean compared with the number of Argo floats. The RMS differences for T and S were calculated using all Argo T–S observations, this, however, can have led to a masking of the corrections caused by trajectory assimilation, since not all T–S profiles were available with a “connecting” trajectory.

Fig. 5
figure 5

Two-week Mediterranean mean RMS misfits for the CTRL (black) and TRAJ (red) experiments, as well as the RMS differences between the CTRL (blue) and the TRAJ (green) analyses and observations during 2005–2007. Top panel, sea-level anomaly RMS differences (in cm); middle panel, temperature RMS differences (°C) at 350 m depth; lower panel, salinity RMS differences at 350 m depth

Full size image

4.2 The effects of assimilating Argo float trajectories

4.2.1 Cumulative effects of trajectory assimilation

The effects of assimilating Argo float trajectories on the velocity fields were studied 4 weeks after the numerical experiment had been launched, and the differences between the velocity fields on 25 January 2005 are given in Fig. 6. The speed differences between the CTRL and TRAJ background velocity fields were calculated as:

$$ \Delta{\bf u} = \sqrt{u_{\rm CTRL}^2+v_{\rm CTRL}^2} - \sqrt{u_{\rm TRAJ}^2+v_{\rm TRAJ}^2}, \label{eq:vel} $$
(9)

where u is the horizontal velocity field at 365 m depth. In this case, the differences between the two velocity fields are related to the corrections due to previously assimilated float positions in the TRAJ experiment (the observations on 25 January have not been assimilated yet, cf. Fig. 2). The negative (blue) areas indicate where the CTRL velocity field was weaker compared with the TRAJ velocity field, and vice versa for the positive (red) areas.

Fig. 6
figure 6

An example of propagation of velocity corrections (in cm − 1) in the model fields, here shown as: a speed differences between the CTRL and TRAJ background fields on 25 January 2005 and b speed differences between the CTRL and TRAJ background fields from the 25 to the 26 of January.The black dots mark in (a) the positions of all previously assimilated Argo floats and in (b) the Argo float positions assimilated on 25 January

Full size image

The largest velocity differences (generally around ± 3 cm − 1) were found in the vicinity of previously assimilated Argo floats (1–24 January), cf. the black dots in Fig. 6a. However, in many of these regions trajectory data had not been assimilated for several days, which indicates that the corrections made by OceanVar remain in the modeled velocity field memory at least on the order of the “Argo cycle” (Δt).

On January 25, surfacing coordinates from Argo floats were obtained in the AC (float 50769) and off the Spanish North-East coast (float 35504) as shown by the black dots in Fig. 6b. These positions were subsequently assimilated by OceanVar and in order to evaluate the impact of trajectory assimilation from one day to another, the speed differences between the CTRL and TRAJ background velocity fields on 25–26 January were calculated as detailed below:

$$ \begin{array}{rll} \Delta{\bf u} &\! =\! & \sqrt{[u_{\rm CTRL}(26)\!- \!u_{\rm CTRL}(25)]^2 \!+\! [v_{\rm CTRL}(26)-v_{\rm CTRL}(25)]^2} \nonumber\\[6pt] &&-\sqrt{[u_{\rm TRAJ}(26)\!-\!u_{\rm TRAJ}(25)]^2\!+\![v_{\rm TRAJ}(26)\!-\!v_{\rm TRAJ}(25)]^2}.\\ \label{eq:vel2} \end{array} $$
(10)

The largest changes in the velocity fields in Fig. 6b were found close to the two assimilated Argo floats, however, less evident alterations were found in most areas where trajectories had previously been assimilated. These findings further corroborate the assumption that velocity field corrections from previous assimilation cycles propagate both in time (on the order of ∼Δt) and across the model domain.

4.2.2 Local impact on the modeled velocity fields

The statistical analysis of the modeled trajectories gave at hand that the trajectory assimilation make the forecasted intermediate velocity fields more accurate. Here, the representativity of the modeled trajectories during the float drift at parking depth are to be examined. Observed float positions were compared with the CTRL- and TRAJ-modeled trajectories, and their corresponding velocity fields. Figure 7 offers an example of these comparisons for Argo float (50769) that was observed off the Algerian coast in January 2005. The CTRL and TRAJ trajectory predictions of the float drift, made on 20 January, were added to Fig. 7a, b.

Fig. 7
figure 7

Zoom in on the Algerian current velocity fields (in cm − 1) at 365 m depth: a the CTRL background fields on 25 January, b the TRAJ background fields on 25 January, and c the TRAJ background fields on 26 January. The observed “50769 float” positions on 10, 15, 20, and 25 January are marked as black (not assimilated) and red (assimilated) dots; the first position is indicated with a larger marker. The predicted “50769 float” trajectories were marked with red thick lines in (a) and (b). The analyzed “50769 float” trajectory was added in (c)

Full size image

Moreover, the TRAJ velocity background field on 26 January is provided in Fig. 7c, this to be interpreted as the “analyzed” fields for the previous day, since initial conditions had been corrected by the 25th January observations (cf. Fig. 2). The analyzed trajectory on 25 January has been superimposed on the “analyzed” velocity field.

The well-developed (but erroneous) eddy at approximately 2°40′ E, 37°40′ N, that causes the northward direction of the predicted trajectory in the CTRL velocity fields, is somewhat reduced in the TRAJ background velocity fields. However, both predicted trajectories (made on 20 January) failed to arrive to the observed float position on 25 January. This is probably addressable to the high day-to-day variability of the meandering AC system which may be responsible of a decrease in the reliability of the 5-day trajectory forecasts in this area.

The direction and strength of the “analyzed” velocity field is representative of the analyzed trajectory, and in good agreement with the observed float positions. In this case, these results suggest that trajectory assimilation can improve the quality of the modeled velocity fields both in terms of forecasts and analyses.

4.2.3 Vertical propagation of corrections

The temperature and salinity RMS misfits at intermediate depth in Fig. 5b, c showed no significant influence of the trajectory assimilation, however, changes in the T and S vertical distributions due to the altered velocity fields are likely to occur.

This plausible vertical propagation of OceanVar corrections in the velocity and mass fields is here examined in the framework of float 50769 in the AC. The CTRL and TRAJ zonal velocity fields as well as the differences between them along 3°23′ E are given in Fig. 8 on 25 January 2005. The black dots indicate the position of the float before surfacing on 25 January. It was noted that the largest differences in the velocity fields were found in the upper 300 m. Moreover, the positioning of the AC core was shifted ∼0.1° north and the westward surface flow around 37°30′–38° N was weaker in the TRAJ fields (2–6 cm − 1). This fact is further illustrated in Fig.  7b, where the small gyre at 3°40′ E, 37° N in the TRAJ fields forces a northward meandering of the AC.

Fig. 8
figure 8

Meridional transect of the Algerian Current (3°23′ E) showing the vertical propagation of corrections due to trajectory assimilation on 25 January 2005: a the CTRL zonal velocity fields (in cm − 1), b the TRAJ zonal velocity fields (in cm − 1), and c the differences between the CTRL and TRAJ zonal velocity fields. The black dots indicate the position of the “50769 float” before surfacing on 25 January

Full size image

Similar conclusions could be drawn for the mass fields along this transect. The most important discrepancies in the T fields (cf. Fig. 9) were related to the shift of the AC and found in the upper 300 m of the water column. The changes of the temperature gradients due to trajectory assimilation yielded T differences on the order of 0.2°C in the coastal area. The core of slightly warmer water (ΔT ~0.05°C at 400 m, probably of LIW origin) compared with the surrounding water mass at approximately 200–500 m depth, was to a large extent reduced in the TRAJ fields.

Fig. 9
figure 9

Meridional transect of the Algerian Current (3°23′ E) showing the vertical propagation of corrections due to trajectory assimilation on 25 January 2005: a the CTRL temperature fields (°C), b the TRAJ temperature fields (°C), and c the differences between the CTRL and TRAJ temperature fields. The black dots indicate the position of the “50769 float” before surfacing on 25 January

Full size image

The transect of the salinity distribution in Fig. 10 suggested changes in the upper 200 m due to the shift of the meandering AC, and the largest differences were found in the coastal region where the TRAJ fields were generally fresher (ΔS∼0.04–0.1) than the CTRL output. In this case, the trajectory assimilation tended to make the coastal current both colder and fresher, hence changing its buoyancy properties. Next, the accuracy of the alterations in the vertical is to be evaluated by comparing the in situ T–S profile from Argo float 50769 with the corresponding CTRL and TRAJ model profiles.

Fig. 10
figure 10

Meridional transect of the Algerian Current (3°23′ E) showing the vertical propagation of corrections due to trajectory assimilation on 25 January 2005: a the CTRL salinity fields, b the TRAJ salinity fields, and c the differences between the CTRL and TRAJ salinity fields. The black dots indicate the position of the “50769 float” before surfacing on 25 January

Full size image

The RMS misfits and the RMS differences between the analyses and the observed profile were calculated for both model outputs on 25 January. To establish if trajectory assimilation has a long-term influence on the vertical T and S structure, these diagnostics were also calculated as 3-year mean values using all available float profiles. The results are presented in Fig. 11, and it was found that trajectory assimilation appear to not have a significant influence on the T and S analyses on either synoptic or longer time scales. This statement holds true also for the 3-year mean T and S RMS misfits, however, the corresponding values based on the 25 January fields showed notable discrepancies between the CTRL and TRAJ outputs in the upper 350 m of the water column.

Fig. 11
figure 11

Vertical profiles of RMS differences between model values and observations. Upper panels show snapshots from float 50769 on 25 January 2005 of the a temperature and b salinity in the Algerian Current. Lower panels show 3-year mean c temperature and d salinity estimates in the Mediterranean Sea. The lines are color coded as follows: CTRL misfits (black), TRAJ misfits (red), CTRL analyses-observations residuals (blue), and TRAJ analyses-observations (green)

Full size image

It can be deduced from Fig. 11a, b that, on this occasion, the TRAJ experiment yields locally less accurate temperature and salinity background fields than CTRL. The TRAJ temperature and salinity proved to be approximately 0.01–0.02°C and ∼0.08 units less accurate in the upper 300 m, respectively, compared with the CTRL results. Although, as the temperature observational error is ∼0.01°C, it is noteworthy that the representativity of the temperature observation is at its limit. The salinity observational error is ∼0.01 units which implies that in this case, the trajectory assimilation has introduced errors locally in the salinity fields.

These overall results indicate that the assimilation of Argo trajectories not only affects the intermediate currents, but can indeed change the local properties of the state variables above parking depth. For example, a local change of ±6 cm − 1 in the zonal velocity fields corresponds to an increase/decrease in transport of approximately 0.13 Sv for a current cross-section area of 100 m depth and ∼0.2° in width, hence on the order of 10% of the net transport in the Algerian current (1.7 Sv according to Benzohra and Millot 1970).

4.3 Independent comparisons with surface drifter trajectories

Daily observations of surface drifter positions (lowpass filtered with a cut-off frequency of 36 h) were made available by OGS through the Mediterranean Surface Velocity Programme (MedSVP) for validation of the model velocity fields. The surface buoy of the drifter is attached to a drogue which is centered at 15 m depth (cf. Gerin et al. 2009) hereby making the observations representative of the near-surface circulation.

During 2006, a large number of Argo floats were located in the Eastern Mediterranean, and in particular within the Levantine basin. It was noted that Argo float 50761 was drifting in the Mersa–Matruh gyre in December 2006 and, at this moment, surface drifter 59770 was located in the proximity of this gyre at a distance of approximately 50 km from the float. Hence, an opportunity was found to assess the influence of (Argo float) trajectory assimilation on meso-scale gyre systems and the accuracy of the surface velocity fields.

The CTRL and TRAJ velocity fields at 15 and 365 m depth are provided in Fig. 12 on 16 December 2006. The Argo float positions on 1, 6, 11, and 16 December are marked with black (not assimilated) and red (assimilated) dots, while the surface drifter positions during 16–21 December are marked with black stars. The CTRL and TRAJ predicted float trajectories during 11–16 December were added in red to their corresponding 365-m velocity fields in the lower panels. Furthermore, the analyzed float trajectory was superimposed on the TRAJ 365-m velocity field in black.

Fig. 12
figure 12

Validation of TRAJ velocity fields using independent surface drifter data in the Levantine on 16 December 2006. Top panels show velocity fields (in cm − 1) at 15 m depth deduced from the a CTRL and b TRAJ experiments. The lower panels present the velocity fields at 365 m depth: c CTRL with forecasted “50761 float” trajectory marked in red and d TRAJ with the corresponding predicted (red) and analyzed (black) “50761 float” trajectories. The observed “50761 float” positions (1–16 December) are marked with black (not assimilated) and red (assimilated) dots. The surface drifter (59770) positions on 16–21 December are marked with black stars. The first float and drifter positions are indicated with larger markers

Full size image

The Mersa–Matruh gyre, located at 28°30′–29°30′ E and 33–34° N, is one of the strongest sub-basin scale features in the Eastern Mediterranean (Golnaraghi 1993; Taupier-Letage 2008), here showing typical upper thermocline velocities on the order of 30 cms − 1. The structure and strength of the gyre varies between the CTRL and TRAJ velocity fields (cf. Fig. 12c and d), although its horizontal dimensions of 1°×1° are roughly maintained. Both the CTRL and TRAJ float trajectory predictions underestimated the gyre velocity, while the analyzed trajectory was capable of reproducing the observed float position on 16 December.

The corrections of the intermediate velocity fields proved to propagate vertically towards the surface layers in the TRAJ experiment, this due to the barotropic linear model in V η (cf. Eq. 5). Moreover, these corrections altered the structure and the strength of the gyre at 15 m depth, and reduced the eddy meandering near 29°45′ and 34° N, cf. Fig. 12a, b. In this case, the velocity field corrections appear to have improved the representation of the upper 400 m velocity fields. During this period, the CTRL experiment forecasted a distinct northward flow around 30°16′E, 34°30′ N while the TRAJ results suggested a somewhat weaker north-eastward flow (shown only for 16 December). The observed drifter positions on 16–21 December indicated that the TRAJ velocity fields were in better agreement with the true state of the ocean on this occasion. Ultimately, this example shows that trajectory assimilation can in some cases improve the forecast quality of the velocity fields above the Argo float parking depth.

5 Conclusions

Basin-wide trajectory assimilation experiments have been undertaken for a 3-year period (2005–2007) in the Mediterranean Sea. It has been established that the extended OceanVar assimilation scheme is capable of improving the quality of the intermediate velocity fields based upon analyses. Indeed, statistical studies of the root mean square differences between the observed and forecasted float positions showed that ∼15% more accurate velocity fields are obtained when trajectory assimilation was performed. The accuracy of the trajectory analyses was found to be on the order of the model horizontal grid resolution.

Furthermore, it was demonstrated that OceanVar manages to minimize the differences between the velocity fields and the Argo float positions, this without introducing spurious values in the modeled velocity, sea-surface height and mass fields. Occasionally, the velocity field corrections caused local degeneration (in the upper 300 m) of the temperature and salinity quality on the order of the observational error. This problem was probably caused by an inaccurate representation of the local T and S error variability in the EOF-based covariance error matrix. However, this inconsistency proved to be statistically insignificant as the 3-year mean RMS misfit profiles indicated no negative impact on the vertical mass structure in the Mediterranean Sea in general.

MFS proved capable of yielding reasonable dimensions, strengths and directions of the AC boundary current and the Mersa–Matruh gyre in the trajectory assimilation experiment. Independent validation of the surface currents using drifter trajectories gave at hand that the corrections of the intermediate velocity fields can in some cases propagate vertically and improve the velocity field forecast quality at levels above the Argo float parking depth.

Occasionally, OceanVar experienced difficulties during the minimization procedure of the cost function, and it was noted that this tended to occur when float data had been missing, hence leaving gaps in the trajectory datasets. It could be of importance to study the effects on the model fields of assimilating noncontinuous trajectory time series, and to clarify the details of these irregularities within the scope of a future study.