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

%% 
%% 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}



 
 
 
 
 
 \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{Modeling Surface Currents in the Eastern Levantine Mediterranenan 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, Pierre-Marie Poulain, Laurent Mortier, Dan Hayes, Milad Fakhri}

\address{}

\begin{abstract}
%% Text of abstract

\end{abstract}

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword

%% PACS codes here, in the form: \PACS code \sep code

%% 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{}
Recurrent marine pollution like the ones observed near the coastal regions of Lebanon is a major threat to the marine environment and halieutic resources. Polluting agents do not only have a immediate local effect but they are also transported through ocean currents to deep waters far away from the coast thereby having a long term, large scale effect. Modeling ocean currents is thus a major topic that mobilizes multidisciplinary researchers.  

The objective of this project is to combine two types of observations (drifters and satellite data) in order to reconstruct mesoscale features of the surface currents in the Levantine Mediterranean including coastal regions. The figure below shows positions of drifters and altimetry data in the Eastern Mediterranean from the 27th to the 29th of November 2009: the altimetric data gives the velocity field in the whole basin with a resolution of $10$ km smoothing considerably some mesoscale features with errors especially near the coastal regions. The buoys on the other hand allow a better reconstruction of the velocity field but only along their trajectories. These two types of data are therefore complimentary and we can then formulate the problem as follows: find the optimal velocity field from these two sources of information.   

The mathematical problem consists of ``inverting'' the Lagrangian trajectories (Drifters positions) to correct the Eulerian velocity field given my the altimetric data. The literature on the subject is abundant: ranging from optimal interpolation  \cite{ozgokmen2003assimilation}, to statistical approaches based on extended Kalman filter methods \cite{salman2006method}, to variational techniques \cite{kamachi1995continuous}, \cite{nodet2006variational}, \cite{taillandier2006variational}.
We propose to use a variational assimilation technique where we correct the Eulerian velocity field by minimizing the distance between a model solution and observations of the buoys. We use the incremental approach described in \cite{courtier1994strategy}, \cite{dimet2010variational}. The advantage of this method is that while being midway in complexity between OI (lower) and Kalman filter (higher), it permits a simple mathematical formulation of the dynamical constraints we wish to impose. 


\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/RealvsSimulatedTraj.pdf}
\vspace{-30mm}
\caption{Real CNRS observations of 3 drifters starting Aug 28 2013 -x , vs simulated -- with aviso same date. On top of averaged Aviso for 6 days starting aug 28  }
\end{center}
\end{figure}



\section{Data}
\subsection {Altimetry data}
\subsection{Drifter data}
\subsection {Model data}
%Wind + Dan


\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 each $\Delta t$, over a period $[0,T_f]$.
Our objective is to determine an estimate of the two-dimensional Eulerian velocity field  \begin{equation}\notag
\mathbf{u}(x,y,t)=(u(x,y,t),v(x,y,t)) 
\end{equation}
characterized by a length scale $R$ [Refs], 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 $\Delta x,$ in the time frame $[0, T_f].$ 

The estimation is done following a variational assimilation approach [Refs], 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 done using a sliding time window of size $T_w$, where we assume $\Delta t<T_w<T_L,$ and where $T_L$ is the Lagrangian time scale associated with the drifters in the concerned region [Refs]. 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 due the wind [Refs]. The details of this procedure are given in Section 3.3.




