RIA
laboria

laboratoire de recherche
en informatique
et automatique

On the generation
of currents by winds:
an identification method
to determine oceanic parameters.

<figure><img src="image_1.png" /></figure>

D. Begis, M. Crepon

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

Rapport de Recherche n° 118
Mai 1975
ON THE GENERATION OF CURRENTS BY WINDS :
AN IDENTIFICATION METHOD TO DETERMINE OCEANIC PARAMETERS

D. Begis, M. Crepon

RESUME :

Lorsqu'il y a une thermocline bien connue et si l'observation est faite assez loin de la côte, nous observons que le courant dans la couche de surface dépend principalement du vent local. Il y a une formulation déterministe entre le courant et le vent, qui est une équation différentielle dont le modèle est connu mais dont les coefficients sont à déterminer. Par une méthode de gradient, on calcule ces coefficients par la minimisation d'une fonctionnelle quadratique. Ces coefficients sont le coefficient de frottement $C_{10}$, le paramètre de Coriolis $f$, et le coefficient d'amortissement. Ceci peut conduire à la construction d'une méthode de prévision de courants océaniques.

ABSTRACT

When there is a well defined thermocline and if the observation point is far enough from the shore, one shows that the current in the mixed layer depends mainly on the local wind. There is a deterministic relationship between current and wind which is of the form of a differential equation whose structure is known, but the coefficients unknown. From a gradient method, one calculates these coefficients by minimizing a quadratic functional. They are the drag coefficient $C_{10}$, the Coriolis parameter $f$ and the damping law. This leads to the possibility of building an oceanic current forecast.
On the generation of currents by winds ;

an identification method to determine oceanic parameters.

D. BEGIS* and M. CREPON **

INTRODUCTION

As far as we know the only attempt to compare observed and wind simulated currents was made by R. POLLARD and R. MILLARD (1970). In the present paper we deal with models which do not depend on spatial variations of winds and which are extensions of POLLARD and MILLARD'S model. A fundamental question was still unsolved : how can spatial variations of winds and the presence of a shore affect the validity of these calculations ?

In the first part we give an analytical justification of the foregoing assumption using a geostrophic adjustment process. We find that its validity depends on the ratio of the distance of the observation point from the shore over the internal radius of deformation.

<footnote>* IRIA, Laboratoire de Recherche en Informatique et Automatique, Rocquencourt, 78150, France.</footnote>

<footnote>** Laboratoire d'Oceanographie physique du Museum National d'Histoire Naturelle PARIS 5ème</footnote>
If this ratio is far greater than one, the model is valid.

In the second part, we check different models. First we deal with a model similar to that of POLLARD and MILLARD but our aim is different. We identify the coefficients of the differential equation which connects the wind and the current together by minimizing a quadratic functional which represents the difference between observed and wind simulated currents. Hence we get the best values for the drag coefficient $C_{10}$, for the Coriolis parameter $f$ and for the damping law. We also apply the same technique to a depth dependent model with a no slip condition on the thermocline.

# First Part. Theoretical justification of models.

## 1.1. - Formulation of the problem

Let a two layer ocean be. We deal with the linearized equations of hydrodynamics. The pressure is supposed to be hydrostatic. We assume the coupling between the upper layer and the lower layer is only due to pressure forces. In other words, the upper layer slips without any stresses on the lower layer. This model is well-known and schematizes very well the true ocean when there is a well defined thermocline and for meso-scale motions ; see for example G. VERONIS and H. STOMMEL (1956), M. RATTRAY (1964), J. GONELLA, M. CREPON and F. MADELAIN (1969), J. GONELLA (1971), R. POLLARD, P. RHINES and R. THOMPSON (1973).

We have the following equations in the upper layer.

$$ (1.1.) \quad \begin{cases} \left(\frac{\partial}{\partial t} + v_1\right) \vec{u}_1 + f \vec{z} \times \vec{u}_1 = -h_1 \nabla \phi_1 + \frac{\vec{\tau}}{\rho_1} \\ \\ \frac{\partial}{\partial t} (\zeta_1 - \zeta_2) + \nabla \vec{u}_1 = 0 \end{cases} $$
with $\varphi_1 = g \zeta_1$

and in the lower layer

