Changeset 172 for altifloat/doc


Ignore:
Timestamp:
09/11/15 13:49:28 (9 years ago)
Author:
leila_ocean
Message:

on-going changes

Location:
altifloat/doc/ocean_modelling
Files:
16 added
10 edited

Legend:

Unmodified
Added
Removed

  • altifloat/doc/ocean_modelling/Draft1.tex

    r171 r172  
    5151 \usepackage{wrapfig}  
    5252 \usepackage{mathtools}  
     53\usepackage{caption} 
     54\usepackage{subcaption} 
     55 
     56 
     57 
     58  
     59  
     60  
     61  
     62  
    5363 \newcommand{\vectornorm}[1]{\left|\left|#1\right|\right|} 
    5464 \newcommand{\bo}[1]{\mathbf{#1}} 
    … …  
    126136\section{Introduction} 
    127137\label{} 
     138Recurrent 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.   
     139 
     140The 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.    
     141 
     142The 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}. 
     143We 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.  
     144 
     145 
     146\begin{figure}[htbp] 
     147\begin{center} 
     148\includegraphics[scale=0.5]{./fig/RealvsSimulatedTraj.pdf} 
     149\vspace{-30mm} 
     150\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  } 
     151\end{center} 
     152\end{figure} 
     153 
     154 
     155 
    128156\section{Data} 
    129 \subsection {altimetric data} 
    130 \subsection{Drifter data} 
    131 \subsection {Model data} 
    132 %Wind + Dan 
    133  
    134 \section{\label{sec:method}Method} 
     157 
     158 
     159 
     160 
     161\section{Method} 
     162 
     163%%%%%%%%%%%%%%%% 
     164\subsection{Statement of the problem} 
     165 
     166We consider $N_f$ Lagrangian drifters released at the time $t=0$ at various locations. These drifters provide their positions each $\Delta t$, over a period $[0,t_f]$, where $t_f=N \Delta t.$  
     167Our objective is to determine an estimate of the Eulerian two-dimensional velocity field  \begin{equation}\notag 
     168\mathbf{u}(x,y,t)=(u(x,y,t),v(x,y,t))  
     169\end{equation} 
     170on a specified grid, in a time frame $[0, t_f],$ given the observations  
     171\begin{equation} 
     172\bo{r}^{obs}_i(n\Delta t), \,\,\, i=1,2, \cdots, N_f, \,\,\, n=1, 2, \cdots N 
     173\end{equation} 
     174Based on a variational assimilation approach (REFS), where we  
     175match observations   
     176to a model $\mathbf{r}(t)=(x(t),y(t))=\mathcal{M} (\bo{u}, \bo{r}$ that simulates drifters' trajectories.  
     177The proposed improvements are: 
     178 \begin{enumerate} 
     179 \item Accounting for the wind: we shall consider that $\bo{u}=\bo{u_{geo}}+\bo{u_{wind}}$  
     180 \item Physical constraint of the geostrophic part to be divergence free. 
     181 using a first guess or background made of AVISO and wind fields. details in section. 
     182different is: all at once, ,moving window, time varying, wind and divergence. added to the objective function 
     183\item smoothed correction with the moving window. 
     184 \end{enumerate} 
     185 
     186 
     187 
     188\subsection{Model for Lagrangian data} 
     189 
     190 
     191The position of a specific drifter is the solution of the non-linear advection equation 
     192\begin{equation} \label{advection} 
     193 \odif{\mathbf{r}}{t}=\mathbf{u}(\mathbf{r}(t),t), \,\,\,\,\bo{r}(0)=\bo{r}_0, \mathbf{u}(x,y,0)=\bo{u}_0. 
     194 \end{equation} 
     195This equation \eqref{advection} is integrated numerically using an Euler scheme for example. Also, since the drifters may not coincide the Eulerian velocity's grid, a spatial interpolation of $\mathbf{u}$ to these positions is needed. The numerical model is therefore: 
     196 
     197\begin{equation} \label{euler_advection} 
     198\bo{r}(k\delta t)=\bo{r}((k-1)\delta t)+2\delta t \, interp(\bo{u}, \bo{r}((k-1)\delta t)), \,\,\,\,\, k=1,2, \cdots 
     199\end{equation} 
     200%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 
     201where $\delta t$ the time step of the scheme, typically a fraction of $\Delta t$. We can choose bilinear interpolation  
     202\begin{align} 
     203interp(\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 
     204&+ (\mathbf{u}_1-\mathbf{u}_2-\mathbf{u}_3 + \mathbf{u}_4 )\frac{(x-x_1 )(y -y_1 )}{\Delta x \, \Delta y}  
     205\end{align} 
     206 
     207\begin{align} \notag 
     208\mathbf{u}_1 &= \mathbf{u} (x_1 , y_1 ), \\  \notag 
     209\mathbf{u}_2 &= \mathbf{u} (x_1 + 1, y_1 ), \\ \notag 
     210\mathbf{u}_3 &=\mathbf{u}(x_1 , y_1 + 1), \\  \notag 
     211 \mathbf{u}_4 &= \mathbf{u}(x_1 + 1, y_1 + 1). \notag 
     212\end{align} 
     213 
     214 
     215\subsection{The Linear tangent model} 
     216%We write $\mathbf{r}=\mathbf{r}^b+\delta \bo{r}$ and $\mathbf{u}=\mathbf{u}^b+\delta \mathbf{u}$ and obtain 
     217%\begin{align} \label{Evolve} 
     218%\bo{r}^b_k&=\bo{r}^b_{k-2}+2\delta t \, interp(\mathbf{u}^b, \bo{r}^b_{k-1})  \\  
     219%\delta \bo{r}_k &= \delta \bo{r}_{k-2} + 2 \delta t \, \{ interp(\delta \bo{u},\bo{r}^b_{k-1})+ \delta \bo{r}_{k-1} \cdot \partial _{(x,y)} interp (\bo{u}^b,\bo{r}^b_{k-1})\} \notag 
     220%\end{align} 
     221%where $  \partial _{(x,y)} interp$ is the derivative of the ÔinterpÕ function with respect to $(x,y)$ (not differentiable at $(x_1,y_1)$).  
     222%We initialize this with $\bo{r}^b_{0,1}=\bo{r}^{obs}_{0,1}$ and $\delta \bo{r}_{0,1}=\bo{0}.$ 
     223% 
     224% 
     225%\textbf{Remark}: Technically here, $\delta \bo{u}$ and $\bo{u}^b$ are time independent.  
     226%If we were to include the full physical model (denote it by $\mathcal{M}$), we would have \[\mathcal{M}(t_i,t_0)(\bo{u_0})=\mathcal{M}(t_i,t_0)(\bo{u^b+\delta \bo{u}})=\bo{u}_i\approx \mathcal{M}(t_i,t_0)(\bo{u}^b)+\mathcal{M}'(t_i,t_o) \delta \bo{u},\]   that is we would have time dependent $\bo{u}^b$ and  $\delta \bo{u}$. 
     227% 
     228%\subsection{Discussion on the choice of $\delta t$ } 
     229%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:  
     230%\begin{enumerate}[(i)] 
     231%\item If $\Delta s >> \Delta x$, then we can take $\delta t \sim \Delta t$ (that should simplify equation \eqref{Evolve})  
     232%\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$ 
     233%\end{enumerate} 
     234% 
     235% 
     236%\subsection{Transformation} 
     237%We can transform the state vector using $\delta \bo{v}=\bo{B}^{-1/2}\delta \bo{u}$. We can define this matrix directly using smoothing operators \cite{weaver2006correlation}. This formulation has the advantage that $\delta \bo{v}$ has now the identity matrix as a covariance matrix.  
     238% 
     239% 
     240%\subsection{SLIDES} 
     241%Minimize objective function with respect to the state $\vec{U_0}$  
     242%\[ 
     243%\mathcal{J}(\vec{U_0})= \frac{1}{2} 
     244% \left \{ 
     245% \sum_{i,m} 
     246% \vectornorm{ \mathcal{G} \mathcal{M}(\vec{U_0})-\vec{r}_i^{\,o}(m\Delta t) } ^2+ \vectornorm{ \vec{U_0}-\vec{U}_b}^2_{\mathbf{B}} 
     247% \right \}\] 
     248%\vspace{2mm} 
     249%where $\vectornorm{\psi}^2_{\bo{B}} \equiv \psi^T \mathbf{B}^{-1} \psi.$ 
     250% 
     251%\begin{itemize} 
     252%\item $\vec{U}^b$ is background state of  $\vec{U_0}$, with error covariance matrix $\mathbf{B}^{-1}$ 
     253%\small{ 
     254%\begin{itemize} 
     255%\item  dual purpose: regularization and information spreading (smoothing) 
     256%%\item in extremely idealized assimilation contexts, for example with floaters "everywhere", can use identity 
     257%\item use the diffusion filter method [Weaver and Courtier, 2001] 
     258%%need R  
     259%\end{itemize} 
     260%} 
     261%\item If $\mathcal{J}$ is not quadratic in $\vec{U}_0$, optimization is very costly and may not converge: work with linear tangents (around the background) of model and advection  
     262%\item  Incremental approach takes into account weak nonlinearities, as models are updated several times. 
     263%\end{itemize} 
     264% 
     265% 
     266%%%%%%--- 
     267% 
     268% 
     269% 
     270 
     271\subsection{Optimization} 
     272 
    135273\begin{itemize}  
    136274 
    … …  
    139277\small{ 
    140278\begin{align} \notag 
    141 &\vec{r}^{\,b}_i(k\delta t)=\vec{r}^{\,b}_i((k-1)\delta t)+2\delta t \, interp(\vec{U}^b, \vec{r}^{\,b}_i((k-1)\delta t)),\,\,\,\,\, \text{\textcolor{blue}{model}}  \\  \notag 
    142 &\vec{\delta r}_i(k\delta t) = \vec{\delta r}_i((k-1)\delta t) + 2 \delta t \, \{ interp(\textcolor{green}{\vec{\delta U}},\vec{r}^{\,b}_i((k-1)\delta t)) \\ \notag  
    143 &+ \vec{\delta {r}}_i((k-1)\delta t) \cdot \partial _{(x,y)} interp (\vec{U}^b,\vec{r}^{\,b}_i((k-1)\delta t))\}  \,\,\, \text{\textcolor{blue}{tangent}} \notag 
     279&\vec{r}^{\,b}_i(k\delta t)=\vec{r}^{\,b}_i((k-1)\delta t)+2\delta t \, interp(\vec{U}^b, \vec{r}^{\,b}_i((k-1)\delta t)),\,\,\,\,\, \text{model}  \\  \notag 
     280&\vec{\delta r}_i(k\delta t) = \vec{\delta r}_i((k-1)\delta t) + 2 \delta t \, \{ interp(\vec{\delta U},\vec{r}^{\,b}_i((k-1)\delta t)) \\ \notag  
     281&+ \vec{\delta {r}}_i((k-1)\delta t) \cdot \partial _{(x,y)} interp (\vec{U}^b,\vec{r}^{\,b}_i((k-1)\delta t))\}  \,\,\, \text{tangent} \notag 
    144282\end{align} 
    145283initialized with observations, for $k=1,2,3, \cdots $ and $i=1, \cdots N_{drift}$  
    … …  
    147285%\vec{r}(0)= \vec{r}^{o}(0)$ 
    148286\item Perform sequences of optimization 
    149 \[\mathcal{J}(\textcolor{green}{\delta \vec{U}})=  \sum _i \sum_{m=1}^{P} \vectornorm{\vec{r}^{\,b}_i+\delta \vec{r}_i(\textcolor{green}{\delta \vec{U}}) -\vec{r}_i^{\,o}(m\Delta t) }^2 +\alpha_1 \,\vectornorm{ \textcolor{green}{\vec{\delta U}} }^2_{\bo{B}} +\alpha_2 \,div (\vec{U}_{geo}) \] 
     287\[\mathcal{J}(\delta \vec{U})=  \sum _i \sum_{m=1}^{P} \vectornorm{\vec{r}^{\,b}_i+\delta \vec{r}_i(\delta \vec{U}) -\vec{r}_i^{\,o}(m\Delta t) }^2 +\alpha_1 \,\vectornorm{ \vec{\delta U} }^2_{\bo{B}} +\alpha_2 \,div (\vec{U}_{geo}) \] 
    150288$P$ depends on $T_{L}/\Delta t$. $T_L$ characteristic of autocorrelation of drifters (typical 1-3 days)  
    151289 
    … …  
    156294 
    157295 
     296 
     297 
     298 
     299 
     300 
    158301\section{Twin Experiments} 
     302In the twin experiment approach, the observations are simulated using a dynamic model for the Eulerian velocity field. \textcolor{red}{Need details}. 
     303 
     304 
     305\begin{figure}[htbp] 
     306\begin{center} 
     307\includegraphics[scale=0.5]{./fig/RegionErroronAviso.pdf} 
     308\vspace{-30mm} 
     309\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, } 
     310\end{center} 
     311\end{figure} 
     312 
     313 
     314 
     315\begin{figure}[htbp] 
     316\begin{center} 
     317\includegraphics[scale=0.5]{./fig/Shifts_win24_dt1_f14_tf72.pdf} 
     318\vspace{-30mm} 
     319\caption{The effect of the moving window for a window of size $24$hs with 14 drifters and $dt=2hs$ } 
     320\end{center} 
     321\end{figure} 
     322 
     323 
     324\begin{figure}[htbp] 
     325\begin{center} 
     326\includegraphics[scale=0.5]{./fig/Wins_optshift_dt1_f14_tf72.pdf} 
     327\vspace{-30mm} 
     328\caption{The effect of the window size for 14 drifters and $dt=2$hs} 
     329\end{center} 
     330\end{figure} 
     331 
     332\begin{figure}[htbp] 
     333\begin{center} 
     334\includegraphics[scale=0.5]{./fig/Nfs_win24_dt1_tf72.pdf} 
     335\vspace{-30mm} 
     336\caption{The effect of the number of drifters} 
     337\end{center} 
     338\end{figure} 
     339 
     340\begin{figure}[htbp] 
     341\begin{center} 
     342\includegraphics[scale=0.5]{./fig/Dts_win24_f14_tf72.pdf} 
     343\vspace{-30mm} 
     344\caption{The effect of the time sampling} 
     345\end{center} 
     346\end{figure} 
     347 
     348 
     349 
     350 
     351 
     352\begin{figure}[htdp] 
     353%\centering 
     354\begin{subfigure}{0.55\textwidth} 
     355 % \entering 
     356  \includegraphics[width=1.2\linewidth]{./fig/Before_L2pointw_win24_MEAN_color.pdf} 
     357%  \caption{Pressure drop.} 
     358  \label{bla} 
     359\end{subfigure}% 
     360\hspace{-10mm} 
     361\begin{subfigure}{0.55\textwidth} 
     362 % \centering 
     363  \includegraphics[width=1.2\linewidth]{./fig/After_L2pointw_win24_MEAN_color.pdf} 
     364%  \caption{Mass flowrate.} 
     365  \label{blo} 
     366\end{subfigure} 
     367\vspace{-30mm} 
     368\caption{the point-wise L2 error before and after in m/s} 
     369\label{fblo} 
     370\end{figure} 
     371 
     372 
     373\begin{figure}[htdp] 
     374%\centering 
     375\begin{subfigure}{0.55\textwidth} 
     376 % \entering 
     377  \includegraphics[width=1.2\linewidth]{./fig/Before_CORRANGLEpointw_win24_MEAN_color.pdf} 
     378%  \caption{Pressure drop.} 
     379  \label{bla} 
     380\end{subfigure}% 
     381\hspace{-10mm} 
     382\begin{subfigure}{0.55\textwidth} 
     383 % \centering 
     384  \includegraphics[width=1.2\linewidth]{./fig/After_CORRANGLEpointw_win24_MEAN_color.pdf} 
     385%  \caption{Mass flowrate.} 
     386  \label{blo} 
     387\end{subfigure} 
     388\vspace{-30mm} 
     389\caption{the point angle correlation error before and after in m/s} 
     390\label{fblo} 
     391\end{figure} 
     392 
     393 
     394 
     395 
    159396\section{Experiment with Real Data} 
    160 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. 
    161 \subsection{Improvement of velocity field near the coast} 
    162 %lebanese drifters 
    163  
    164 \subsection{Improvement of velocity field in an eddy} 
    165 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. 
    166  
    167 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. 
    168 \begin{figure}[h] 
    169 \centering 
    170 \includegraphics[scale =0.6]{./fig/Eddy_velocity.png} 
    171 \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.} 
    172 \end{figure} 
     397 
     398 
     399 
     400 
     401 
     402\begin{figure}[htbp] 
     403\begin{center} 
     404\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_6days_average.pdf} 
     405\vspace{-30mm} 
     406\caption{Real Exp starting Aug 28 for 6 days, window=24h, move=6hs, } 
     407\end{center} 
     408\end{figure} 
     409 
     410 
     411\begin{figure}[htbp] 
     412\begin{center} 
     413\includegraphics[scale=0.5]{./fig/ReconstructedCNRSExp_bmtogreen_2days_average_zoom.pdf} 
     414\vspace{-30mm} 
     415\caption{Real 2-1 Exp starting Aug 28 for 2 days, window=24h, move=24hs, } 
     416\end{center} 
     417\end{figure} 
     418 
     419 
    173420 
    174421%% The Appendices part is started with the command \appendix; 
Note: See TracChangeset for help on using the changeset viewer.