\subsection{Linearized 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 using an Euler scheme for example. Since the drifters positions may not coincide with the Eulerian velocity's grid points, a spatial interpolation of $\mathbf{u}$ to these positions is also needed. 


The observation operator, denoted it schematically by  $\bo{r}=\mathcal{M} (\bo{u}, \bo{r}),$ consists then of numerical advection and interpolation, and it is given by  
\begin{equation} \label{euler_advection}
\bo{r}(k\delta t)=\bo{r}((k-1)\delta t)+\delta t \, interp(\bo{u}((k-1)\delta t), \bo{r}((k-1)\delta t)), \,\,\,\,\, k=1,2, \cdots
\end{equation}
%Let us denote by $\bo{r}_k = (x_k ,y_k)$ the position of the drifter at time $t_k=k\delta t$ , $\mathbf{u}_k$ the Eulerian velocity at position $\mathbf{r}_k$ and time $t_k$, and
where $\delta t$ the time step of the scheme, typically a fraction of $\Delta t$. We choose bilinear interpolation 
\begin{align}
interp(\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)$. 



%This has the advantage that the cost function becomes quadratical; it has a unique minimum and this minimum is assumed to be close to that of the full non-quadratical cost function (below). 

Using the incremental approach [Refs], the nonlinear observation operator $\mathcal{M}$ is linearized 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 linearized equations become 
\begin{align} \label{REquations}
&\bo{r^b}(k\delta t)=\bo{r^b}((k-1)\delta t)+\delta t \, interp\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 \, \{ interp(\bo{\delta u},\bo{r^b}((k-1)\delta t)) \\ \notag 
&+ \bo{\delta r}((k-1)\delta t) \cdot \partial _{(x,y)} interp \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 initialized with observations, and where $k=1,2,3, \cdots    \left \lfloor{T_w/\delta t}\right \rfloor .$ 
Here, $ \partial _{(x,y)} interp$ is the derivative of the interpolation function 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 parametrized by two parameters as described in PPM. So we have  
\[\bo{u}^b=\bo{u}_{geo}+\bo{u}_{wind}
\]
The data for both velocities are provided daily as described in the Data Section. This means that an interpolation in time of these velocities may be needed.



%\subsection{Discussion on the choice of $\delta t$ }
%The observations $\bo{r}^{obs}$ are available to us for each $\Delta t$, that is we have $\bo{r}^{obs}(i \Delta t), \,\,  i=0,1,2, \cdots $. Let us denote by $\Delta x$ the typical distance traveled by a buoys during $\Delta t$ and let $\Delta s$ denote the spatial resolution of the correction $\delta \bo{u}$. We have two scenarios: 
%\begin{enumerate}[(i)]
%\item If $\Delta s >> \Delta x$, then we can take $\delta t \sim \Delta t$ (that should simplify equation \eqref{Evolve}) 
%\item If $\Delta s << \Delta x$, we need to iterate many $\delta t$ to arrive at an observation. We take $N \delta t=\Delta t$
%\end{enumerate}

\subsection{Algorithm for assimilation}

We perform sequences of optimizations, where we minimize 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 regularization and information spreading or smoothing. To obtain $\bo{B}$. we use the diffusion filter method of [Weaver and Courtier, 2001], 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 regularization 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. More here?.We note here that the total velocity may have a divergent component due to the wind. 


Inside a time window $[0, T_w]$, we assimilate a whole trajectory of drifters, which gives a constant correction $\bo{\delta u}.$ We refine the method by reconstructing a smooth time dependent velocity inside $[0, T_f]$: we achieve this by sliding the window by a time shift $\sigma$ to obtain other corrections $\bo{\delta u}_k$ in $[k\sigma, k\sigma +T_w], \, \, k=0, 1, 2 \cdots$. The reconstructed velocity is then obtained as \[\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 at specific instant $t_i$ only takes into accounts $N^i_w$ windows sliding through $t_i$ and weight $w_k$ given to correction $\bo{\delta u}_k$ is inversely proportional to the ``distance" between time $t_i$ and the window's position. 
%Here $i=1,2, \cdots T_f/\delta t$




DIAGRAM IS USEFUL!

\subsection{Implementation of the assimilation method in YAO}

\section{Twin Experiment}
In the twin experiment approach, the observations are simulated using a known velocity field referred to as 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 the true fields. This will be based on the time-dependent RMS error 
\begin{equation} \label {RMSError}
error (u, t)=\bigg(    \frac{\sum_{i,j} |\bo{u}_{true} (i,j,t)-\bo{u} (i,j,t)|^2}{\sum_{i,j} |\bo{u}_{true} (i,j,t)|^2 } \bigg),
\end{equation}
where $\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 the corrected velocity field and the observed trajectories. 



In our experiment, we obtain the true velocity field from a dynamic model \textcolor{red}{Need details from Dan}. We start the optimisation with the background velocity field as described in the method section. The region we study is located off the shore of Beirut Explain that CNRS drifters were launched south of Beirut, so we put drifters close to there. As shown in Fig. 2, we compute the error in a box between 33.7 and 34.25 North and 34.9 E and the coast. There we ``deploy" a maximum of 14 drifters. The duration of the experiment is 3 days (nothing forbids of going longer, but some drifters hit the coast). We shall study the density with respect to the number of drifters (respecting coverage), to the time sampling $\Delta t,$ to the window size $T_w$ and to the moving parameter $\sigma.$



\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/RegionErroronAviso.pdf}
\vspace{-30mm}
\caption{Region for RMS Error Computation surrounding the observations generated by Dan's model of Sept 1st 2013. on top of the average over 3 days Aviso Field starting  Sept 1st 2013, }
\end{center}
\end{figure}



\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/Shifts_win24_dt1_f14_tf72.pdf}
\vspace{-30mm}
\caption{The effect of the moving window for a window of size $24$hs with 14 drifters and $dt=2hs$ }
\end{center}
\end{figure}


\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/Wins_optshift_dt1_f14_tf72.pdf}
\vspace{-30mm}
\caption{The effect of the window size for 14 drifters and $dt=2$hs}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/Nfs_win24_dt1_tf72.pdf}
\vspace{-30mm}
\caption{The effect of the number of drifters}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/Dts_win24_f14_tf72.pdf}
\vspace{-30mm}
\caption{The effect of the time sampling}
\end{center}
\end{figure}