$$ (1.2.) \quad \left\{ \begin{array}{l} \left(\frac{\partial}{\partial t} + v_2\right) \vec{u}_2 + f \vec{z} \times \vec{u}_2 = -h_2 \nabla \varphi_2 \\ \\ \frac{\partial}{\partial t} \zeta_2 + \nabla \vec{u}_2 = 0 \end{array} \right. $$

with $\varphi_2 = g \frac{\rho_1}{\rho_2} \zeta_1 + g \frac{\Delta \rho}{\rho_2} \zeta_2$

where $\vec{u}_i(x, y, t)$ is the integrated velocity from the bottom to the interface (Suffix 2) or from the interface to the surface (Suffix 1).

$\zeta_i(x, y, t)$ is the elevation of the sea surface (Suffix 1) or of the interface (Suffix 2) from the rest.

$f$ is the Coriolis parameter assumed to be constant.

$\vec{z}$ the unit vector along the vertical, positive upwards.

$h_i(x, y)$ the depth of the upper layer (Suffix 1) or of the lower layer (Suffix 2) at rest.

$\rho_i$ the density of layer i.

$\vec{\tau}(x, y, t)$ the wind stress.

$\nabla$ the bi-dimensional operator $\frac{\partial}{\partial x}$, $\frac{\partial}{\partial y}$

$v_i$ a dissipation coefficient in layer i.
The first equation (1.1.) yields :

$$ (1.3.) \quad \left[ \left( \frac{\partial}{\partial t} + v_i \right)^2 + f^2 \right] \vec{u}_1 = \frac{1}{\rho_1} \left( \vec{\tau} - h_1 \rho_1 \nabla \varphi_1 \right) \left[ \left( \frac{\partial}{\partial t} + v_i \right) + \times \vec{z} \, f \right] $$

The integrated velocity in the upper layer $\vec{u}_1$ depends both on the local wind $\vec{\tau}$ and on the spatial structure of the wind through $\nabla\varphi_1$. If $h_1\rho_1 \nabla\varphi_1$ is small with respect to $\vec{\tau}$, the integrated velocity is a function of the local wind only. Therefore the main problem is to determine $\varphi_1$. The resolution of the time dependent system formed by (1.1) and (1.2.) is difficult. We solve it by the Laplace Transform with respect to the time.

We assume the motions start from the rest, ie $\vec{u}_i, \zeta_i, \vec{\tau}$ are equal to zero for $t < 0$. Let $\vec{U}_i, Z_i, \Phi_i, \vec{T}$ be the Laplace Transform respectively of $\vec{u}_i, \zeta_i, \varphi_i, \vec{\tau}$ such as.:

$$ \vec{U}_i(x, y, p) = \int_0^\infty e^{-pt} \vec{u}_i(x, y, t) dt $$

From (1.1) and (1.2.) we obtain the following equation for $\Phi_1$ and $\Phi_2$:

$$ (1.4.) \quad \left\{ \begin{array}{l} \Phi_1 - \Phi_2 - g \frac{\Delta \rho}{\rho} \left[ \frac{h_1}{p^2 k_1^2} \nabla^2 \Phi_1 - \frac{1}{\rho_1 p^2 k_1^2} (\nabla \vec{T} + \frac{f \vec{z}}{p+v_1} \nabla \times \vec{T}) \right] = 0 \\ \Phi_2 - \frac{\rho_1}{\rho_2} \Phi_1 - g \frac{\Delta \rho}{\rho} \frac{h_2}{p^2 k_2^2} \nabla^2 \Phi_2 = 0 \end{array} \right. $$

$$ \text{with } k_i = \left[ (p + v_i)^2 + f^2 \right] / \left[ p (p + v_i) \right] $$

In order to simplify the calculations, we assume $v_1 = v_2$ in the next part.

From (1.4.) and the boundary conditions we can obtain $\Phi_1$ and $\Phi_2$.

Thus, the inverse of Laplace Transform yields $\varphi_1$.
<figure><img src="image_2.png" /><figcaption>Fig. 1 - Typical geometry for the two-layers model with a free surface and an interface.</figcaption></figure>
Let us now solve a simple case in order to understand the
mechanism of the problem

1.2. -Constant wind in space blowing over a semi-infinite ocean

Let a semi-infinite ocean be bounded by a shore of equation
x = 0 (fig. 1)

Let us consider a wind stress of the form.

$$
\vec{\tau} = (\tau_x \vec{x} + \tau_y \vec{y}) Y(x) Y(t),
$$

where $Y(x)$ is the step function

$\tau_x(t)$ (resp. $\tau_y(t)$) the component of the wind stress on
$x$ axis (resp. $y$ axis), is a function of $t$ only. (Wind constant
in space).

We have

$$
\operatorname{div} \vec{\tau} = \nabla \cdot \vec{\tau} = \tau_x \delta(x)
$$

$$
\operatorname{curl} \vec{\tau} = \nabla \times \vec{\tau} = \tau_y \delta(x)
$$

δ ( ) is the Dirac function

As we can see on eq (1.4)the effects of the curl ' of the wind
stress are more important than those of its divergence, for mesoscale
motions (if the duration of the meteorological perturbation is greater
than the Coriolis period 2π/f). Thus, we shall only deal in the following
part with a wind of the form :

$$
(1.5.) \quad \vec{\tau} = \tau_y(t) \vec{y} Y(x) Y(t)
$$

which wind is parallel to the shore.

The equations (1.3.) and (1.4.) yield for the integrated
velocity in the upper layer :

$$
\begin{equation}
\begin{split}
(1.6.) \quad \vec{u}_1(x,t) = {}& \vec{u}_L Y(t) + \frac{h_1}{h_1+h_2} \vec{u} (x/C_1, t) Y(t-x/C_1) + \\
& \qquad \frac{h_2}{h_1+h_2} \vec{u} (x/C_2, t) Y(t-x/C_2)
\end{split}
\end{equation}
$$
$\vec{u}_L(t)Y(t)$ is the integrated velocity due to the local wind.

$\vec{u}(x/C_i, t)Y(t-x/C_i)$ is an integrated velocity due to the vertical boundary. This integrated velocity is generated at the shore and its front propagates from there with a velocity $C_i$ where $C_i$ is the velocity of long barotropic waves, $C_2$ the velocity of long baroclinic waves.

$$ (C_1 \approx \sqrt{g(h_1+h_2)} , C_2 \approx \sqrt{\frac{h_1 h_2}{g(h_1+h_2)} \frac{\Delta \rho}{\rho}} $$

For a wind stress of the form (1.5), we obtain, if we neglect the damping coefficient ($\nu=0$).

$$ (1.7.) \qquad \begin{aligned} \vec{u}_L(t)Y(t) &= \frac{1}{\rho f} \tau_y(t)Y(t) \quad {}^* \text{(ft)} \left[ \sin(ft) \vec{x} + \cos(ft) \vec{y} \right] \\ \vec{u}(x/C_i, t)Y(t-x/C_i) &= \frac{1}{\rho f} \tau_y(t)Y(t) \quad {}^* \text{(ft)} \left[ -\sin(f(t-x/C_i))Y(t-x/C_i)\vec{x} \right. \\ & \qquad \left. + (1-\cos(f(t-x/C_i)))Y(t-x/C_i)\vec{y} + \frac{\frac{J_1(f\sqrt{t^2-x^2/C_i^2})}{C_i}}{f\sqrt{t^2-x^2/C_i^2}} Y(t-x/C_i) \quad {}^* \text{(ft)} \right. \\ & \qquad \left. \sin(ft)Y(t)\vec{x} + (1-\cos(ft))Y(t)\vec{y} \right] \end{aligned} $$

where the symbol $*$ represents the convolution product, $J_n$ is the Bessel function of the first kind of order n.

From the properties of the Laplace Transform, we can deduce the asymptotic behaviour of $\vec{u}(x/C_i, t)$ as $t$ tends towards infinity by applying the method developed by W. SUTTON (1934). If we suppose $\tau_y(t) = \tau_o Y(t)$, this yields:
$$ (1.8.) \quad \vec{u}(x/C_i, t) = \frac{\tau_0}{\rho f} \left[ \begin{aligned} & \left(e^{-\frac{xf}{C_i}} + \cos(\frac{\pi}{2}t)\right) \vec{x} + \\ & \left(e^{-\frac{xf}{C_i}} (\text{ft}) - \sin(\frac{\pi}{2}t)\right) \vec{y} + O\left(\frac{\frac{xf}{C_i}}{(\text{ft})^{1/2}}\right) \end{aligned} \right] $$

The integrated velocity increases as $e^{-\frac{xf}{C_i}}$ (ft) as (ft) increases and it decreases as we go farther from the shore. There is a generation of a costal jet stream (G. CSANADY - 1972, J.O.BRIEN and H. HURLBURT 1972).

For deep ocean, the most important term of (1.6.) is the first one. The second one which is of order $O(\frac{h_1}{h_1+h_2}) = 10^{-2}$, can be neglected. The third one is equal to zero as long as $t < x/C_2$, and thereafter increases very slowly. The greater the ratio $xf/C_2$ is, the more slowly $\vec{u}(x/C_2, t)Y(t-x/C_2)$ increases, as we see on (1.8.).

A complete calculation of eq (1.7.) shows that $\vec{u}(x/C_2, t)Y(t-x/C_2) \ll \vec{u}_L(t)Y(t)$ as long as $ft < m \frac{xf}{C_2}$, $m$ being a number of order of some units and an increasing function of $\frac{xf}{C_2}$.

We conclude that the current will depend on local wind if.

1) the ratio $xf/C_2 \gg 1$, ie if the ratio of the distance $x$ between the observation point and the shore over the internal radius of deformation $C_2/f$ is far greater than one.

2) the duration of the storm is not too long with respect to the Coriolis period.

**Remark 1.1. :**

A long internal wave is associated to the mass transport $\vec{u}(x/C_2, t)Y(t-x/C_2)$. It is of the form

$$ (1.9.) \quad \zeta_2(x, t) = - \frac{1}{\rho c f} \sqrt{\frac{h_2 p}{h_1 \Delta \rho}} \tau_y(t) Y(t) \quad (\text{ft})^* $$

$$ \left[ J_o f \sqrt{t^2 - x^2/C_2^2} Y(t-x/C_2) \quad (\text{ft})^* Y(t) \right] $$
<figure><img src="image_3.png" /><figcaption>Fig. 2 - Position of the "Bouée-Laboratoire".</figcaption></figure>
Remark 1.2. :

The equations (1.7.) and (1.9.) are similar to these obtained
by A. CAHN (1945) for solving the time dependent geostrophic adjustment
problem. In fact, we have solved a baroclinic geostrophic adjustment
problem.

Second part. Identification method : theory and results

During the summer 1969, currents and wind measurements were done on
the French "bouée Lab". By that time, the "bouée Lab" was anchored between
France and Corsica (fig. 2), at a distance d of about 100 km off the shore.
The depth of the thermocline h₁ was about 25 m., the depth of the sea
h₂ 2500 m., and the difference of density Δρ/ρ above and below the thermo-
cline of order 10⁻³. Thus, the internal radius of deformation C₂/f = 0(5 km)
and the ratio d/C₂/f = 0 (20). From the results of the first part, we can
state the conditions of validity of a model which depends only on local wind
are fully filled.

We first deal with the following model, which is a slab-like model.

$$
(2.1) \quad \left\{
\begin{aligned}
& \left[ \frac{\partial}{\partial t} + v_{\alpha} (u^2 + v^2)^{\alpha} \right] u -fv = C_{10} \frac{\rho_a}{\rho_w} |w| w_x \\
& \left[ \frac{\partial}{\partial t} + v_{\alpha} (u^2 + v^2)^{\alpha} \right] v + fu = C_{10} \frac{\rho_a}{\rho_w} |w| w_y
\end{aligned}
\right.
$$

$$
\text{with } \vec{u} = \int_{-h}^{0} \vec{u} * dz
$$

where $\vec{u} *$ is the velocity at depth z.

h is the depth of the thermocline.

$\vec{W}(t)=(w_x, w_y)$ is the wind vector, $|w|$ the wind modulus.
The term $v_{\alpha}(u^2+v^2)^{\alpha}$ takes into account the loss of energy of the mean flow into turbulence and dissipation. The case $\alpha = 0$ corresponds to a damping proportional to the mass transport (POLLARD and MILLARD'S model). The case $\alpha = \frac{1}{2}$ corresponds to a quadratic damping (MANNING formula, cf C. BRETSCHNEIDER (1966)).

The unknowns of the problem are $C_{10}$, $f$, and $v_{\alpha}$. We solve it by a gradient method.

## 2.1. - Statement of the problem

Let us denote :

$$ -\mathcal{U}_{\emptyset}(t) = (u_{\emptyset}(t), v_{\emptyset}(t)) \text{ the observed current on the period of time } ]0,T[ $$

$$ - \mathcal{V} = (v, f, C_{10}) \text{ the control vector.} $$

$$ - \mathcal{V}_{\text{ad}} \text{ the set of the boundary constraints on } \mathcal{V}. $$

$$ - \mathcal{U} = (u(t), v(t)) \text{ the calculated current which is the solution of the equation (2.1). This equation can be generalized by (2.1')} $$

$$ (2.1') \quad \begin{cases} \left[ \frac{d}{dt} + v\phi(u^2+v^2) u(t) \right] - fv(t) = \frac{C_{10}}{\rho} |w| w_x(t) ; u(0) = u_{\emptyset}(0) \\ \left[ \frac{d}{dt} + v\phi(u^2+v^2) v(t) \right] + fu(t) = \frac{C_{10}}{\rho} |w| w_y(t) ; v(0) = v_{\emptyset}(0) \end{cases} $$

where $\phi(x)$ is a positive, differentiable and increasing function.

$$ - C_{v,f,C_{10}}(\mathcal{U}) = (C) \text{ the cost functional defined by} $$

$$ (2.2.) \quad C_{v,f,C_{10}}(\mathcal{U}) = \int_{0}^{T} \{ (u_{\emptyset}(t) - u(t))^2 + (v_{\emptyset}(t) - v(t))^2 \} dt. $$
We consider the following identification problem :
we compute in $V_{ad}$ the best estimate of $(v, f, C_{10})$ so that
$C_{v,f,C_{10}}(u)$ is minimum. Hence we calculate the gradient of
$C(u)$ with respect to $v, f$ and $C_{10}$.

2.2. - Necessary conditions of optimality

We cannot directly get the gradient in V of C. For this we
need an intermediate function (adjoint state) and we compute
the gradient as a function of the adjoint state and of the
state **u**.

We use the optimal control theory (see LIONS).

Let us call

$$
- \mathcal{P} = (p(t), q(t)) \text{ the adjoint state}
$$

$$
- \mathcal{L}(u, \mathcal{P}, v) \text{ the hamiltonian defined by :}
$$

$$
\begin{align*}
(2.3.) \quad \mathcal{L}(u, \mathcal{P}, v) &= \int_0^T (u(t) - u_\phi(t))^2 dt + \int_0^T (v(t) - v_\phi(t))^2 dt \\
&\quad - 2 \int_0^T p(t) (u'(t) + v\phi(u^2 + v^2) u(t) - fv(t) - C_{10} |\frac{w}{\rho}| w_x(t)) dt \\
&\quad - 2 \int_0^T q(t) (v'(t) + v\phi(u^2 + v^2) v(t) + fu(t) - C_{10} |\frac{w}{\rho}| w_y(t)) dt
\end{align*}
$$

the optimality in **u** gives the equation of **P**

$$
\left(\frac{\partial \mathcal{L}}{\partial u} = 0\right)
$$

$$
\begin{equation}
\begin{aligned}
& \qquad - p'(t) + v_p(t) (\phi(u^2 + v^2) + 2u^2 \phi'(u^2 + v^2)) + (2 v u(t) v(t) \phi'(u^2 + v^2) + f) q(t) \\
& \qquad - q'(t) + v_q(t) (\phi(u^2 + v^2) + 2v^2 \phi'(u^2 + v^2)) \\
& \qquad \qquad + (2 v u(t) v(t) \phi'(u^2 + v^2) - f) p(t) \\
& = u(t) - u_{\phi}(t) \\
& \qquad - v(t) - v_{\phi}(t)
\end{aligned}
\tag{2.4.}
\end{equation}
$$
with $p(t) = q(t) = 0$

the optimality in V :

$$ \{- \int_{0}^{T} \varphi (u^2 + v^2) (p(t)u(t) + q(t)v(t))dt\} \quad (\mu - \nu) \geq 0 $$

$$ \begin{array}{l} \left\{ - \int_{0}^{T} (q(t)u(t) - p(t)v(t))dt \right\} (g-f) \geq 0 \quad \forall (\mu,g,d) \in \text{Vad} \\ \text{and } (v,f,C_{10}) \in \text{Vad} \end{array} $$

$$ \left\{ \int_{0}^{T} (p(t)w_{x}(t) + q(t)w_{y}(t))dt \right\} (d - C_{10}) > 0 $$

## 2.3. - Approximation

Let $k$ be the time step $-k = T/N$ where $N \in \mathbb{N}^* = \{k|k = 1,2,3,4,...\}$

We suppose that we have one observation every $k/I$ minutes with $I \in \mathbb{N}^*$

### 2.3.1. - Approximation of the cost functional.

Let us write

$$ u_i = u((i-1)k), \quad u_{\emptyset}^i = u_{\emptyset}((i-1)k), \quad \tau_x^i = |w| w_x((i-1)k), \quad \tau_y^i = |w| w_y((i-1)k) $$

$C_{v,f,C_{10}}(u)$ will be approximated by

$$ (2.5.) \qquad C_h(u) = k \sum_{i=1}^{N+1} \alpha_i \{ (u_{\emptyset}^i - u_i)^2 + (v_{\emptyset}^i - v_i)^2 \} $$

where $\alpha_i = 1$, except $\alpha_1 = \alpha_{N+1} = .5$

### 2.3.2. - Approximation of the state equation (2.1').

We have chosen a discretization scheme which is always stable.
Its accuracy is $O(k^2)$. The discretization of the term $\frac{du}{dt}$ is such as it is satisfied for the polynomials $t^\alpha$ where $\alpha \le 2$.

We can schematize (2.1) by $u' + Au = f$.

Knowing $u_{n-1}$, $u_n$, we compute $u_{n+1}$ with

$$ (2.6.) \qquad \frac{3/2 u_{n+1} - 2u_n^{1/2} u_{n-1}}{k} + Au_{n+1} = f_{n+1} $$

**Remark 2.1.** : There is a difficulty for starting the procedure because we only know $u_1 = u_\phi(0)$ - We compute $u_2$ with a two step scheme either the Crank-Nicholson scheme (accuracy $O(k^2)$) or an implicit one (accuracy $O(k)$) (see RICHTMEYER and MORTON).

We get the following system

$$ (2.7.) \qquad \left\{ \begin{aligned} & \frac{u_2 - u_1}{k} + v((1-\theta)\varphi(u_1^2+v_1^2)u_1 + \varphi(u_2^2+v_2^2)u_2) - (1-\theta)f v_1 - \epsilon f v_2 = C_{10}((1-\theta)\tau_x^1 + \theta\tau_x^2) \\ & \frac{v_2 - v_1}{k} + v((1-\theta)\varphi(u_1^2+v_1^2)v_1 + \varphi(u_2^2+v_2^2)v_2) + (1-\theta) f u_1 + \epsilon f u_2 = C_{10}((1-\theta)\tau_y^1 + \theta\tau_y^2) \\ & \frac{3/2 u_i - 2u_{i-1} + 1/2 u_{i-2}}{k} + v\varphi(u_i^2+v_i^2)u_i - f v_i = C_{10} \tau_x^i \\ & \frac{3/2 v_i - 2v_{i-1} + 1/2 v_{i-2}}{k} + v\varphi(u_i^2+v_i^2)v_i + f v_i = C_{10} \tau_y^i \end{aligned} \right. $$

for $i = 3$ to $N+1$.
$$ \text{and } \Theta = \frac{1}{2} \quad \text{(Crank Nicholson scheme)} $$

= 1 (implicit scheme)

we must solve a non linear system of the following type

$$ (2.8.) \quad \left\{ \begin{array}{l} \alpha u + \nu \varphi (u^2 + v^2) u - \nu v = b_1 \\ \alpha v + \nu \varphi (u^2 + v^2) v + \alpha u = b_2 \end{array} \right. $$

For solving such a system we have chosen two methods

1). The alternating direction iteration - ADI

2). Newton-Raphson iterative method.

1) ADI (see VARGA)

Let $\rho$ be positive, knowing $\mathcal{U}^m$, we compute $\mathcal{U}^{m+1}$ with

$$ (2.9.) \quad \left\{ \begin{array}{l} \rho(u^{m+1/2} - u^m) + \nu \varphi \left[ (u^{m+1/2})^2 + (v^{m+1/2})^2 \right] u^{m+1/2} + \alpha u^m - \nu v^m = b_1 \\ \rho(v^{m+1/2} - v^m) + \nu \varphi \left[ (u^{m+1/2})^2 + (v^{m+1/2})^2 \right] v^{m+1/2} + \alpha v^m + \alpha u^m = b_2 \end{array} \right. $$

$$ (2.10.) \quad \left\{ \begin{array}{l} \rho(u^{m+1} - u^m) + \nu \varphi \left[ (u^{m+1/2})^2 + (v^{m+1/2})^2 \right] u^{m+1/2} + \alpha u^{m+1} - \nu v^{m+1} = b_1 \\ \rho(v^{m+1} - v^m) + \nu \varphi \left[ (u^{m+1/2})^2 + (v^{m+1/2})^2 \right] v^{m+1/2} + \alpha v^{m+1} + \alpha u^{m+1} = b_2 \end{array} \right. $$

For solving (2.9.) we set $u^{m+1/2} = r \sin \theta$, $v^{m+1/2} = r \cos \theta$

$$ d_1 = b_1 + (\rho - \alpha) u^m + \nu v^m \quad \text{and} \quad d_2 = b_2 + (\rho - \alpha) v^m - \alpha u^m $$

we get

$$ (2.11.) \quad \nu r \varphi (r^2) + \rho r - \sqrt{d_1^2 + d_2^2} = 0 $$
If $\varphi(r^2)$ is equal to 1 or to r we have a simple equation

If $\varphi(r^2)$ is any suitable function we solve (2.11.) with
a Newton iterative method for one variable.

knowing r we get

$$
u^{n+1/2} = \frac{d_1}{\rho + \nu \varphi(r^2)} \quad \text{and} \quad v^{m+1/2} = \frac{d_2}{\rho + \nu \varphi(r^2)}
$$

(2.10.) being a linear system (2-2) we get easily
$u^{m+1}$ and $v^{m+1}$.

-We iterate until we get $\sup_i \{|u_i^{m+1} - u_i^m|, |v_i^{m+1} - v_i^m|\} \le \varepsilon$

And the solution of (2.8.) is $(u^{m+1}, v^{m+1})$

**Remark 2.2.** : ρ can be adjusted to improve the convergence
of the algorithm.

2) *Newton-Raphson method* (see DURAND)

Let $\Delta^S$ be the (2 - 2) - matrix

$$
\Delta^S =
\begin{pmatrix}
\alpha + \nu\varphi((u^S)^2 + (v^S)^2) + 2\nu(u^S)^2 \varphi'((u^S)^2 + (v^S)^2) & 2\nu u^S v^S \varphi'((u^S)^2 + (v^S)^2) - f \\
2\nu u^S v^S \varphi'((u^S)^2 + (v^S)^2) + f & \alpha + \nu\varphi((u^S)^2 + (v^S)^2) + 2\nu\varphi'((u^S)^2 + (v^S)^2)(v^S)^2
\end{pmatrix}
$$

we get the following algorithm

$$
(2.12.) \quad u^{s+1} = u^s - (\Delta^s)^{-1}
\begin{pmatrix}
\alpha u^s + \nu \varphi((u^s)^2 + (v^s)^2) u^s - f v^s - b_1 \\
\alpha v^s + \nu \varphi((u^s)^2 + (v^s)^2) v^s + f u^s - b_2
\end{pmatrix}
$$
We iterate until we get $\sup_i \{|u_i^{s+1} - u_i^s|, |v_i^{s+1} - v_i^s|\} \le \varepsilon$

**Remark 2.3.** : the advantage of this method is that an optimal parameter is not required. We only solve a sequence of (2-2) linear systems.

If $\varphi(x)$ is constant we obtain the linear system again and the method converges in one step

### 2.3.3. - Necessary conditions of optimality (we take $\theta = 1$)
The approximated Hamiltonian $\mathcal{L}_h$ of problem is defined by

$$ 
\begin{aligned}
\mathcal{L}_h(u,v,p,q,v,f, C_{10}) = C_h(\mathcal{U}) & - 2k \left\{ p_2 \frac{u_2^{-u_1}}{k} + v u_2 \varphi(u_2^2+v_2^2) -fv_2 - C_{10} \tau_x^2 \right\} \\
& + q_2 \left( \frac{v_2^{-v_1}}{k} + v v_2 \varphi(u_2^2+v_2^2) + f u_2 - C_{10} \tau_y^2 \right) \\
& + \sum_{i=3}^{N+1} \alpha_i \left( p_i \frac{3}{k} \frac{u_i^{-2u_{i-1}}}{u_{i-1}^{-1}} + \frac{1}{2} \frac{u_{i-2}}{u_i} \right. \\
& \qquad \left. + v u_i \varphi(u_i^2+v_i^2) - f v_i - C_{10} \tau_x^i \right) + q_i \left( \frac{3}{k} \frac{v_i^{-2v_{i-1}}}{v_{i-1}^{-1}} + \frac{1}{2} \frac{v_{i-2}}{v_i} \right) \\
& + v v_i \varphi(u_i^2+v_i^2) + f u_i - C_{10} \tau_y^i \right\}
\end{aligned}
$$

- $(u_i, v_i)$ is the solution of the non linear system (2.7.)

- $(p_i, q_i)$ is the solution of the linear system

$$ (2.14) \qquad \left( \frac{\partial \mathcal{L}_h}{\partial u_i}, \frac{\partial \mathcal{L}_h}{\partial v_i} \right) = 0 $$
- We have the optimality condition.

$$ (\frac{\partial \mathcal{L}}{\partial v})(U, \mathcal{P}, V)(W - v) \ge 0 \quad \forall W \in \text{Vad and } V \in \text{Vad} $$

i.e. i) optimality in $v$: $\left(\frac{\partial \mathcal{L}}{\partial v}\right)(\mu - v) \ge 0$

$$ \left(-2k \begin{array}{c} \Sigma \\ i=2 \end{array}\right) \alpha_i \varphi\left(u_i^2 + v_i^2\right)\left(p_i u_i + q_i v_i\right) (\mu - v) \geq 0 \quad \forall \mu $$

ii) optimality in $f$: $\left(\frac{\partial \mathcal{L}}{\partial f}\right)(g-f) \ge 0$

$$ \left(-2k \begin{array}{c} \Sigma \\ i=2 \end{array}\right) \alpha_i (q_i u_i - p_i v_i) (g-f) \geq 0 \quad \forall g $$

iii) optimality in $C_{10}^*$: $\left(\frac{\partial \mathcal{L}}{\partial c_{10}}\right)(d - c_{10}) \ge 0$

$$ \left(2k \begin{array}{c} \Sigma \\ i=2 \end{array}\right) \alpha_i (p_i \tau_x^i + q_i \tau_y^i) (d - c_{10}) \geq 0 \quad \forall d $$

where $(\mu, g, d) \in V_{ad}$

## 2.4. Gradient algorithms to solve the approximated problem

See GLOWINSKI - LIONS - TREMOLIERES

### 2.4.1. Gradient algorithm with projection.

$V^\circ$ given in $V_{ad}$ and $V^m$ being known we compute

i) $\mathcal{U}^{m+1}$ (by means of (2.7.))

ii) $\mathcal{P}^{m+1}$ (by (2.4.))

iii) $$ V^{m+1} = \operatorname*{Pr}_{V_{ad}} \left( V^m - \rho^m \frac{\partial \mathcal{L}}{\partial v} \left( \mathcal{U}^{m+1}, \mathcal{P}^{m+1}, V^m \right) \right) $$
Where $\overline{\text{Pr}}$ is the projection operator on $V_{ad}$ (there
$V_{ad}$
is no difficulty because $V_{ad}$ includes only boundary constraints) and $\rho^m$
is determined by dichotomy.

2.4.2. - Relaxation algorithm

It is obvious that $\mathcal{U}$ is a function of $V$ and $\mathcal{P}$ is a
function of $V$ and $\mathcal{U}$.

We shall denote $\mathcal{U}(v, f, c_{10})$ and $\mathcal{P}(\mathcal{U}, v, f, c_{10})$.

$V^\circ$ given on $V_{ad}$ and $V^m$ being known we compute :

- $\mathcal{U}(v^m, f^m, c_{10}^m)$ denoted $\mathcal{U}^m$.

- $\mathcal{P}(\mathcal{U}^m, v^m, f^m, c_{10}^m)$ denoted $\mathcal{P}^m$.

- $v^{m+1} = \text{Pr} \left( \nu^m - \rho^m \frac{\partial \mathcal{L}}{\partial \nu} h(\mathcal{U}^m, \mathcal{P}^m, \nu^m, f^m, c_{10}^m) \right)$.

- $\mathcal{U}(v^{m+1}, f^m, c_{10}^m)$ denoted $\mathcal{U}^{m+1/3}$.

- $\mathcal{P}(\mathcal{U}^{m+1/3}, v^{m+1}, f^m, c_{10}^m)$ denoted $\mathcal{P}^{m+1/3}$.

- $f^{m+1} = \text{Pr} \left( f^m - \rho^{m+1/3} \frac{\partial \mathcal{L}}{\partial f} h(\mathcal{U}^{m+1/3}, \mathcal{P}^{m+1/3}, v^{m+1}, f^m, c_{10}^m) \right)$.

- $\mathcal{U}(v^{m+1}, f^{m+1}, c_{10}^m)$ denoted $\mathcal{U}^{m+2/3}$.

- $\mathcal{P}(\mathcal{U}^{m+2/3}, v^{m+1}, f^{m+1}, c_{10}^m)$ denoted $\mathcal{P}^{m+2/3}$.

- $c_{10}^{m+1} = \text{Pr} \left( c_{10}^{m} - \rho^{m+1/2/3} \frac{\partial \mathcal{L}}{\partial c_{10}} h \left( \mathcal{U}^{m+2/3}, \mathcal{P}^{m+2/3}, v^{m+1}, f^{m+1}, c_{10}^{m} \right) \right)$.

Where $\rho^{m+i}/3$ is determined by dichotomy and

$V_{ad} = V_v \times V_f \times V_{c10}$

**Remark 2.4.** : At each step m+i/3 we minimize the
functional with respect to a single variable. We note that only the sign
of the derivative is required (and not its magnitude).
**Remark 2.5.** : For starting the procedure we can use the first algorithm.

**Remark 2.6.** : We can get accurate results by considering the initial value $\mathcal{U}(0)$ not as a data but as a parameter to be identified. We use the same algorithm, but now $\mathcal{U}^{m+1}$ depends on $\mathcal{U}^m(0)$ and by an additional gradient step we obtain.

$$ \mathcal{U}^{m+1}(0) = \text{Pr}(\mathcal{U}^m(0) - \rho \frac{\partial \mathcal{L}}{\partial \mathcal{U}_1} (\mathcal{U}^{m+1}, \mathcal{P}^{m+1}, \nu^{m+1}, f^{m+1}, c_{10}^{m+1})). $$

## 2.5. - Numerical results

We have applied the foregoing method to "COFRASOV" and "COBLAMED" data. Current measurements were done from the "bouée Lab" by Plessey current meters.

During "COFRASOV" experiment, there were records at 5m., 10 m., 15 m., 30 m. and 60 m.. The thermocline was at about 25 m.. The complex cross-correlation coefficient between data above the thermocline are high (over 0.8) and small for data on each side the thermocline (M. CREPON and D. CHOUCHAN), showing a slab like motions for the upper layer.

During "COBLAMED" experiment, there were records at 10 m., 25 m., and 30 m.. The thermocline was below 30 m.

### 2.5.1. - Integrated velocity model (2.1.)

We have determined for each experiment the drag coefficient $C_{10}$, the Coriolis parameter $f$, the damping parameter $\nu_\alpha$ for $\alpha = 0$ and $\alpha = \frac{1}{2}$ and also the best value of $\alpha$. We have calculated a non-dimensional parameter $\text{Cri}$ which tests the validity of the method.

$$ \text{Cri} = \frac{\int_{0}^{T} (u_{\emptyset} - u)^2 + (v_{\emptyset} - v)^2 \, dt}{\int_{0}^{T} (u^2 + v^2) \, dt} $$

(when $u_{\emptyset} = \mathcal{U}$ , $\text{Cri} = 0$)
We also give the ratio of the calculated kinetic energy (calcu. E) over the observed kinetic energy (obs. E). All these parameters are presented in Tab. 1.

We remark that the parameter $C_{ri}$ is smaller and the fit better when the initial value $\mathcal{U}(0)$ is calculated as a parameter ($C_{ri1}$ in Tab. 1), and not taken equal to the first data ($C_{ri2}$ in Tab. 1). Therefore, every parameter presented in Tab. 1 is calculated using the calculated initial value technique (Remark 2-6).

The only parameter we can get theoretically with accuracy is the Coriolis parameter ($f = .980 \ 10^{-4} \ s^{-1}$, $T = 17^{th} \ 50$ at the Bouée-Laboratoire). The comparison between the theoretical and the calculated Coriolis parameter gives us an idea of the accuracy of the method. The worst discrepancy between these two parameters is of about 5%.

The drag coefficient $C_{10}$ seems to be an increasing function of the wind speed or more precisely of the strength of the gust of wind. The highest value of $C_{10}$ was obtained during the gust of mistral of "COFRASOV". During "COBLAMED" the wind speed was smaller, and so was the drag coefficient.

Another point of interest is that the dissipations are important, specially during "COFRASOV" when the mixed layer is well formed. In this case turbulence is fully developed and an important part of the energy given by the wind is transformed into turbulent kinetic energy. There is no significant improvement of the accuracy as $\alpha$ increases. Therefore we can conclude that the linear damping is a good representation of the physics of dissipations.

We have drawn the E-W and N-S components of observed (dotted line) and calculated (full line) mean currents for "COFRASOV" in fig. 3, and for "COBLAMED" in fig. 4. We see the calculated current fits the observed current very well when the wind speed is high, and poorly when the wind speed is small (from 07 - 23 - 69 to 07 - 29 - 69 $12^h$ for "COFRASOV" experiment). We have observed a long trend on "COBLAMED" E - W component. This trend is not included in our model. In fig. (5) we show the same graphs for the gust of mistral of "COFRASOV" with $\alpha = 0$, (linear damping).
COFRASOV (07 - 18 - 69 - to 07 - 31 - 69)

<table><thead><tr><th></th><th>C<sub>10</sub></th><th>f</th><th>α</th><th>α̂</th><th>Calcul E./Obs.E.</th><th>Cri<sub>1</sub></th><th>Cri<sub>2</sub></th></tr></thead><tbody><tr><td>07.18.69</td><td>1.01 10<sup>-3</sup></td><td>1.002 10<sup>-4</sup>s<sup>-1</sup></td><td>0</td><td>0.9 10<sup>-6</sup>s<sup>-1</sup></td><td>0.68</td><td>0.33</td><td>0.36</td></tr><tr><td>to</td><td></td><td>(T=17h25)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>07.31.69</td><td>1.01 10<sup>-3</sup></td><td>1.005 10<sup>-4</sup>s<sup>-1</sup></td><td>0.5</td><td>0.25 10<sup>-6</sup></td><td>0.67</td><td>0.33</td><td>0.39</td></tr><tr><td>gust of</td><td>1.62 10<sup>-3</sup></td><td>0.9462 10<sup>-4</sup>s<sup>-1</sup></td><td>0</td><td>0.6 10<sup>-5</sup>s<sup>-1</sup></td><td>0.94</td><td></td><td>0.09</td></tr><tr><td>tristral</td><td></td><td>(T=18h27)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>07.26.69</td><td>1.62 10<sup>-3</sup></td><td>0.9467 10<sup>-4</sup>s<sup>-1</sup></td><td>0.5</td><td>0.7 10<sup>-6</sup></td><td>0.94</td><td>0.06</td><td>0.09</td></tr><tr><td>(12h)</td><td></td><td>(T=18h26)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>to</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>07.31.69</td><td>1.69 10<sup>-3</sup></td><td>0.9530 10<sup>-4</sup>s<sup>-1</sup></td><td>5</td><td>0.18 10<sup>-14</sup></td><td>0.93</td><td>0.065</td><td>0.08</td></tr><tr><td>(12h)</td><td></td><td>(T=18h19)</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

COBLAMED (09 - 17 - 69 - to 10 - 02 - 69)

<table><caption>Tab. 1</caption><thead><tr><th></th><th>C<sub>10</sub></th><th>f</th><th>α</th><th>α̂</th><th>Calcul E./Obs.E.</th><th>Cri<sub>1</sub></th><th>Cri<sub>2</sub></th></tr></thead><tbody><tr><td>0.72 10<sup>-3</sup></td><td>1.018 10<sup>-4</sup>s<sup>-1</sup></td><td>0</td><td>0.135 10<sup>-5</sup>s<sup>-1</sup></td><td>0.66</td><td>0.34</td><td>0.36</td></tr><tr><td>0.84 10<sup>-3</sup></td><td>1.02010<sup>-4</sup>s<sup>-1</sup></td><td>0.5</td><td>0.43 10<sup>-6</sup></td><td>0.66</td><td>0.34</td><td>0.35</td></tr></tbody></table>
COFRASOV

JULY 69

N. S. COMPONENT

CURRENT

<figure><img src="image_4.png" /></figure>

CURRENT

E.W. COMPONENT

<figure><img src="image_5.png" /></figure>

$v = 0.25 \ 10^{-6}$

$f = 1.005 \ 10^{-4}$

$C_{10} = 1.01 \ 10^{-3}$

$C_{r1} = 0.33$

Fig. 3 - Comparison of the observed mean current (dotted line) with the calculated mean current (full line) during the COFRASOV experiment.
COBLAMED

SEPTEMBER 69

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4 - Comparison of the observed mean current (dotted line) with the calculated mean current (full line) during the gust of Mistral of COFRASOV.</figcaption></figure>
COFRASOV

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5 - Comparison of the observed mean current (dotted line) with the calculated mean current (full line) during the gust of Mistral of COFRASOV.</figcaption></figure>
2.5.2. - Depth dependent model (2.15.)

We have also considered another model where the current
was not integrated in z but a function of the depth : hence **u** and **u**<sub>φ</sub> are
function of z and t.

$\mathcal{U}(z,t)$ is the solution of the following state equation.

$$
(2.15.) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t}(z,t) - fv(z,t) = \nu \frac{\partial^2}{\partial z^2} u(z,t) \\[1em]
\displaystyle \frac{\partial v}{\partial t}(z,t) + fu(z,t) = \nu \frac{\partial^2}{\partial z^2} v(z,t)
\end{array}
\right.
$$

