source: altifloat/doc/ocean_modelling/Draft1.tex @ 202

%% 
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%% 
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%% 
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%% 
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%% 
%% Template article for Elsevier's document class `elsarticle'
%% with harvard style bibliographic references
%% SP 2008/03/01

\documentclass[preprint,12pt,authoryear]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times,authoryear]{elsarticle}
%% \documentclass[final,1p,times,twocolumn,authoryear]{elsarticle}
%% \documentclass[final,3p,times,authoryear]{elsarticle}
%% \documentclass[final,3p,times,twocolumn,authoryear]{elsarticle}
%% \documentclass[final,5p,times,authoryear]{elsarticle}
%% \documentclass[final,5p,times,twocolumn,authoryear]{elsarticle}

%% For including figures, graphicx.sty has been loaded in
%% elsarticle.cls. If you prefer to use the old commands
%% please give \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
\usepackage{multimedia}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{multimedia}
\usepackage{graphics}
\usepackage{bbm}
\usepackage{pstricks}
\usepackage{amsthm}
 \usepackage{graphics}
 \usepackage{pgf}
 \usepackage{wrapfig} 
 \usepackage{mathtools} 
\usepackage{caption}
\usepackage{subcaption}
%\usepackage{hyperref}



 
 
 
 
 
 \newcommand{\vectornorm}[1]{\left|\left|#1\right|\right|}
 \newcommand{\bo}[1]{\mathbf{#1}}
%\newcommand{\vectornorm}[1]{\left|\left|#1\right|\right|}
\newcommand{\dif}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\diff}[2]{\frac{\partial^2 #1}{\partial {#2}^2}}
\newcommand{\odif}[2]{\frac{d #1}{d #2}} 
\newcommand{\odiff}[2]{\frac{d^2 #1}{d #2^2}}

%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}

\journal{Ocean Modellling}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{Modelling Surface Currents in the Eastern Levantine Mediterranean Using Surface Drifters and Satellite Altimetry Data}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \address[label1]{}
%% \address[label2]{}

\author{Leila Issa, Julien Brajard, Milad Fakhri, Daniel Hayes, Laurent Mortier and Pierre-Marie Poulain. }

\address{ Department of Computer Science and Mathematics\\
        Lebanese American University\\
        Beirut, Lebanon\\
    Email: leila.issa@lau.edu.lb
}


\begin{abstract}
We present a new and fast method for blending altimetry and surface drifters data in the Eastern Levantine Mediterranean. The method is based on a variational assimilation approach for which the velocity is corrected 
by matching real drifters positions with a simple advection model simulation,
%after drifters data are matched to a simple advection model for their positions, 
taking into account the effect of the wind. The velocity correction is done in a time-continuous fashion by assimilating at once a whole trajectory of drifters in a time window, and by moving that window to exploit correlations between observations. We show that with few drifters, our method improves the estimation of velocity in two typical situations : an eddy between the Lebanese coast and Cyprus, and velocities along the Lebanese coast. 
\end{abstract}

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword
Altimetry \sep lagrangian data \sep data assimilations \sep drifters \sep surface velocity fields
%% PACS codes here, in the form: \PACS code \sep code
Eastern Levantine Mediterranean
%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

%% \linenumbers

%% main text
\section{Introduction}
\label{}
An accurate estimation of mesoscale to sub-mesoscale surface dynamics of the ocean is critical in several applications in the Eastern Levantine Mediterranean basin. For instance, this estimation can be used in the study of pollutant dispersion, which is important in this heavily populated region. A good knowledge of the surface velocity field is challenging, especially when direct observations are relatively sparse.

Altimetry has been widely used to predict the mesoscale features of the ocean resolving typically lengths on the order of $100$ km \citep{chelton2007global}. There are, however, limitations to its usage. It is inaccurate in resolving short temporal and spatial scales of some physical processes, like eddies, which results in blurring these structures. Further errors and inaccuracies occur near the coastal areas (within 20-50 km from land),
where satellite information is degraded; this is due to various factors such as land contamination, inaccurate tidal and geophysical
corrections and incorrect removal
of high frequency atmospheric effects at the sea surface \citep{caballero2014validation}.

To improve geostrophic velocities, especially near the coast, in situ observations provided by drifters can be considered (e.g. \citet{bouffard2008, ruiz2009mesoscale}). %[Bouffard et al., 2010; Ruiz et al., 2009] . 
Drifters follow the currents and when numerous, they allow for an extensive spatial coverage of the region of interest. They are not very expensive, easily deployable and provide accurate information on their position and other
environmental parameters \citep{lumpkin2007measuring}.

To illustrate the information provided by drifters data, we show in Figure~\ref{fig:cnrs} the real-time positions of three drifters launched south of Beirut on August 28 2013. These positions can be compared to the positions that would have been obtained if the drifters were advected by the altimetric velocity field. We observe that unlike the corresponding positions simulated by the altimetric field provided by \textit{Aviso} (see section~\ref{sec:aviso}), the drifters stay within 10-20 km from the coast. The background velocity field shown in that figure is the geostrophic field, averaged over a period of 6 days. The drifters' in situ data render a more precise image of the local surface velocity than the altimetric one; however, this only possible along the path following their trajectory. These types of data are therefore complementary. 


Numerous studies aim at exploiting the information provided by drifters (Lagrangian data) to assess the Eulerian surface velocity. A large number of these rely on modifying a dynamical model of this velocity by minimising the distance between observed and model simulated drifters trajectories. This variational assimilation approach, which was classically used in weather predictions \citep{courtier1994strategy,dimet1986variational}, was tested successfully in this context, by using several types of models for the velocity, such as idealised point vortex models \citep{kuznetsov2003method}, General Circulation Models with simplified stratification (e.g. \cite{kamachi1995continuous};  \cite{molcard2005lagrangian}; \cite{ozgokmen2003assimilation}, \cite{nodet2006variational}). However, in a lot of applications involving pollutant spreading such as the ones we are interested in, a fast diagnosis of the velocity field is needed in areas which are not a priori known in details. This prompts the need for a simple model that is fast and easy to implement, but that keeps the essential physical features of the velocity. In this work, we propose a new algorithm that blends geostrophic and drifters data in an optimal way. The method is based on a simple advection model for the drifters, that takes into account the wind effect and that imposes a divergence free constraint on the geostrophic component. 
The algorithm is used to estimate the surface velocity field in the 
Eastern Levantine basin, in particular in the region between Cyprus and the Syrio-Lebanese coast, a part of the Mediterranean basin that has not been so well studied in the literature before. 
 

%In the present work, we are interested in Near- Real time applications, simple model for combining drifters with altimetry. challenge remains to keep it physical  


%Realistic models involving highly complex physics are of course desirable but 

From the methodological point of view, combining altimetric and drifters data has been done using statistical approaches, with availability of extensive data sets. A common approach is to use regression models to combine geostrophic, wind and drifters components, with the drifters' velocity component being computed from drifters' positions using a pseudo-Lagrangian approach. When large data sets are available, this approach produces an unbiased refinement of the geostrophic circulation maps, with better spatial resolution. (e.g. \citet{poulain2012surface,menna2012surface,uchida2003eulerian,maximenko2009mean,niiler2003near,stanichny2015parameterization}). Another approach relies on variational assimilation: the work of  \citet{taillandier2006variational} is based on a simple advection model for the drifters' positions that is matched to observations via optimisation. The implementation of this method first assumes the time-independent approximation of the velocity correction, then superimposes inertial oscillations on the mesoscale field. 
These variational techniques had
led to the development of the so called ``LAgrangian Variational Analysis" (LAVA) algorithm,  initially tested and applied to correct model velocity fields using drifter trajectories \citep{taillandier2006assimilation,taillandier2008variational} and later 
 customised to several other applications such as model assimilation \citep{chang2011enhanced,taillandier2010integration} and more recently blending drifters and altimetry to estimate surface currents in the Gulf of Mexico \citet{berta2015improved}. 
 %applied it 
 %, where they also added a measure of performance consisting of skill scores, that compare 
%the separation between observed and hindcast trajectories to the observed absolute dispersion.


From the application point of view, blending drifters and altimetric data has been successfully applied to several basins, for example in: the Gulf of Mexico \citep{berta2015improved}, the Black Sea \citep{kubryakov2011mean,stanichny2015parameterization} the North Pacific \citep{uchida2003eulerian}, and the Mediterranean Sea \citep{taillandier2006assimilation,poulain2012surface,menna2012surface}. In \citet{menna2012surface}, there was a particular attention to the levantine sub-basin, where large historical data sets from 1992 to 2010 were used to characterise surface currents. 
The specific region which lies between the coasts of Lebanon, Syria and Cyprus, is however characterised by a scarcity of data. In the present work, we use in addition to the data sets used in \citet{menna2012surface}, more recent data from 2013 (in the context of Altifloat project) to study this particular region.  


Our contribution focuses on the methodological aspect, and it can be considered an extension of the variational approach used in \citet{taillandier2006variational}. The purpose is to add physical considerations to the surface velocity estimation, without making the method too complex, in order to still allow for Near Real Time applications. We do that by constraining the geostrophic component of that velocity to be divergence-free, and by adding a component due to the effect of the wind, in the fashion done in \citet{poulain2009}. We also provide a time-continuous correction by: (i) assimilating a whole trajectory of drifters at once and (ii) using a moving time window where observations are correlated. 

We show that with few drifters, our method improves the estimation of an eddy between the Lebanese coast and Cyprus, and predicts real drifters trajectories along the Lebanese coast. 




%One of the particularities of the Levantine bassin is an ensemble of 'tourbillons",? where?



This manuscript is organised as follows. We begin in Section 2. by describing the data sets used in the method and the validation process. In Section 3., we  provide a thorough description of the method including the definition of parameters involved, the model, and the optimisation procedure. We validate the method by conducting a twin experiment and a set of sensitivity analysis in Section 4., followed by two real experiments in Section 5., one in a coastal configuration  and another in an eddy.  
% le tourbillon au sud de Chypre et le tourbillon de Shikmona ˆ peu prs ˆ la mme latitude ˆ lÕouest de la c™te du Liban. Cet ensemble, parfois aussi appelŽ complexe tourbillonaire de Shikmona 12, est une structure permanente au sud de Chypre avec une variabilitŽ saisonnire. CÕest sur cet ensemble et sur son lien avec la topographie, notamment le mont sous-marin ƒratosthne sur lequel nous nous pencherons en particulier dans cette Žtude, comme dÕaprs la figure En effet, les monts sous-marins sont considŽrŽs comme une des causes des Žvolutions marines de mŽso-Žchelle.13. Dans le cas du mont ƒratosthne, il est possible que son influence sur la masse dÕeau environnante soit augmentŽe par son intŽraction avec les tourbillons quasi-permanents du complexe de Shikmona 14. La circulation dans le secteur du mont ƒratosthne est dominŽe par un anticyclone correspondant au tourbillon de Chypre. Here say that coastal configuration ??











\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/RealvsSimulatedTraj.pdf}
%\vspace{-30mm}
\caption{Altifloat drifters deployed on Aug 28 2013 (shown in  $-$x) versus  trajectories simulated using the \textit{Aviso} field (shown in $\tiny{--}$).  The velocity field shown is the  \textit{Aviso} field, averaged over 6 days.}
\label{fig:cnrs}
\end{center}
\end{figure}



\section{Data}
All the data detailed in this section were extracted for two target period : first from 25 August 2009 to 3 September 2009, and second from 28 August 2013 to 4 September 2013.
\subsection {\label{sec:aviso}Altimetry data}
Geostrophic surface velocity fields used as a background in the study were produced by Ssalt/\textit{Duacs} and distributed by \textit{Aviso}. Altimetric mission used were  Saral, Cryosat-2, Jason-1\&2. The geostrophic absolute velocity fields were deduced from Maps of Absolute Dynamic Topography (MADT) using the regional Mediterranean Sea product~\footnote{www.aviso.altimetry.fr}.

Data were mapped daily at a resolution of 1/8$^o$. Data were linearly interpolated every hour at the advection model time step.

\subsection{\label{sec:drifters}Drifters data}
Drifters were deployed during two target periods, 2 drifters were selected for the first period in 2009 and 3 in the second period in 2013. Table~\ref{tab:drifters} presents a summary of the 5 drifters used in this study. Drifter models were SVP with a drogue at a depth of 15m. Drifter positions were filtered with a low-pass filter in order to remove high-frequency current component especially inertial currents. The final time series were obtained by sampling every 6h. A more complete description of the drifters and the data processing procedure can be found in~\citet{poulain2009}.
\begin{table}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Project & Deploy Date & Lat & Lon & Last Date & Lat & Lon \\
\hline
NEMED & 29 Jul. 2009 & 31.90 & 34.42 & 28 Oct. 2009 & 34.1 & 31.77 \\
\hline
NEMED & 03 Aug. 2009 & 32.59 & 32.63 & 26 Dec. 2009 & 32.92 & 34.28 \\
\hline
Altifloat & 27 Aug. 2013 & 33.28 & 34.95 & 22 Sep. 2013 & 36.77 & 35.94 \\
\hline
Altifloat & 27 Aug. 2013 & 33.28 & 34.98 & 04 Sep. 2013 & 34.13 & 35.64 \\
\hline
Altifloat & 27. Aug. 2013 & 33.28 & 35.03 & 17 Sep. 2013 & 34.88 & 35.88 \\
\hline
\end{tabular}
\caption{\label{tab:drifters} List of drifters used to illustrate the methodology presented in this study, 
2 drifters deployed in 2009 (results were detailed in section~\ref{sec:cyprus}) and 
3 drifter were deployed in 2013 (results were detailed in sections~\ref{sec:lebanon})}
\end{table}

\subsection{Wind Data}
ECMW ERA-Interim wind products~\citep{Dee2011} were extracted in order to estimate wind-driven currents. Wind velocities closest to the surface (10 m) were extracted at a resolution of 1/8$^o$ at the same grid point as the \textit{Aviso} data. The data were resampled on a hourly time step.

Wind velocities were used to estimate wind-driven effect on drifter velocity.
The Eulerian velocity field being in the advection model(Eq.~\ref{advection}) is the sum of the geostrophic velocity and the wind induced velocity (Eq.~\ref{euler_vel}) given by the formula~\citep{poulain2009}:
\begin{equation}
\mathbf{U_{wind}} = 0.01exp(-28^oi)\times \mathbf{U_{10}}
\end{equation}
where $\mathbf{U_{wind}} = u_{wind}+iv_{wind}$ is the velocity induced by the wind and $ \mathbf{U_{10}} = u_{10}+iv_{10}$ is the wind velocity above the surface (10m) expressed as complex numbers.

\subsection {Model data}
Modeled surface velocity fields for September 2013 were used to calibrate the assimilation method presented in section~\ref{sec:method}. The model selected was the CYCOFOS-CYCOM high resolution model~\citep{zodiatis2003,zodiatis2008} that covers 
northeast Levantine basin (1km resolution, west and south boundaries extended to 31$^o$00'E and 33$^o$00'N and north and east reach land).
%the North-East Levantin Bassin 
%(31$^o$ 30’E - 36$^o$ 13’E  and 33$^o$ 30’N – 36$^o$ 55’N). 
The model forecasts were used without assimilation and were reinteroplated on a 1.8$^o$ grid point with an time step of one hour.
% The model forecast used for calibration purpose on September 2013.


\section{\label{sec:method}Method}

%%%%%%%%%%%%%%%%
\subsection{Statement of the problem}

We consider $N_f$ Lagrangian drifters released at time $t=0$ at various locations. These drifters provide their positions every $\Delta t$, over a period $[0,T_f]$.
Our objective is to determine an estimate of the two-dimensional Eulerian surface velocity field  \begin{equation}\notag
\mathbf{u}(x,y,t)=(u(x,y,t),v(x,y,t)) 
\end{equation}
characterized by a typical length scale $R$, given observations of the drifters' positions 
\begin{equation}
\bo{r}^{obs}_i(n\Delta t), \,\,\, i=1,2, \cdots, N_f, \,\,\, n=1, 2, \cdots N, \,\,\, \text{where}\,\,\, N \Delta t= T_f.
\end{equation}
The velocity shall be estimated on a specified grid with resolution of $1/8^{\circ}$ in both longitude and latitude, and in the time frame $[0, T_f].$ 

The estimation is done following a variational assimilation approach \citep{courtier1994strategy,dimet1986variational}, whereby the first guessed velocity, or background $\bo{u_b}$, is corrected by matching the observations with a model that simulates the drifters' trajectories. This correction is obtained using a sliding time window of size $T_w$, where we assume $\Delta t<T_w \leq T_L,$ and where $T_L$ is the Lagrangian time scale associated with the drifters in the concerned region. The background field is considered to be the sum of a geostrophic component (provided by altimetry) on which we impose a divergence free constraint, and a velocity component due to the wind. The details of this procedure are given in Section 3.3.




\subsection{Linearised model for Lagrangian data}

The position of a specific drifter $\mathbf{r}(t)=(x(t),y(t))$ is the solution of the non-linear advection equation
\begin{equation} \label{advection}
 \odif{\mathbf{r}}{t}=\mathbf{u}(\mathbf{r}(t),t), \,\,\,\,\bo{r}(0)=\bo{r}_0, \mathbf{u}(x,y,0)=\bo{u}_0.
 \end{equation}
This equation is integrated numerically, for example, using an Euler scheme. Since the drifters positions do not coincide with the Eulerian velocity's grid points, a spatial interpolation of $\mathbf{u}$ to these positions is needed. 


The observation operator, denoted it schematically by  $\bo{r}=\mathcal{M} (\bo{u}, \bo{r}),$ consists then of numerical advection and interpolation $\mathcal{I}$, and it is given by  
\begin{equation} \label{euler_advection}
\bo{r}(k\delta t)=\bo{r}((k-1)\delta t)+\delta t \, \mathcal{I}(\bo{u}((k-1)\delta t), \bo{r}((k-1)\delta t)), \,\,\,\,\, k=1,2, \cdots
\end{equation}
where $\delta t$ the time step of the scheme, typically a fraction of $\Delta t$. We choose bilinear interpolation 
\begin{align}
\mathcal{I}(\mathbf{u}, (x , y )) &= \mathbf{u}_1 + (\mathbf{u}_2 -\mathbf{u}_1 )\frac{(x - x_1)}{\Delta x} + (\mathbf{u}_3 -\mathbf{u}_1 )\frac{(y -y_1 )}{ \Delta y} \\ \notag
&+ (\mathbf{u}_1-\mathbf{u}_2-\mathbf{u}_3 + \mathbf{u}_4 )\frac{(x-x_1 )(y -y_1 )}{\Delta x \, \Delta y}, 
\end{align}
where
\begin{align} \notag
\mathbf{u}_1 &= \mathbf{u} (x_1 , y_1 ), \\  \notag
\mathbf{u}_2 &= \mathbf{u} (x_1 + 1, y_1 ), \\ \notag
\mathbf{u}_3 &=\mathbf{u}(x_1 , y_1 + 1), \\  \notag
 \mathbf{u}_4 &= \mathbf{u}(x_1 + 1, y_1 + 1). \notag
\end{align}
Here, $x_1=\left \lfloor{x}\right \rfloor $ is the floor function and $(x_1, y_1), (x_1 + 1, y_1), (x_1, y_1 + 1)$ and $(x_1 + 1, y_1 + 1)$ are the grid points which are nearest to $(x, y)$. 





Using the incremental approach \citep{courtier1994strategy}, the nonlinear observation operator $\mathcal{M}$ is linearised around a reference state. In a specific time window, we consider time independent perturbations $\delta \bo{u}$ on top of the background velocity field, that is
\begin{align}\label{totalR}
\bo{r}&=\bo{r^b}+\delta \bo{r} \\ \notag
\bo{u}&=\bo{u^b}+\delta \bo{u}. 
\end{align}
The linearised equations become 
\begin{align} \label{REquations}
&\bo{r^b}(k\delta t)=\bo{r^b}((k-1)\delta t)+\delta t \, \mathcal{I}\bigl(   \bo{u^b}((k-1)\delta t)), \bo{r^b}((k-1)\delta t     \bigr) ,\,\,\,\,\, \text{background}  \\  \notag
&\bo{\delta r}(k\delta t) = \bo{\delta r}((k-1)\delta t) + \delta t \, \{ \mathcal{I}(\bo{\delta u},\bo{r^b}((k-1)\delta t)) \\ \notag 
&+ \bo{\delta r}((k-1)\delta t) \cdot \partial _{(x,y)} \mathcal{I} \bigl(\bo{u^b}((k-1)\delta t),\bo{r^b}((k-1)\delta t)\bigr)\},  \,\,\, \text{tangent} 
\end{align}
where the drifters' positions are initialised with observations, and where $k=1,2,3, \cdots    \left \lfloor{T_w/\delta t}\right \rfloor .$ 
Here, $ \partial _{(x,y)} \mathcal{I}$ is the derivative of the interpolation operator with respect to $(x,y)$. 


The background velocity used in the advection of the drifters is the aggregate of a geostrophic component $\bo{u}_{geo}$ provided by altimetry and a component driven by the wind $\bo{u}_{wind}$, which is parametrised by two parameters as described in Section 2.3 \citep{poulain2009}. So we have  
\begin{equation}\label{euler_vel}
\bo{u}^b=\bo{u}_{geo}+\bo{u}_{wind}
\end{equation}
The data for both velocities are provided daily as described in the Data Section. This means that an interpolation in time of these velocities is needed.





\subsection{Algorithm for velocity correction}

We perform sequences of optimisations, where we minimise the following objective function with respect to the time independent correction $\delta \bo{u},$ in a specific time window $[0,T_w]$
\begin{equation}
\mathcal{J}(\delta \bo{u})=  \sum _{i=1}^{N_f} \sum_{m=1}^{\left \lfloor{T_w/\Delta t}\right \rfloor} \vectornorm{\bo{r}^{\,b}_{i}(\bo{u^b})+\delta \bo{r}_i(\delta \bo{u}) -\bo{r}_i^{\,obs}(m\Delta t) }^2 +\alpha_1 \vectornorm{ \bo{\delta u} }^2_{\bo{B}} +\alpha_2 \,div (\bo{u}_{geo}). 
\end{equation}

The first component of the objective function quantifies the misfit between the model 
obtained by iterations of \eqref{REquations}, and observations $\bo{r}^{\,obs}(m\Delta t)$. 
We highlight the dependence of $\bo{r}^b$ on the background velocity only, whereas $\bo{\delta r}$ depends on both background and correction.
The second component states that the corrected field is required to stay close to the background velocity. 
Here the $B$-norm is defined as $\vectornorm{\psi}^2_{\bo{B}} \equiv \psi^T \mathbf{B}^{-1} \psi,$ where $\bo{B}$ is the error covariance matrix.  This term serves the dual purpose of regularisation and information spreading or smoothing. To obtain $\bo{B}$, we use the diffusion filter method of \citet{weaver2001correlation}, where a priori information on the typical length scale $R$ of the Eulerian velocity can be inserted.
The parameter $\alpha_1$ represents the relative weight of this regularisation term with respect to the other terms.
The last component is a constraint on the geostrophic part of the velocity, required to stay divergence free. We note here that the total velocity may have a divergent component due to the wind. 


Inside a specific time window $[0, T_w]$, a whole trajectory of drifters contribute to give a constant correction in time $\bo{\delta u}.$ 
We now refine the method in a way to produce a smooth time-dependent velocity field in $[0, T_f]$, that accounts for temporal correlations in these various corrections $\delta u.$
One way to achieve this is to slide the window by a time shift $\sigma$ to obtain corrections $\bo{\delta u}_k$ in a window \[[k\sigma, k\sigma +T_w], \, \, k=0, 1, 2 \cdots. \]
The reconstructed velocity is then obtained as a superposition of the time dependent background field and weighted corrections   \[\bo{u}_{corrected}(t_i)=\bo{u^{b}}(t_i)+\sum_{k=0}^{N^i_w} w_k \bo{\delta u}_k.\]  A correction at a specific instant $t_i$ takes into account only $N^i_w$ windows sliding through $t_i$. The weight given to correction $\bo{\delta u}_k$ is chosen as \[w_k=\frac{1}{\vectornorm{t_i-tc_k}},\] where $tc_k$ is the centre of window $k,$ which means that it is inversely proportional to the ``distance" between time $t_i$ and the window's position.


We end this section by pointing out that we implement the algorithm described above in YAO~\citep{badran2008},
a numerical tool very well adapted to variational assimilation problems that simplifies the computation and implementation of the adjoint needed in the optimisation. 

The solution was found by using the M1QN3 minimizer linked with the YAO tool. The convergence of the assimilation on a typical time window of 24h is made 20 seconds on a sequential code compiled on a CPU Intel(R) Core(TM) at 3.40GHz.


\section{Twin Experiment}
In the twin experiment approach, the observations are simulated using a known velocity field, the ``true" velocity $\bo{u}_{true}.$ In this context, we are able to assess the validity of our approach by comparing the corrected, $\bo{u}_{corrected}$, and true fields. This is based on the time-dependent RMS error 
\begin{equation} \label {RMSError}
error (u, t)=\bigg(    \frac{\sum_{i,j} \big | \big |\bo{u}_{true} (i,j,t)-\bo{u} (i,j,t)\big | \big |^2}{\sum_{i,j} \big | \big |\bo{u}_{true} (i,j,t)\big | \big |^2 } \bigg)^{1/2},
\end{equation}
where $\big | \big |. \big | \big |$ refers the the $L_2$ norm of a vector, and $\bo{u}$ could be $\bo{u}_b$, giving the error before assimilation or $\bo{u}_{corrected}$, giving the error after assimilation. We shall also compare the trajectories of the drifters obtained by simulation with corrected and true velocity fields.



The configuration of our twin experiment is the following: we put ourselves in the same context as that of the real experiment conducted by the CNRS (refer to Altifloat drifters in Table~\ref{tab:drifters}), where the drifters are launched south of Beirut starting the end of August 2013. As shown in Fig.~\ref{fig:synth},  we deploy ``synthetic'' drifters in the region located  between 33.7 $^{\circ}$ and 34.25 $^{\circ}$ North and 34.9 $^{\circ}$ E and the coast. This is the same box in which the computation of the RMS error \eqref{RMSError} is done. The initial positions of the drifters shown in red coincide with the positions of the Altifloat drifters on September 1 2013.  The drifters' positions are simulated using a velocity field $\bo{u}_{true}$ obtained from the CYCOM model described in Section 2.4. The experiment lasts for a duration of $T_f=3$ days. In principle, nothing forbids us of conducting longer experiments, but in this coastal region, the drifters hit land after 3 days, as shown in Fig.~\ref{fig:synth}.
 The background velocity field is composed of the geostrophic component obtained from \textit{Aviso} and the wind component as described in the method section. These fields are interpolated to $\delta t$. A sensitivity analysis yields the optimal choice of $R=20$ km, which is consistent with the range of values found in the Northwestern Mediterranean \citep{taillandier2006variational}.
 
 
Using the relative RMS error before and after assimilation as a measure, we study the sensitivity of our method to the number of drifters $N_f$, the time sampling $\Delta t,$ the window size $T_w$ and to the moving parameter $\sigma.$ (Figures~\ref{fig:wsize},\ref{fig:mwin},\ref{fig:numb}).
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/RegionErroronAviso.pdf}
%\vspace{-30mm}
\caption{Region of RMS error computation for the twin experiment. Observations generated by CYCOM model starting on Sept 1st 2013 (for 3 days) are shown on top of the background field. The red locations correspond to Altifloat drifters' locations.}
\label{fig:synth}
\end{center}
\end{figure}
We first show the effect of the window size $T_w.$ On one hand, this parameter has to be within the Lagrangian time scale $T_L,$ estimated here to be $1-3$ days, but it also cannot be too large because we consider corrections that are time independent in each window. In Fig.~\ref{fig:wsize}, we show the results corresponding to various window sizes (fixing $N_f=14$ and $\Delta_t=2$ h), by displaying the relative RMS error (computed in the aforementioned box), before and after the correction. We first see that the error curve (after correction) tends to increase generally as time goes by. This is due to this special coastal configuration where the first three drifters hit the shore after $48$ h and also to the effect of the spatial filter, correcting the region behind the initial positions of the drifters. We then observe that the optimal window size for this configuration is $24$ h, which is within the range mentioned above. The error in this case is reduced by almost a half from the original one. We mention here that with window sizes close to three days, the algorithm becomes ill conditioned, which is expected due to the fact that the correction is fixed in a specific window, as mentioned before. 
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{./fig/Wins_optshift_dt1_f14_tf72.pdf}
%\vspace{-30mm}
\caption{The effect of the window size. Error before correction is shown with a solid line. Errors after are shown with symbols for several window sizes. Here $N_f=14$ and $\Delta t=2$ h}
\label{fig:wsize}
\end{center}
\end{figure}
We then show the effect of shifting the window by varying the parameter $\sigma.$ The window size here is fixed to $T_{w}=24$ h, time sampling to $\Delta t=2$ h, and $N_f=14.$ We start with $\sigma=0$, which amounts to doing separate corrections, then slide every $\sigma=12, 8, 6$ h. In Fig.~\ref{fig:mwin}, we show the results by displaying the relative RMS error before and after the correction. 
We observe that if the corrections are done separately, the correction is not smooth; in fact smaller values of $\sigma$ yield not only smoother, but better corrections \textcolor{red}{put an interpretation, we exploit the correlation...}. 
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{./fig/Shifts_win24_dt1_f14_tf72.pdf}
%\vspace{-25mm}
\caption{The effect of the moving window: smaller shifts $\sigma$ yield smoother and better corrections. Here $N_f=14$ and $\Delta t=2$ h.
 }
\label{fig:mwin}
\end{center}
\end{figure}
The effect of the number of drifters is shown next in Fig.~\ref{fig:numb}. Respecting coverage, we start with $N_f=14$ drifters (positioned as shown in Fig.~\ref{fig:synth}), then reduce to $N_f=10,6,3.$ Naturally more drifters yield a better correction but we notice that even with three drifters, the error is still reduced by $20\%$ and much more so close to the beginning of the experiment. We also show in this figure the effect of removing the drifters that fail before the end of the experiment: the corresponding error curve is shown in the dashed curve of Fig.~\ref{fig:numb}, and it is evenly distributed in time as expected. Finally, we show the effect of the time sampling parameter $\Delta t$ of the observations in Fig.~\ref{fig:time}. Curves after correction correspond to $\Delta t=6, 4$ and $2$ hours and as we see from the figure, the difference between these cases is not too large. The realistic scenario of  $\Delta t=6$ h still yields a good correction.
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{./fig/Nfs_win24_dt1_tf72_nofail.pdf}
%\vspace{-30mm}
\caption{The effect of the number of drifters. More drifters yield better corrections but corrections are possible with 3 drifters only. The dashed line shows the effect of just taking drifters that do not hit the shore before the end of the experiment. Here $T_w=24$ h and $\Delta t=2$ h.}
\label{fig:numb}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{./fig/Dts_win24_f14_tf72.pdf}
%\vspace{-30mm}
\caption{The effect of the time sampling $\Delta t$ of the observations. Here $T_w=24$ h, and $N_f=14.$ The realistic scenario of  $\Delta t=6$ h is not too far from the smallest $\Delta t=2$ h.}
\label{fig:time}
\end{center}
\end{figure}

For the twin experiment with the optimal choice of parameters ( $T_w=24$ h, $\sigma=6$ h, $N_f=14$ and $\Delta t=2$), we now show the trajectories of the drifters simulated with the corrected velocity field on top of the ``true" observations. We also compare background and corrected fields in the region of interest. In Fig.~\ref{fig:lerror}, we display the point-wise $L_2$ error, defined as 
\begin{equation} \label {L2Error}
error (u,i,j,t)=\big | \big |\bo{u}_{true} (i,j,t)-\bo{u} (i,j,t) \big| \big|,
\end{equation}
between the true field and either the background or corrected fields. In the left panel, we show that error between the background and true fields, averaged in time over the duration of the experiment. The right side shows this same error after correction. On top of that, we observe the excellent agreement between the positions of the drifters simulated with the corrected field and the ``true" observations. Next, the correction in terms of direction is shown in Fig.~\ref{fig:cerror}: we display the cosine of the angle between the background and true field on the left side versus the  cosine of the angle between the corrected and true fields on the right. Note that a cosine of one indicates a strong correlation in direction between the two fields. We see this strong correlation between true and corrected fields by observing how the blue color (left pannel of Fig.\ref{fig:cerror} ) turns into deep red (right pannel of Fig.\ref{fig:cerror} ) in the region where the drifters were deployed. Finally, in Fig. \ref{fig:summary}, we show the actual current maps before and after correction. We clearly see that the drifters corrected the poorly simulated coastal meander.


\begin{figure}[htdp]

\begin{subfigure}{0.6\textwidth}
\includegraphics[width=1.08\linewidth]{./fig/Before_L2pointw_win24_MEAN_color.pdf}
\end{subfigure}%
%\hspace{-20mm}
\begin{subfigure}{0.6\textwidth}
 % \centering
 \includegraphics[width=1\linewidth]{./fig/After_L2pointw_win24_MEAN_color.pdf}
\end{subfigure}
\vspace{-25mm}
\caption{Point-wise $L_2$ error averaged over time, before (left) and after (right) correction. In the right frame, drifters' positions obtained by simulation with corrected field (magenta) versus ``true" observations(black) are shown on top of the error.}
\label{fig:lerror}
\end{figure}


\begin{figure}[htdp]
%\centering
\begin{subfigure}{0.6\textwidth}
 % \entering
  \includegraphics[width=1.05\linewidth]{./fig/Before_CORRANGLEpointw_win24_MEAN_color.pdf}
%  \label{bla}
\end{subfigure}%
%\hspace{-10mm}
\begin{subfigure}{0.6\textwidth}
 % \centering
  \includegraphics[width=1.05\linewidth]{./fig/After_CORRANGLEpointw_win24_MEAN_color.pdf}
%  \label{blo}
\end{subfigure}
\vspace{-5mm}
\caption{Correction in terms of direction. Left: $\cos{(\bo{u}_{b},\bo{u}_{true})}$, right: $\cos{(\bo{u}_{corrected},\bo{u}_{true})}.$ }
\label{fig:cerror}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.4]{./fig/summary_twin.pdf}
\caption{Background velocity field (blue) versus corrected velocity field (red) for the twin experiment with the optimal choice of parameters.}
\label{fig:summary}
\end{center}
\end{figure}


\section{Experiments with Real Data}


The methodology described in section~\ref{sec:method} was applied to two case studies : one along the Lebanese coast and one in an eddy south-east of Cyprus.
\subsection{\label{sec:lebanon}Improvement of velocity field near the coast}

Three drifters were launched on August 28 2013 from the South of Beirut, at the positions shown in circles in Fig.~\ref{fig:leb1}. They provide their position very $\Delta t= 6$ h and stay within 20 km of the coast for the duration of the experiment. 
The experiment considered here lasts for six days (a time frame where the three drifters are still spatially close before two of them hit the shore). The window size is $T_w=24$ h. The smoothing parameter $\sigma=6$ h. 
Fig.~\ref{fig:leb1}, shows the trajectories simulated with corrected field on top of the observed ones, 
in good agreement, even for small scale structures near the coast.
Average correction over 6 days are shown on the figure, but the actual corrections are time-dependent.

As expected, the velocity field is modified in the neighbourhood of the drifters trajectories. It can be noticed that the main effect of the correction is to increase the velocity parallel to the coast, and decrease the velocity normal to the coast. The background field was determined using altimetric data and is expected to have significant bias close to the coast~\citep{bouffard2008}, and the consequence is that the method is able to correct some of this bias.

To validate more quantitatively the corrected velocities, sensitivity study was carried out. Only two drifters (the easternmost magenta drifter and the westernmost black drifter) were assimilated in order to correct the velocity field. The third drifter is used only to validate the corrected field by comparing its actual trajectory with the simulated trajectory using the velocity field.

Figure~\ref{fig:lebzoom} shows the results of this experiment. The real drifter trajectory (empty circle with thin line) was compared to the simulated trajectory using either the background field (bold cyan line) or the corrected field (bold green line). 

It can be noticed that the trajectory was greatly improved using the corrected field. It shows that the corrected field can be used to simulate realistic trajectories in the neighbourhood of the assimilation positions, even in a coastal region. 
It can be a decisive point for application such as pollutant transport estimation.

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_6days_average.pdf}
%\vspace{-30mm}
\caption{\label{fig:leb1} Prediction of the positions of 3 Altifloat drifters, launched on August 28 2013. $T_f=6$ days.  $T_w=24$ h and $\sigma=6$ h. Positions of drifters simulated with corrected field (cross markers) are shown on top of observed positions (circle markers). Corrected field is shown in red whereas background field is shown in blue. }
\end{center}
\end{figure}


\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_bmtogreen_2days_average_zoom.pdf}
%\vspace{-30mm}
\caption{\label{fig:lebzoom} Prediction of the position of the green drifter using the observed black and magenta drifters. $T_f=2$ days.  $T_w=24$ h and $\sigma=6$ h. Position of the green drifter simulated with corrected field is shown in green squares, on top of observed position shown in light green circles. Compare to the position of the drifter obtained with background field only, shown in cyan. Corrected field is shown in red whereas background field is shown in blue.}
\end{center}
\end{figure}



\subsection{\label{sec:cyprus}Improvement of velocity field in an eddy}
In the context of the Nemed deployment (see section ~\ref{sec:drifters}), we selected two drifters trajectories from 25 August 2009 to 3 September 2009. Assimilating the successive positions of the drifters every six hours, the \textit{Aviso} velocity field was corrected.

In this experiment the window size $T_w$ was chosen to be 72 hours as the velocity field was more stable in this case than in coastal areas. The shifting of the time window was of 18 hours.

In Fig.~\ref{fig:eddy-velocity}, the trajectory of the drifters were represented in gray, the mean \textit{Aviso} surface geostrophic velocity field in blue and the mean corrected geostrophic field in red. 

The real trajectory of the drifters and the simulated trajectory using the total corrected field (sum of corrected field in red and the wind-induced velocity) are very close. 
The mean position error 
%(expressed in arc length) 
is 
0.96 km
%$8.6\times 10^{-3}$ degrees 
with a maximum of 
%$0.06$ degrees. 
6.7km.
The real trajectory and simulated trajectory would be indiscernible in Fig.~\ref{fig:eddy-velocity}. 



\begin{figure}[h]
\centering
\includegraphics[scale =0.6]{./fig/Eddy_velocity.png}
\caption{\label{fig:eddy-velocity} Corrected surface velocity field (in red) compared to \textit{Aviso} background field (in blue). The assimilated drifter trajectories are represented in gray. The North-West coast in the figure is Cyprus.}
\end{figure}

In this case, it can be seen that the drifter trajectories were situated in an eddy. The \textit{Aviso} field is produced by an interpolation method which tends to overestimate the spatial extent of the eddy and underestimate the intensity. In order to estimate the effect of the assimilation on the eddy characteristics, we computed the Okubo-Weiss parameter~\citep{isern2004} on the mean velocity fields before correction (background) and after correction. Eddies are characterized by a negative Okubo-Weiss parameter, the value of the parameter is an indicator of the intensity of the eddy. Results are shown in Fig.~\ref{fig:okubo-weiss}. As expected, it can be noticed that the Okubo-Weiss parameter had greater absolute values and a slightly smaller spatial extent which indicated a improvement of the Aviso processing bias. This results constitutes a validation of the assimilation method presented in this paper showing that eddies were better resolved after assimilating drifter trajectories.

\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{./fig/okubo_weiss_aviso.png}
\includegraphics[scale=0.5]{./fig/okubo_weiss_analyse.png}
\caption{\label{fig:okubo-weiss} Okubo-Weiss parameter calculated on background field (upper panel) and corrected field (lower panel). The negativity of this parameter characterised eddies, and the absolute value corresponds to the intensity of the eddy. It can be noticed that eddy is smaller in size and more intense after the correction process.}
\end{figure}

\section{Conclusion}
We presented a simple and efficient algorithm to blend drifter lagrangian data with altimetry Eulerian velocities in the Eastern Levantine Mediterranean. The method has a cheap implementation and is quick to converge, so is is well fitted for near-real time applications. Assimilating two successive drifter positions produces a correction of the velocity field within a radius of 20km and for approximatively 24h before and after the measurement.

This algorithm was able to correct some typical weakness of alimetry fields, in particular the estimation of velocity near the coast and accurate estimations of eddies dimensions and intensity.

\section{Acknowledgement}
The altimeter products were produced by Ssalto/Duacs and distributed by Aviso, with support from Cnes (http://www.aviso.altimetry.fr/duacs/).

Wind data were produced by ECMWF and downloaded from (http://apps.ecmwf.int/datasets/data/interim-full-daily/).

This work was funded by the ENVI-MED program in the framework of the Altifloat project.

 The Lebanese CNRS funded through "CANA project", the campaigns of drifters' deployment using the platform of the Lebanese research vessel "CANA-CNRS". These MetOcean Iridium drifters (SVP) were provided through the Istituto Nazionale di Oceanografia e di Feofisica Sperimentale (OGS), Italy and LOCEAN institute of "Pierre et Marie Curie University", France
%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% If you have bibdatabase file and want bibtex to generate the
%% bibitems, please use
%%
\section{Bibliography}

\bibliographystyle{elsarticle-harv} 
\bibliography{mybib.bib}

%% else use the following coding to input the bibitems directly in the
%% TeX file.

%\begin{thebibliography}{00}

%% \bibitem[Author(year)]{label}
%% Text of bibliographic item

%\bibitem[ ()]{}

%\end{thebibliography}
\end{document}

\endinput
%%
%% End of file `elsarticle-template-harv.tex'.