\begin{figure}[htdp]
%\centering
\begin{subfigure}{0.55\textwidth}
 % \entering
  \includegraphics[width=1.2\linewidth]{./fig/Before_L2pointw_win24_MEAN_color.pdf}
%  \caption{Pressure drop.}
  \label{bla}
\end{subfigure}%
\hspace{-10mm}
\begin{subfigure}{0.55\textwidth}
 % \centering
  \includegraphics[width=1.2\linewidth]{./fig/After_L2pointw_win24_MEAN_color.pdf}
%  \caption{Mass flowrate.}
  \label{blo}
\end{subfigure}
\vspace{-30mm}
\caption{the point-wise L2 error before and after in m/s}
\label{fblo}
\end{figure}


\begin{figure}[htdp]
%\centering
\begin{subfigure}{0.55\textwidth}
 % \entering
  \includegraphics[width=1.2\linewidth]{./fig/Before_CORRANGLEpointw_win24_MEAN_color.pdf}
%  \caption{Pressure drop.}
  \label{bla}
\end{subfigure}%
\hspace{-10mm}
\begin{subfigure}{0.55\textwidth}
 % \centering
  \includegraphics[width=1.2\linewidth]{./fig/After_CORRANGLEpointw_win24_MEAN_color.pdf}
%  \caption{Mass flowrate.}
  \label{blo}
\end{subfigure}
\vspace{-30mm}
\caption{the point angle correlation error before and after in m/s}
\label{fblo}
\end{figure}




\section{Experiment 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{Improvement of velocity field near the coast}
%lebanese drifters
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_6days_average.pdf}
\vspace{-30mm}
\caption{Real Exp starting Aug 28 for 6 days, window=24h, move=6hs, }
\end{center}
\end{figure}


\begin{figure}[htbp]
\begin{center}
\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_bmtogreen_2days_average_zoom.pdf}
\vspace{-30mm}
\caption{Real 2-1 Exp starting Aug 28 for 2 days, window=24h, move=24hs, }
\end{center}
\end{figure}

\subsection{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 AVISO velocity field was corrected.

In Fig.~\ref{fig:eddy-velocity}, the trajectory of the drifters were represented in gray, the AVISO surface geostrophic velocity field in blue and the corrected geostrophic field in red. In this case, it can be seen that the drifter trajectories were situated in an eddy. The AVISO field is produced by an interpolation method which tends to overestimate the spatial extent of the eddy and underestimate the intensity.
\begin{figure}[h]
\centering
\includegraphics[scale =0.6]{./fig/Eddy_velocity.png}
\caption{\label{fig:eddy-velocity} Corrected field (in red) compared to AVISO background fields (in blue). The assimilated drifter trajectories are represented in gray. The North-West coast in the figure is Cyprus.}
\end{figure}

%% 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
%%
%%  \bibliographystyle{elsarticle-harv} 
%%  \bibliography{<your bibdatabase>}

%% 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'.