with the boundary conditions : $\nu \frac{\partial \mathcal{U}}{\partial z}(0,t) = C_{10} \frac{\rho_a}{\rho_w} |W| W(t)$

$$
\frac{\partial \mathcal{U}}{\partial z}(-H, t) = 0.
$$

Where $v$ is kinematic turbulent viscosity

$\rho_a$ the air density

$\rho_w$ the water density

the initial condition is the following : At each gradient step
we build a solution $\mathcal{U}_0(z)$ of the form.

i) $\mathcal{U}_0(z_\phi) = \mathcal{U}_\phi(z_\phi)$

where $\mathcal{U}_{\phi}(z_{\phi})$ is the observed function at the observation point $z_{\phi}$.

$$
ii) \frac{\partial \mathcal{U}_0(0)}{\partial z} = C_{10} \frac{\rho_a}{\rho_w} |W| W(0)
$$

$$
\text{iii)} \quad \frac{\partial \mathcal{U}_0}{\partial z}(-H) = 0
$$

we have $\mathcal{U}(z,0) = \mathcal{U}_0(z) \forall z \in ]0, -h[$. We have used the same method as for the integrated problem. The time discretization is the same one, We have approximated the continuous problem in z by a finite difference method.
The results are presented in Tab. 2.

The drag coefficients $C_{10}$ are smaller than those found in the preceding case (2.5.1.). They are still an increasing function of the strength of the gust of wind.

The kinematic eddy viscosity coefficient $ν$ seems to increase as the wind speed increases. This agrees with the behaviour of $ν_α$ (2.5.1.).

The calculated Coriolis periods are the same as those calculated previously.

We have drawn the E-W and N-S component of observed (dotted line) and calculated (full line) currents for the gust of Mistral of COFRASOV in fig. 6.

This method is very interesting; but it requires a great number of current-meters to determine with a sufficient accuracy $\frac{\partial u}{\partial z}$, which it is not the case in these experiments.

## 2.6. - Conclusions

The method we have been dealing with in this paper is particularly efficient when the conditions of validity defined in the first part are filled.

There are only a few constraints on the structure of the partial differential equations (2.1.) or (2.15.). Thus the model could be extended without any difficulty. For example it might be possible to check the best law for the stress with respect to the wind.

Furthermore it must be pointed out that this method deals with the totality of the energy of the data.

Knowing the oceanographic parameters we can build an ocean - current forecast from the wind.

## 2.7 - Acknowledgements

This work was supported by the Centre National pour l'Exploitation des Oceans (C.N.E.X.O.) and by IRIA - LABORIA. The COFRASOV and COBLAMED experiment data were obtained by Dr.J.GONELLA and his team.
COFRASOV

<figure><img src="image_10.png" /><figcaption>Fig. 6 - Comparison of the observed currents (dotted line) and the calculated mean currents (full line) for the depth dependent model during the gust of Mistral of COFRASOV.</figcaption></figure>
REFERENCES

BRETSCHNEIDER (C.L.) (1966). On wind tides and longshore currents over the continental shelf due to wind blowing at an angle to the coast. National Engineering Science Co. SN 134-13.

CAHN A. (1945) An investigation of the free oscillations of a simple current system. J. Meteorology 2-113-119.

CSANADY G. (1972) Response of large stratified lake to wind. J. Physc. Oceanog 1. - p 3-13.

CREPON M. and D. CHOUCHAN (1974). Complex cross correlation coefficients and observations of oceanic current made on the bouée Lab. (Internal Report).

DURAND E. (1960) Solutions numériques des équations algébriques, Paris, Masson.

GLOWINSKI R. - LIONS J.L. - TREMOLIERES R. (1975) Application des méthodes d'optimisation de différences et d'éléments finis appliquées aux inéquations variationnelles. Paris, Dunod.

GONELLA J., M. CREPON and F. MADELAIN (1969). Observations de courants, vents et températures à la bouée Laboratoire. Cah. Oceanog. 21 p. 845-850.

GONELLA J. (1971). The drift current from observations made on the bouée Laboratoire. Cah. Oceanog. 23 p 1-15.

LIONS J.L. (1968). Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles. Paris, Dunod.

O'BRIEN J.J. and H.E. HURLBURT (1972). A numerical Model of Costal Upewelling J.Physi. Oceanog. I p.14 - 26.

POLLARD R.T. and R.C. MILLARD (1970) Comparison between observed and simulated wind - generated inertial oscillations - Deep-sea Res. 17 813-821.
POLLARD R.T., P. RHINES and R.O.R.Y. THOMPSON (1973). The deepening of the wind mixed layer. -Geophys. Fluid dynamics. V3. p. 381-404.

RATTRAY M. (1964). Time dependent motion in an ocean ; A unified two-layer, beta-plane approximation. Studies in Oceanog. p. 19-29.

R.D. RICHTMYER and K.W. NORTON. (1967). Difference methods for initial-value problems, Interscience publishers.

SUTTON W.G. (1934). The asymptotic expansion of a function whose operational equivalent is known J. LONDON Math. Soc. 9 131-137.

VARGA. R.S. (1962). Matrix Iterative Analysis. Prentice Hall, Englewood cliffs. N.J.

VERONIS G. and H. STOMMEL (1956). The action of variable wind stress on a stratified ocean. Mar. Res. 15 p. 43-75.