# Estimating the pycnocline depth from the SAR signature of internal waves in the Alboran Sea

Dessert Morgane, Marc Honnorat, Jean-Marc Le Caillec, Christophe

Messager, Xavier Carton

► To cite this version:

Dessert Morgane, Marc Honnorat, Jean-Marc Le Caillec, Christophe Messager, Xavier Carton. Estimating the pycnocline depth from the SAR signature of internal waves in the Alboran Sea. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2022, 15, pp.9048 - 9061. 10.1109/JSTARS.2022.3214298. [hal-03832959](https://hal.science/hal-03832959)

HAL Id: hal-03832959

https://hal.science/hal-03832959

Submitted on 4 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Estimating the Pycnocline Depth From the SAR Signature of Internal Waves in the Alboran Sea

Morgane Dessert ID, Marc Honnorat ID, Jean-Marc Le Caillec ID, Senior Member, IEEE, Christophe Messager ID, and Xavier Carton

**Abstract**—In the Alboran Sea, west of the Straits of Gibraltar, the pycnocline depth has been assessed from the signature of a large amplitude internal wave (LAIW) captured by a synthetic aperture radar (SAR) image. First, the coefficients of the extended Korteweg-deVriès model were expressed using two different models of ocean stratification: an interfacial model and a continuously stratified ocean model. Then, via a backscattering model, the same extended Korteweg-deVriès model was computed. This latter calculation was performed along several transects extracted from the LAIW surface signature and through an improved CARMA-derived method. Using the values of the coefficients for the ocean stratification model and those calculated from the SAR image, we obtained solutions for the pycnocline depth and thickness. This method was applied to an LAIW event on 1 October 2008, for which SAR data and stratification measurements by an in-situ experiment were available jointly. The results are that the interfacial stratification model provides only few solutions for the pycnocline depth, while the continuous stratification model allows an interval of solutions for the pycnocline depth and thickness. These models are nevertheless complementary. Extra applications of this method on other ocean regions would be of interest.

**Index Terms**—Alboran Sea, ocean internal waves, synthetic aperture radar (SAR).

## I. INTRODUCTION

OCEANS are often modeled as a two-layer ocean with lighter surface waters above the pycnocline and heavier waters below. This idealization is called the interfacial model of the ocean. In each layer, the water density is relatively homogeneous and the stratification (the vertical density gradient) is weak. The pycnocline acts as a boundary that mitigates the

vertical motions. It reduces exchanges between upper turbulent surface waters influenced by the atmospheric fluxes and the quieter, more nutrient-rich, and deeper waters. Since the pycnocline is the lower boundary of surface waters, its depth controls the upper ocean heat content for the ocean-atmosphere coupling processes. It also controls the vertical flux of nutrients necessary for phytoplankton blooms. These blooms play a key role in ecosystem regulation (first link in the food chain) and in climate variability (via carbon export). Besides, the strong stratification in the pycnocline interferes with acoustic waves through reflection/refraction and alters sonar measurements. All these points highlight the importance of a good assessment of the pycnocline depth.

Large amplitude internal waves (hereafter referred to as LAIWs) are key events in the modification of the pycnocline depth. Indeed, they interact with the pycnocline depth and they can be observed through remote sensing. Their surface signatures can be captured by satellite borne synthetic aperture radars (SARs) [1] but can be also noticed on optical sensors [2] using sunglint or in ocean color images [3], [4], [5]. Assuming an interfacial model of ocean stratification, i.e., a two-layer model, several authors [6], [7], [8], [9], [10], [11] interpreted the LAIW surface signature to extract ocean dynamics information. However, the actual ocean stratification sometimes differs from this two-layer model. Recently, a continuous stratification model (CSM) was used to predict the LAIW velocity and location in the Gulf of Maine [12]. In addition, most previous studies [9] focused on LAIWs with a soliton shape, whereas various shapes of LAIW exist. Moreover, some studies approximated the LAIW velocity with a linear approximation [11] while others imposed the pycnocline depth to assess the LAIW amplitude [9]. This article aims at circumventing the drawbacks and approximations of the aforementioned papers to estimate the pycnocline depth.

Considered as a “hot spot” of LAIWs, the Alboran Sea is the westernmost and one of the most biologically productive basins of the Mediterranean Sea [13]. It is connected to the Atlantic Ocean by the Straits of Gibraltar. The Alboran Sea is generally regarded as a two-layer system [14] with North Atlantic central waters at the surface and deep waters of the western Mediterranean Sea, separated by a thin pycnocline [15]. When specific stratification, currents, and tidal conditions are met [16], [17], LAIWs are generated in the Gibraltar Straits and propagated into the Alboran Sea. Besides, this region is regularly flown over by SAR satellites (RADAR-SAT, Sentinel) and a large database of SAR images has been available for years.

Manuscript received 20 May 2022; revised 26 September 2022; accepted 6 October 2022. Date of publication 13 October 2022; date of current version 26 October 2022. This work was supported by the Agence Innovation Defense (AID) from the Direction Générale de l'Armed (DGA); French Ministry of Defense Research Agency through a Ph.D grant. This paper was produced by the IEEE Publication Technology Group. They are in Piscataway, NJ, USA. (Corresponding author: Morgane Dessert.)

Morgane Dessert is with the Laboratory for Ocean Physics and Satellite Remote Sensing (LOPS), UBO, 29200 Plouzané, France, and also with the Extreme Weather Expertises (EXWEXs), 29200 Brest, France (e-mail: morgane.dessert@laposte.net).

Marc Honnorat and Christophe Messager are with the Extreme Weather Expertises (EXWEXs), 29200 Brest, France (e-mail: marc.honnorat@exwexs.fr; christophe.messager@exwexs.fr).

Jean-Marc Le Caillec is with the Image & Information Processing Department, IMT Atlantique Bretagne-Pays de la Loire, 29280 Plouzané, France (e-mail: jm.lecaillec@imt-atlantique.fr).

Xavier Carton is with the Laboratory for Ocean Physics and Satellite Remote Sensing (LOPS), UBO, 29200 Plouzané, France (e-mail: xcarton@univ-brest.fr).

Digital Object Identifier 10.1109/JSTARS.2022.3214298
<figure><img src="image_2.png" /><figcaption>Fig. 1. Normalized backscattered cross section (NRCS) σ₀ of the radar in dB from an SAR image of the eastern mouth of Gibraltar Straits. This image was captured on 1 October 2008 at 10:32 A.M. by ENVISAT (with a vertical polarization), with 13 m × 13 m nadir resolution. The localization of each CTD station is annotated with the time they were collected on 30 September 2008 (and the last CTD stations were acquired on 1 October 2008 at 00:20 A.M.). The five transects "T1" (northernmost) to "T5" (southernmost) are drawn as dashed black lines. Only the wavefronts of LAIW captured at 30 September 2008 at 22:17 P.M. are reported on the figure in thin dashed white lines (the entire image is not shown).</figcaption></figure>

All these elements make the Alboran Sea a privileged area for LAIW study.

From 30 September 2008 to 1 October 2008, the GIBRAL-TAR08 experiment at sea, aboard B.O. Sarmiento de Gamboa, performed in-situ measurements of conductivity-temperature-density (CTD) (see Fig. 1 where CTD locations are shown in white crosses and Fig. 2 where the corresponding measurements are depicted). The ENVISAT satellite acquired two SAR images of propagating LAIWs (see Fig. 1) just before and after the measurements. The first occurred on 30 September 2008 at 10:17 P.M. and the second occurred on 1 October 2008 at 10:32 A.M. The first image has a 13 × 13 m nadir-resolution and the second image has a 75 × 75 m nadir-resolution. The two images have vertical polarization (VV). These remote and in-situ measurements within a short time-lag allow us to study the relations between the SAR signal and the LAIW structure and dynamics.

As previously stated, this article aims at assessing the well-known two-layer interfacial model widely used for estimating the pycnocline depth from LAIW SAR signature. For this purpose, the coefficients of the equation governing the LAIW deformation/propagation are derived from the SAR images, through the following two different water column models:

1) the well-known interfacial two-layer ocean model with an infinitesimal pycnocline;

2) a continuous two-layer ocean model with a relatively thick pycnocline, for comparison.

In Section II, the two stratification models are presented. They are then applied to the following two situations:

1) a simulated surface signature of a theoretical LAIW in order to validate our approach;

2) the LAIW surface signature observed on the SAR image on 1 October 2008 in Section III.
<figure><img src="image_3.png" /><figcaption>Fig. 2. σ-Density, (the density ρ(z) minus 1000, blue line) and stratification N² (thin black line) derived from CTD profiles acquired from 30 September 2008 to 1 October 2008, during the GIBRALTAR08 experiment at sea, aboard B.O. Sarmiento de Gamboa. The hours of acquisition are indicated above each profile and permit to locate the station position on the map on Fig. 1. The last station (lower right) was acquired two times: at 23:24 on 30 September (solid line) and at 00:20 on 1 October (thin dashed line).</figcaption></figure>

Finally, Sections IV and V offer, respectively, a discussion and a conclusion. The main variables and their definitions are summarized for convenience in Table I.

## II. MODELS TYING THE PYCNOCLINE DEPTH AND THE LAIW SAR SIGNATURE

### A. Relation Between NRCS and the LAIW-Induced Surface Currents

SAR images provide the normalised radar cross section (NRCS) denoted σ₀ which is the ratio between the received and backscattered intensities (per surface unit). The NRCS value

results from several complex processes. The backscattered intensity depends on the frequency of the incident wave. Through resonance phenomena, the Bragg backscattering, the incident wave is backscattered preferentially by the ocean surface with wavelengths equal to half the incident SAR wavelength. These wavelengths correspond to the capillary waves range, from several centimeters to decimeters. Moreover, for Bragg waves generation, wind speed has to be higher than a threshold (approximately 2–3 ms⁻¹) and lower than roughly 10 ms⁻¹ [18]. Besides, the SAR incidence angle θᵢ must be larger than 10°. If these conditions are not met, then the wind and the LAIW patterns can no longer be separated. Under satisfying conditions,
<table><caption>TABLE I
MAIN VARIABLES AND DEFINITIONS</caption><thead><tr><th>Variables</th><th>Definitions</th></tr></thead><tbody><tr><td colspan="2"><strong>Backscattering model</strong></td></tr><tr><td>θ<sub>i</sub></td><td>SAR incidence angle</td></tr><tr><td>σ<sub>0</sub></td><td>Normalized Radar Cross Section (NRCS)</td></tr><tr><td>k<sub>0</sub></td><td>SAR wavenumber</td></tr><tr><td>k<sub>B</sub>, ω<sub>B</sub>, T<sub>B</sub></td><td>Bragg wavenumber, pulsation and period</td></tr><tr><td>c<sub>g</sub></td><td>Bragg waves group velocity</td></tr><tr><td>Ψ(⃗k)</td><td>surface sea spectrum</td></tr><tr><td>U<sub>0</sub>, U<sub>0,LAIW</sub></td><td>surface current, LAIW induced surface current</td></tr><tr><td>x<sub>r</sub>, x<sub>a</sub></td><td>range and azimuthal direction axis coordinate</td></tr><tr><td>x<sub>P</sub></td><td>LAIW propagation direction axis coordinate</td></tr><tr><td>x<sub>T</sub></td><td>axis coordinate in the direction perpendicular to the LAIW propagation direction</td></tr><tr><td colspan="2"><strong>Extended KdV model</strong></td></tr><tr><td>ξ</td><td>LAIW isopycnals elevation field</td></tr><tr><td>c</td><td>phase velocity under linear approximation</td></tr><tr><td>α<sub>1</sub></td><td>first order nonlinearity coefficients</td></tr><tr><td>α<sub>2</sub></td><td>second order nonlinearity coefficients</td></tr><tr><td>β<sub>1</sub></td><td>first order dispersivity coefficient</td></tr><tr><td>η</td><td>reference isopycnal elevation</td></tr><tr><td>N</td><td>Brunt-Väisälä pulsation</td></tr><tr><td>Φ(z)</td><td>linear weighted vertical function of isopycnal elevation</td></tr><tr><td>F(z)</td><td>first order nonlinear correction of the Φ(z) function</td></tr><tr><td colspan="2"><strong>Ocean stratification model</strong></td></tr><tr><td>ρ<sub>up</sub>, ρ<sub>bot</sub></td><td>mean density of the upper, bottom ocean layer</td></tr><tr><td>h<sub>1</sub>, h<sub>2</sub></td><td>upper, bottom layer thickness</td></tr><tr><td>H</td><td>total ocean depth</td></tr><tr><td>z<sub>p</sub>, t<sub>p</sub></td><td>pycnocline depth, thickness</td></tr><tr><td>η<sub>0</sub></td><td>maximum amplitude of pycnocline elevation</td></tr><tr><td colspan="2"><strong>Improved CARMA-derived model</strong></td></tr><tr><td>s, δ<sub>s</sub></td><td>discrete moving coordinate, sampling step of s</td></tr><tr><td>IP(s)</td><td>Integrated Profile along the transect extracted from the LAIW surface signature</td></tr><tr><td>a<sub>n</sub></td><td>weighting coefficients for the IP(s) decomposition</td></tr><tr><td>N</td><td>order of the IP(s) linear decomposition</td></tr><tr><td>M</td><td>order of the IP(s) quadratic and cubic decomposition</td></tr><tr><td>K</td><td>number of vertical modes</td></tr><tr><td>R</td><td>order up to which is described the vertical mode</td></tr></tbody></table>

$\sigma_0$ can be expressed as in [19]

$$ \sigma_0(\theta_i)_p = 16\pi k_0^4 |g_p(\theta_i)|^2 \Psi(\vec{k}_B) \quad (1) $$

where $\vec{k}_B$ is the Bragg wavevector (in what follows, the vectors are denoted as $\vec{k}$ while their corresponding norm as $k$). The Bragg wavenumber $k_B$, pulsation $\omega_B$, and period $T_B$ can be computed as

$$ k_B = 2 \cdot k_0 \cdot \sin(\theta_i) \quad (2a) $$

$$ \omega_B = \sqrt{g \cdot k_B} \quad (2b) $$

$$ T_B = \frac{2\pi}{\omega_B} \quad (2c) $$

where $k_0$ is the SAR wavenumber.

The $g_p$ functions are the first scattering coefficients for horizontal $p = HH$ or for the vertical $p = VV$ polarization

$$ g_{HH} = \frac{(\varepsilon_r - 1) \cos^2(\theta_i)}{[\cos(\theta_i) + \sqrt{\varepsilon_r - \sin^2(\theta_i)}]^2} \quad (3a) $$

$$ g_{VV} = \frac{(\varepsilon_r - 1) [\varepsilon_r(1 + \sin^2(\theta_i)) - \sin^2(\theta_i)] \cos^2(\theta_i)}{[\varepsilon_r \cos(\theta_i) + \sqrt{\varepsilon_r - \sin^2(\theta_i)}]^2} \quad (3b) $$

and $\varepsilon_r$ is the relative dielectric constant of seawater. $\Psi(\vec{k})$ is the surface sea spectrum. This spectrum can be computed either statistically or via models [20], [21], [22]. These models have been used to analyze the effects of the sea state on the LAIW surface signatures [23], [24], [25], [26] but are difficult to handle.

Since the space and time scales of the LAIWs are larger and slower than the space and time scales of the Bragg waves, the sea spectrum is given by the action balance equation:

$$ \begin{aligned} & \frac{\partial \Psi(\vec{k}_B)}{\partial t} + (c_g + U_0) \nabla \Psi(\vec{k}_B) \\ & = S_w(\vec{k}_B) + S_{w-w}(\vec{k}_B) + S_{c-w}(\vec{k}_B) + S_d(\vec{k}_B) \end{aligned} \quad (4) $$

where $\vec{k}_B$ is the Bragg wavevector, $c_g$ is the group velocity of Bragg waves, and $U_0$ is the surface current velocity. $U_0$ is considered as only induced by the LAIW in our case and will thus be denoted as $U_{0,LAIW}$. The right-hand side of (4) gathers the sources due to the wind $S_w$, the wave-wave interactions $S_{w-w}$, the current-wave interactions $S_{u-c}$, and the dissipation $S_{dis}$. The source functions are then simplified using the relaxation time approximation as in [27]. These authors assumed that the surface current $U_{0,LAIW}$ induces only a small modulation of the sea surface spectrum (in the weak hydrodynamic interaction framework). They simplified (1) and (4), and, thus, they induced a relation between the LAIW-induced backscattering anomaly $\bar{\sigma}_{0,LAIW}$ and surface current modulation as

$$ \bar{\sigma}_{0,LAIW}(x_r, x_a) = -4.5 \tau_r \frac{\partial U_{0,LAIW}(x_r, x_a)}{\partial x_r} - \frac{R}{V} \sin(\theta_i) \frac{\partial U_{0,LAIW}(x_r, x_a)}{\partial x_a} \quad (5) $$

where $\bar{\sigma}_{0,LAIW}(x_r, x_a) = (\sigma_0 - \bar{\sigma}_0)/\bar{\sigma}_0$ is the LAIW-induced NRCS anomaly where $\bar{\sigma}_0$ stands for the background NRCS (in the LAIW area within the SAR image), and $x_r$ is the range (across the radar track) axis coordinate, whereas $x_a$ is the azimuth (along the radar track) axis coordinate. $R$ and $V$ are respectively the ground-carrier distance and the platform velocity ($R/V \approx 130$ s for SEASAT satellite [27]). Finally, $\tau_r$ is the relaxation time conveying the source functions of (4). Physically, $\tau_r$ is the response time of the wave system to the current variation and was set to $\tau_r \approx 30 - 40$ s [27]. Then, as in [27], (5) is rewritten in the $(\vec{x}_P, \vec{x}_T)$ referential where $\vec{x}_P$ is the direction along the LAIW propagation direction, whereas $\vec{x}_T$ is the direction across the LAIW propagation direction. Thus, the rotation and the derivative introduce $\cos^2(\phi)$ and $\cos(\phi) \sin(\phi)$ factors with $\phi$ the angle between the satellite range direction $\vec{x}_r$ and the propagation direction $\vec{x}_P$ of the LAIW. Following the same approach, $\frac{\partial U_0}{\partial x_P}$ is assumed larger than $\frac{\partial U_0}{\partial x_T}$, implying that the $\bar{\sigma}_{0,LAIW}$ variations are negligible along $\vec{x}_T$. Equation (5) becomes

$$ \bar{\sigma}_{0,LAIW}(x_P) = -Q \frac{\partial U_{0,LAIW}(x_P)}{\partial x_P} \quad (6a) $$

$$ \text{with } Q = \left[ 4.5 \tau_r \cos^2(\phi) + \frac{R}{V} \cos(\phi) \sin(\phi) \sin(\theta_i) \right]. \quad (6b) $$
From (6), the surface current variation $\frac{\partial U_{0,\text{LAIW}}}{\partial x_P}$ can be deduced from the NRCS $\sigma_0$. The next sections describe the relation between the LAIW-induced current field $U_{0,\text{LAIW}}(x_P)$ and the density profile $\rho(z)$ in order to connect the NRCS and the density variation due to LAIWs.

## B. Relation Between the LAIW-Induced Surface Current and the Pycnocline Depth

1) **Definitions:** The following three main models are used for modeling LAIWs depending on the ratio between the horizontal scale $\lambda$ and the vertical scale $H$, which is the sea bottom depth:

1) the Korteweg–deVriès model when $\lambda > H$ [28];

2) the Benjamin–Ono model when $\lambda \approx H$ [29], [30];

3) the Joseph–Kubota model when $\lambda < H$ [31], [32].

The maximal depth of the ocean at the mouth of the Gibraltar Straits is around $H = 800$ m (GEBCO 2019). On 30 September 2008, SAR image (Fig. 1), the distance between bright and dark bands is around $\lambda = 2$ km, which is in accordance with typical horizontal length of the Alboran Sea [14]. Under these conditions ($\lambda > H$), the KdV theory is the most appropriate to study such LAIWs. This theory describes the evolution of the LAIW-induced isopycnal deformation. It is based on a balance between the wave dispersion with the wave steepening by nonlinear effects. These effects are conveyed by two parameters

$$ \alpha = \frac{A}{h} \quad \text{and} \quad \beta = \left(\frac{h}{\bar{\lambda}}\right)^2 \qquad (7) $$

where $\alpha$ conveys the nonlinear effects, $\beta$ is the dissipative effects, $A$ is the maximal pycnocline elevation induced by the LAIW, and $h$ is the pycnocline depth at rest. CTD sounding of 30 September 2008 provided in-situ density profiles in the region, and the ocean waters were stratified with a pycnocline depth ranging between 25 and 125 m depth (Fig. 2). These profiles were acquired at the same moment as an LAIW was propagating nearby as observed on the SAR image of Fig. 1. The pycnocline depth at rest $h$ is assumed to have usual values [14], i.e., $h \approx 20-50$ m and the magnitude of the LAIW-induced pycnocline elevation was $A \approx 1-10$ m. When $\alpha^2 \approx \beta$ (our case), among all the KdV-type equations, the Gardner or extended KdV (eKdV) equation is more appropriate than the classical KdV equation, for which $\alpha \approx \beta$. In the eKdV framework, the space–time variation of the LAIW isopycnal elevation field $\xi(x, z, t)$ is given by

$$ \xi_t + (c + \alpha_1 \cdot \xi + \alpha_2 \cdot \xi^2) \cdot \xi_x + \beta_1 \cdot \xi_{3x} = 0 \qquad (8) $$

where $c$ is the phase velocity under linear approximation, $\alpha_1$ and $\alpha_2$ are nonlinearity coefficients at first and second orders, $\beta_1$ is the dispersivity coefficient at the first order, and “t” or “x” indices stand for time or spatial derivatives (“3x” stands for third spatial derivative). The isopycnal depth anomaly or isopycnal excursion $\xi$ depends 1) on time $t$ as the soliton packet propagates, 2) on space $x$, and finally 3) on depth $z$ as the amplitude of the excursion depends on the stratification, which is stronger in the pycnocline. In particular, the LAIW-induced isopycnal elevation

$\xi$ can be expressed as a series as in [33]

$$ \frac{\xi(x, z, t)}{h} = \frac{\eta(x, t)}{h} \cdot \Phi(z) + \frac{\eta^2(x, t)}{h^2} \cdot F(z) + o\left(\frac{\eta^2}{h^2}\right) \qquad (9) $$

where $\eta(x, t)$, the isopycnal elevation, still satisfies the eKdV (8) and

$$ \Phi_{zz}(z) + \frac{N^2(z)}{c^2} \Phi(z) = 0 \qquad (10a) $$

$$ F_{zz}(z) + \frac{N^2(z)}{c^2} F(z) = -\frac{\alpha_1}{c} \cdot \Phi_{zz}(z) - \frac{3}{2} \frac{\partial}{\partial z} \left( \frac{\partial \Phi(z)}{\partial z} \right)^2 \qquad (10b) $$

$$ \text{with } N(z) = \sqrt{-\frac{g}{\rho(z)} \frac{d\rho(z)}{dz}} \qquad (10c) $$

where $N(z)$ is the Brunt-Väisälä pulsation and $g$ is the gravity acceleration. The function $\Phi(z)$ is the solution of Taylor-Goldstein equation [see (10a)], with $F(z)$ being the first-order nonlinear correction of the $\Phi(z)$ function. From the definition of the isopycnal deformation [see (9)], the LAIW-induced surface current can be expressed as [34]

$$ \frac{U_{0,\text{LAIW}}(x, t)}{h} = -c \cdot \Phi'(0) \cdot \frac{\eta(x, t)}{h} + \left[ \frac{\alpha_1}{2} \Phi'(z) + c \cdot F'(z) \right] \frac{\eta(x, t)^2}{h^2} + o\left(\frac{\eta^2}{h^2}\right). \qquad (11) $$

Assuming that $(\frac{\eta}{h})^2 \ll 1$, the isopycnal excursion [see (9)] and velocity [see (11)] will, from here, be expanded only to the first order (the linear one). Thus, combining the (6) and (11) leads to

$$ \text{IP}^*(x_P, t) = c \cdot \Phi'(0) \eta(x_P, t) = \frac{1}{Q} \int_{x_P - \epsilon}^{x_P + \epsilon} \bar{\sigma}_{0,\text{LAIW}}(x_0, t) dx_0 \qquad (12) $$

where IP* stands for (nonnormalized) integrated profile, $\epsilon$ defines the small range of the integration, and $Q$ is defined in (6).

2) **Horizontal Description:** In order to be simulated/analyzed, the eKdV (8) has to be normalized through a horizontal length $\lambda$ and a vertical length $H$ such as $\frac{\delta_{x_P}}{\lambda} \ll 1$ (where $\delta_{x_P}$ is the spatial sampling step along the LAIW propagation direction as discussed again in Section III-B1) and $\frac{\eta}{H} \ll 1$. Then, $\eta(x_P, t)$ function is expressed through a moving coordinate: $s = \frac{x_P - V \cdot t}{\lambda}$ since the LAIW is considered to be propagating as a mature wave, i.e., $\frac{\partial \eta}{\partial t} + V \cdot \frac{\partial \eta}{\partial x} = 0$ with $V$ being the LAIW velocity. Now, the eKdV can be expressed only through spatial derivatives

$$ \eta_s^* [c - V + \alpha_1 \cdot H \cdot \eta^* + \alpha_2 \cdot H^2 \cdot \eta^{*2}] + \frac{\beta_1}{\lambda^2} \eta_{3s}^* = 0 \qquad (13) $$

with $\eta^* = \frac{\eta}{H} \ll 1$.

Injecting (12) into (13) leads to a derived eKdV equation linking the NRCS $\sigma_0$ and the eKdV coefficients as

$$ \text{IP}_s(s)[C + A_1 \cdot \text{IP}(s) + A_2 \cdot \text{IP}^2(s)] + B_1 \text{IP}_3(s) = 0 \qquad (14) $$
with

$$ 
\begin{aligned}
C &= c - V \quad A_1 = \frac{H \cdot \alpha_1}{c \cdot \Phi'(0)} \\
A_2 &= \frac{H^2 \cdot \alpha_2}{c^2 \cdot \Phi'(0)^2} \quad B_1 = \frac{\beta_1}{\lambda^2}
\end{aligned}
$$

with $IP = \frac{IP^*}{H} \ll 1$ being the normalized integrated profile. Up to the first order, the derivatives can be approximated as

$$ D_1[\text{IP}](s) = \frac{\text{IP}(s) - \text{IP}(s - \delta_s)}{\delta_s} + o(\delta_s) \quad (16a) $$

$$ D_3[\text{IP}](s) = \frac{-\text{IP}(s - 2\delta_s)}{\delta_s^3} + 3\frac{\text{IP}(s - \delta_s)}{\delta_s^3} - 3\frac{\text{IP}(s)}{\delta_s^3} + \frac{\text{IP}(s + \delta)}{\delta_s^3} + o(\delta_s^3) \quad (16b) $$

where $\delta_s = \frac{\delta_{x_p}}{H} \ll 1$ is the normalized sampling spatial step (along the LAIW propagation direction) and, thus, (14) becomes a weighted summation of linear and nonlinear (quadratic and cubic) products of IP(s) at several lags of s

$$ 
\begin{aligned}
\text{IP}(s) = & \sum_{n=1}^{N} a_n \text{IP}(s - n\delta_s) \\
& + \text{IP}\left(s - \frac{N}{2}\delta_s\right) \cdot \sum_{m=1}^{M} a_{\frac{N}{2}, m} \text{IP}(s - m\delta_s) \\
& + \text{IP}^2\left(s - \frac{N}{2}\delta_s\right) \cdot \sum_{m=1}^{M} a_{\frac{N}{2}, \frac{N}{2}, m} \text{IP}(s - m\delta_s) + \epsilon_d(s)
\end{aligned}
\quad (17) $$

where $a_n$ are the weighting coefficients and $M$ and $N$ are the summation upper bounds and $\epsilon_d(s) = o(\delta_s) \ll 1$ is the independent error function having a null statistical mean (i.e., $\overline{\epsilon_d} = 0$). In what follows, the bar stands for statistical mean along the LAIW propagation direction. In order to reduce the error $\epsilon_d(s)$ due to integration [see (12)], a diffusive numerical scheme is used. This scheme implies that $M$ is even and $N$ is odd. Moreover, isolating the term $\text{IP}(s)$ [right-hand term in (17)] requires $N > M$.

Then, expression (17) is multiplied successively by

1) $\text{IP}(s - n \cdot \delta_s) > 0$ for $n = 1, \dots, N$;

2) $\text{IP}(s - (N - 1)\delta_s/2) \cdot \text{IP}(s - m \cdot \delta_s)$ for $m = 1, \dots, M$;

3) $\text{IP}^2(s - (N - 1)\delta_s/2) \cdot \text{IP}(s - m \cdot \delta_s)$ for $m = 1, \dots, M$.

The noninteger lags (as $s - (N - 1)\delta_s/2$) are obtained by averaging two adjacent samples in this matrix.

Finally, the average operator $\bar{\cdot}$ is applied to each equation leading to a linear system of $N + 2M$ equations. For writing convenience, the function $h(i, j, \dots, k)$ is introduced:

As a remark, the indices are permutative: $h_{i,j,\dots,k} = h_{k,\dots,i,j}$. Equation (19) shown at the bottom of this page, is the system derived for $N = 3$ and $M = 2$, and the $a$ coefficients are defined as

$$ a_1 = -a_2 = 3 - \frac{c - V}{\beta_1} \cdot \delta_s^2 \lambda^2 \quad (20a) $$

$$ a_3 = 1 \quad (20b) $$

$$ a_{\frac{3}{2},1} = -a_{\frac{3}{2},2} = -\frac{\alpha_1}{\beta_1} \cdot \frac{\delta_s^2 \lambda^2}{c \cdot \Phi'(0)} \quad (20c) $$

$$ a_{\frac{3}{2}, \frac{3}{2}, 1} = -a_{\frac{3}{2}, \frac{3}{2}, 2} = \frac{\alpha_2}{\beta_1} \cdot \frac{\delta_s^2 \lambda^2}{c^2 \cdot \Phi'(0)^2} \quad (20d) $$

On the one hand, coefficients $a$ can be calculated from $\text{IP}(s)$ values by solving (19) shown at the bottom of this page, through an LU decomposition. This method is close to the CARMA method already used for LAIW in 2006 by [11]. The differences are as follows: 1) the coherence between the discretized eKdV equation order and the constraining order $N$ or $M$ in selecting the neighboring IP values and 2) the use of Yule-Walker equations (detailed in [35]) as solving method.

On the other hand, coefficients $a$ [right-hand terms of equations from (20a) to (20d)] can also be expressed from $\Phi'(z)$, $c$, $\alpha_1$, and $\beta_1$, which depend on the vertical description of the problem.

Thus, by solving (19), we can estimate the LAIW propagation parameters (20) and then we can connect them to the water column parameters, in particular, the pycnocline depth. In the next section, we develop this approach to two stratification models.

3) Vertical Description: Oceans are often modeled through two idealized vertical stratification profiles: the two-layer interfacial ocean model and the continuously stratified ocean model.

a) Interfacial Stratification Model: The interfacial model is widely used in LAIW modeling. It assumes that the ocean is made up of two layers of constant density where $\rho_{up}$ is the upper

$$ 
\left( \begin{array}{c} h_{0,1} \\ h_{0,2} \\ h_{0,3} \\ h_{0,\frac{3}{2},1} \\ h_{0,\frac{3}{2},2} \\ h_{0,\frac{3}{2},\frac{3}{2},1} \\ h_{0,\frac{3}{2},\frac{3}{2},2} \end{array} \right) = \left( \begin{array}{cccc} h_{1,1} & h_{1,2} & h_{1,3} & h_{1,\frac{3}{2},1} & h_{1,\frac{3}{2},2} & h_{1,\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},2} \\ h_{2,1} & h_{2,2} & h_{2,3} & h_{2,\frac{3}{2},1} & h_{2,\frac{3}{2},2} & h_{2,\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},2} \\ h_{3,1} & h_{3,2} & h_{3,3} & h_{3,\frac{3}{2},1} & h_{3,\frac{3}{2},2} & h_{3,\frac{3}{2},\frac{3}{2},1} & h_{3,\frac{3}{2},\frac{3}{2},2} \\ h_{1,\frac{3}{2},1} & h_{1,\frac{3}{2},2} & h_{1,\frac{3}{2},3} & h_{1,\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},2} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},2} \\ h_{2,\frac{3}{2},1} & h_{2,\frac{3}{2},2} & h_{2,\frac{3}{2},3} & h_{2,\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},2} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},2} \\ h_{1,\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},2} & h_{1,\frac{3}{2},\frac{3}{2},3} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},2} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{1,\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},2} \\ h_{2,\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},2} & h_{2,\frac{3}{2},\frac{3}{2},3} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},2} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},1} & h_{2,\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2},2} \end{array} \right) \cdot \left( \begin{array}{c} a_1 \\ a_2 \\ a_3 \\ a_{1,\frac{3}{2}} \\ a_{2,\frac{3}{2}} \\ a_{1,\frac{3}{2},\frac{3}{2}} \\ a_{2,\frac{3}{2},\frac{3}{2}} \end{array} \right) \quad (19) $$
layer density, while $\rho_{\text{bot}}$ is the lower layer density

$$ \rho(z) = \begin{cases} \rho_{\text{up}} & \text{if } z < h_1 \\ \rho_{\text{bot}} & \text{otherwise} \end{cases}. \qquad (21) $$

The two layers are separated by an interfacial pycnocline (a null thickness pycnocline) at depth $z = h_1$. Under interfacial approximation, (11) and (20) can be simplified since $\Phi'(0) \to -\frac{1}{h_1}$, thus leading to coefficients

$$ c = \sqrt{\frac{g\sigma_{\rho}h_1h_2}{h_1 + h_2}} \qquad (22a) $$

$$ \alpha_1 = \frac{3}{2}c \frac{h_1 - h_2}{h_1 h_2} \qquad (22b) $$

$$ \alpha_2 = \frac{3c}{(h_1 h_2)^2} \left[ \frac{7}{8} (h_1 - h_2)^2 - \frac{h_1^3 + h_2^3}{h_1 + h_2} \right] \qquad (22c) $$

$$ \beta_1 = \frac{c}{6} h_1 h_2 \qquad (22d) $$

where $h_2 = H - h_1$ is the thickness of the bottom layer, $\sigma_{\rho} = \frac{\rho_{\text{up}} - \rho_{\text{bot}}}{\bar{\rho}}$ is the density jump, related to the stiffness of the pycnocline and $\bar{\rho} = \frac{\rho_{\text{up}} + \rho_{\text{bot}}}{2}$ is the average density over the water column (see [36]). Combining (22b), (22d), and (20c) leads to a relation between the $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$ coefficients and the pycnocline depth $h_1$

$$ a_{\frac{3}{2},1} = -a_{\frac{3}{2},2} = -\delta_s^2 \lambda^2 \cdot h_1 \cdot \left[ \frac{\alpha_1}{\beta_1 c} \right] \qquad (23a) $$

$$ a_{\frac{3}{2}, \frac{3}{2}, 1} = -a_{\frac{3}{2}, \frac{3}{2}, 2} = \delta_s^2 \lambda^2 \cdot h_1^2 \cdot \left[ \frac{\alpha_2}{\beta_1 c^2} \right] \qquad (23b) $$

with $\alpha_1$, $\alpha_2$, $\beta_1$, and $c$ nonlinearly depending on $h_1$. Combining (23) and (19) constitutes the interfacial stratification model (ISM) linking the NRCS $\sigma_0$ (extracted from the SAR image) and $h_1$.

b) Continuous Stratification Model: When the pycnocline thickness is too large to be considered as an interface, $\Phi(z)$ in (10a) must be computed by solving an eigenvalue Sturm-Liouville problem with boundary conditions $\Phi(0) = \Phi(H) = 0$. This system has an infinite number of solutions defined by the pairs of eigenvalues $\varphi^{(k)}$ and eigenfunctions (also called vertical modes) $W^{(k)}(z)$ as

$$ \Phi(z) = \sum_{k=1}^{K} \varphi^{(k)} \cdot W^{(k)}(z) \qquad (24) $$

where $K$ is the order up to which the function $\Phi(z)$ is described (i.e., the number of vertical modes). Each vertical mode $W^{(k)}(z)$ is described as

$$ W^{(k)}(z) = \sum_{r=1}^{R^{(k)}} \gamma_r^{(k)} \cdot \sin\left(\frac{r\pi z}{H}\right) \qquad (25) $$

where $R^{(k)}$ is the order up to which the $k$th mode is described. The $k$th vertical mode has a propagation velocity $c^{(k)}$ that can be computed from the $k$th eigenvalue $\varphi^{(k)}$ through

$$ c^{(k)} = 1/\sqrt{\varphi^{(k)}}. \qquad (26) $$

For this reason, we omit the (k) subscript. This first vertical mode $W(z)$ has only one maximum, equal to 1, the argument of which is denoted as $z_{\text{max}}$. In (25), R is not straightforwardly determined since it depends on the stratification. R is required to be higher for a shallower pycnocline than for a pycnocline closer to the half total depth. R is optimal when

$$ \frac{z_{\text{max}} - z_p}{z_p} \ll 1 \qquad (27) $$

and when the quantities

$$ \left[ \frac{\alpha_1}{\beta_1} \cdot \frac{\delta_s^2 \lambda^2}{c \cdot \Phi'(0)} \right] \quad \text{and} \quad \left[ \frac{\alpha_2}{\beta_1} \cdot \frac{\delta_s^2 \lambda^2}{c^2 \cdot \Phi'(0)^2} \right] \qquad (28) $$

remain almost unchanged when increasing R. In (27), $z_p$ is the center of the pycnocline and has a close geophysical meaning to $h_1$ for the ISM. In other words, the above conditions are met when the depth of maximum $\Phi$ and the depth of the middle of the pycnocline are relatively close [see (27)] and when $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$ parameters are unchanged with R [see (28)]. As a remark, the situations where the middle of the pycnocline depth is either close to $\frac{H}{2}$ or close to 0 are not considered in this model since not realistic in our situation. Since the sum of (25) contains many unknowns (all the $\gamma_r$), then a continuous density profile, clearly separating two layers, is introduced as in [37]

$$ \rho(z) = \rho_{\text{up}} \cdot \left( \frac{\rho_{\text{bot}}}{\rho_{\text{up}}} \right)^{\frac{\varrho(z)}{\varrho(H)}} \qquad (29a) $$

$$ N^2(z) = \frac{g \cdot \log \left( \frac{\rho_{\text{bot}}}{\rho_{\text{up}}} \right) \cdot t_p}{\left( t_p^2 + (z - z_p)^2 \right) \cdot \varrho(H, z_p, t_p)} \qquad (29b) $$

$$ \text{where } \varrho(z) = \arctan \left( \frac{z_p}{t_p} \right) + \arctan \left( \frac{z - z_p}{t_p} \right) \qquad (29c) $$

where $t_p$ expresses the pycnocline thickness.

This approach involves introducing $N(z)$ into (10a) and retrieving $\gamma_r$ and then calculating the pair $(\varphi, W(z))$ by a variational approach [38]. Apart from decreasing the number of unknowns, this density profile has the advantage to be easily fitted to most two-layer oceanic stratification conditions. Moreover, this expression can also be analytically integrated for a continuously stratified model, in which the eKdV coefficients are rewritten as in [33]

$$ \alpha_1 = \frac{3}{2I} c \int_0^H \Phi'(z)^3 dz \qquad (30a) $$

$$ \alpha_2 = \frac{1}{2I} \cdot \int_0^H -6c\Phi'(z)^4 + 5\alpha_1\Phi'(z)^3 - \frac{(\alpha_1)^2}{c}\Phi'(z)^2 dz \qquad (30b) $$

$$ \beta_1 = \frac{1}{2I} c \int_0^H \Phi(z)^2 dz \qquad \text{where } I = \int_0^H \Phi'(z)^2 dz \quad \text{and} \quad c = \frac{1}{\sqrt{\varphi}} \qquad (30c) $$

where $\alpha_1$, $\alpha_2$, $\Phi'(0)$, and $c$ can then be calculated from ($z_p, t_p$); but as seen in Section III-A3, the two pycnocline parameters
<figure><img src="image_4.png" /><figcaption>Fig. 3. Description of the methodology for both models: solid light blue arrows stand for the validation experiment through ISM while dashed orange arrows act for the validation experiment through CSM. The realistic application of the method on extracted transects is indicated through green thin dashed line. Further explanation can be read in the text.</figcaption></figure>

cannot be uniquely estimated from these coefficients. Combining (30a), (30b), (30c), and $\Phi(z)$ and $c$ calculated from (25) and (26) leads to

$$a_{\frac{3}{2},1} = -a_{\frac{3}{2},2} = -\delta_s^2 \lambda^2 \cdot \left[ \frac{\alpha_1}{\beta_1 c \Phi'(0)} \right] \quad (31a)$$

$$a_{\frac{3}{2},\frac{3}{2},1} = -a_{\frac{3}{2},\frac{3}{2},2} = \delta_s^2 \lambda^2 \cdot \left[ \frac{\alpha_2}{\beta_1 c^2 \Phi'(0)^2} \right]. \quad (31b)$$

The $a_{\frac{3}{2},}$ and $a_{\frac{3}{2},\frac{3}{2},}$ coefficients are expressed through $\alpha_1$, $\alpha_2$, $\beta$, $c$, and $\Phi'(0)$ which all depend on the stratification parameters ($z_p$, $t_p$). Combining (31) and (19) constitutes the CSM binding the NRCS $\sigma_0$ and ($z_p$, $t_p$).

### III. RESULTS

The ISM and CSM presented above have been applied to two situations. The methodology is illustrated in Fig. 3. As a first experiment, a theoretical surface signature is simulated from a theoretical stratification situation (the density profile is drawn

<figure><img src="image_5.png" /><figcaption>Fig. 4. Description of the soliton experiment with the conditions $z_p = 56$ m, $t_p = 4$ m, and $\eta_0 = -10$ m. Left panel: density profile (dashed light blue line) and vertical mode (solid orange line) developed to the 70th order ($R = 70$) where maximum lays at $z_{max} = 68.8$ m. Upper right panel: the excursion of one simulated isopycnal with added zero-averaged noise (dark blue line). Middle right panel: simulated NRCS and smoothed simulated NRCS. Lower right panel: IP function.</figcaption></figure>

in light blue on the left panel of Fig. 4). This totally defined experiment acts as a validation of the method presented in the previous sections. From the theoretical profile, the eKdV coefficients are first calculated. Then, from the eKdV coefficients, the (Gardner soliton type) excursion of the referential isopycnal is defined and zero-averaged noise is added to evaluate the stability of our method. Surface signature is then calculated from this simulated isopycnal excursion. This simulated surface signature is first smoothed [through a Gaussian filter, the width of which is iteratively adjusted to satisfy $a_3 \approx 1$, since this coefficient is theoretically equal to 1, (20b), but the solving of system (19) provides only a close value, see Fig. 3 for the algorithm details]. Then, the method is applied to the denoised surface signature to obtain the $a$ coefficients. Finally, the ISM and CSM are used in order to estimate the pycnocline depth. The method is also applied to real NRCS surface signature extracted from the SAR image (Fig. 1) considering both ISM and CSM as detailed in Section III-B. This part of the method is indicated through green thin dashed lines in Fig. 3.

#### A. LAIW Soliton Experiment

The ISM and CSM are first evaluated on a known situation satisfying all the hypotheses detailed above: a LAIW Gardner-type soliton propagating in a two-layer ocean.

1) Stratification Simulation: Fig. 4 illustrated at left panel the density profile $\rho(z)$ [defined as in (29)]. At upper right panel,
<table><caption>TABLE II
PARAMETERS FOR SIMULATED LAIW PROPAGATION</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>ρ<sub>up</sub></td><td>1026 kg/m<sup>3</sup></td></tr><tr><td>ρ<sub>bot</sub></td><td>1030.8 kg/m<sup>3</sup></td></tr><tr><td>h<sub>1</sub></td><td>56 m</td></tr><tr><td>t<sub>p</sub></td><td>4 m</td></tr><tr><td>η<sub>0</sub></td><td>10 m</td></tr><tr><td>H</td><td>800m</td></tr></tbody></table>

This values were chosen to be close to the realistic situation.

<table><caption>TABLE III
EKDV ESTIMATED PARAMETERS FOR THE SIMULATED LAIW PROPAGATION (CSM VALUES ARE COMPUTED FOR z<sub>p</sub> = 56 M, t<sub>p</sub> = 4 M, AND R = 70)</caption><thead><tr><th>Parameters</th><th>ISM values</th><th>CSM values</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>-0.02471</td><td>-0.02255</td></tr><tr><td>α<sub>2</sub></td><td>-0.0001738</td><td>-0.0001689</td></tr><tr><td>β<sub>1</sub></td><td>6926</td><td>7785</td></tr><tr><td>c</td><td>0.9975 m.s<sup>-1</sup></td><td>0.9491 m.s<sup>-1</sup></td></tr><tr><td>z<sub>max</sub></td><td></td><td>68.8 m</td></tr><tr><td>a<sub>3/2</sub></td><td>-0.08010</td><td>-0.06683</td></tr><tr><td>a<sub>3/2,3/2</sub></td><td>-0.03164</td><td>-0.02888</td></tr></tbody></table>

the isopycnal excursion $\eta(s)$ is defined as in [36]

$$ \eta(s) = \frac{\eta_0}{b + (1 - b) \cdot \cosh^2\left(\frac{s}{\Delta}\right)} \quad (32a) $$

$$ \text{where } \Delta = \sqrt{\frac{12\beta_1}{\eta_0\left(\alpha_1 + \frac{\eta_0 \cdot \alpha_2}{2}\right)}} \quad (32b) $$

$$ \text{and } b = \frac{-\eta_0 \cdot \alpha_2}{2\alpha_1 + \alpha_2 \cdot \eta_0} \quad (32c) $$

where $\eta_0$ is the maximum amplitude. The maximum amplitude is set to $\eta_0 = 10$ m (usual value). The experiment parameters are summed up in Table II. Before applying any method, the eKdV parameters of the experiment are first calculated for the ISM and the CSM (Table III).

In order to apply the CSM approach, the function $\Phi(z)$ is expanded up to the optimal order $R_{opt}$. In Fig. 5, $\frac{z_{max}-z_p}{z_p}$, $a_{3/2}$, and $a_{3/2,3/2}$, are shown for eight stratification situations: $z_p = -60$ m, -45 m, -30 m, and -15 m while $t_p = 4$ or 8 m, in order to evaluate the departures between the stratification situations and to estimate $R_{opt}$ [under conditions (27) and (28)]. Moreover, the thicker the pycnocline is, the higher the $R$ must be to satisfy $\frac{z_{max}-z_p}{z_p} \ll 1$. This condition is satisfied only for the deepest pycnocline in the sensitivity experiment (Fig. 5). Beyond $R = 70$, the $a_{3/2}$, and $a_{3/2,3/2}$, parameters vary little (especially for a deep $z_p < -30$ m and a thick pycnocline). This $R$ is set to $R_{opt} = 70$ to compute the coefficients in Table III. The eKdV coefficients exhibit some differences between ISM and CSM approaches:

1) $\frac{|\alpha_1(\text{ISM}) - \alpha_1(\text{CSM})|}{\alpha_1} < 1\%$;

2) $\frac{|\alpha_2(\text{ISM}) - \alpha_2(\text{CSM})|}{\alpha_2} \approx 3\%$;

3) $\frac{|c(\text{ISM}) - c(\text{CSM})|}{c} \approx 5\%$;

4) $\frac{|\beta_1(\text{ISM}) - \beta_1(\text{CSM})|}{\beta_1} \approx 12\%$;

5) $\frac{|h_1(\text{ISM}) - z_{\text{max}}(\text{CSM})|}{z_p} > 20\%$.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Sensitivity of $\frac{z_{max}-z_p}{z_p}$ (upper), $\frac{\alpha_1}{\Phi'(0)\beta_1 c}$ (middle) and $\frac{\alpha_2}{\Phi'(0)^2\beta_1 c^2}$ (lower) regarding the order of the first vertical mode definition $R$ for different stratification situation : ($z_p, t_p$) = (-15, 4) in solid blue line; ($z_p, t_p$) = (-15, 8) in dashed blue line; ($z_p, t_p$) = (-30, 4) in orange solid line; ($z_p, t_p$) = (-30, 8) in orange dashed solid line; ($z_p, t_p$) = (-45, 4) in solid green line; ($z_p, t_p$) = (-45, 8) in dashed green line; ($z_p, t_p$) = (-60, 4) in yellow solid line; ($z_p, t_p$) = (-60, 8) in yellow dashed solid line.</figcaption></figure>

<table><caption>TABLE IV
a COEFFICIENTS FOR THE VALIDATION EXPERIMENT</caption><thead><tr><th>Coefficient</th><th></th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>3.0033</td></tr><tr><td>a<sub>2</sub></td><td>-3.0033</td></tr><tr><td>a<sub>3</sub></td><td>1.0000</td></tr><tr><td>a<sub>3/2,1</sub></td><td>0.06106</td></tr><tr><td>a<sub>3/2,2</sub></td><td>-0.06086</td></tr><tr><td>a<sub>3/2,3/2,1</sub></td><td>0.04711</td></tr><tr><td>a<sub>3/2,3/2,2</sub></td><td>-0.04627</td></tr></tbody></table>

In particular, $\beta_1$ value discrepancies are significant. Thus, the eKdV equation parameters, and then the SAR signature of the LAIW, are sensitive to the stratification model.

2) *Simulated Surface Signature:* From the isopycnal excursion with added zero-averaged noise (Fig. 3 upper right panel), the density profile (Fig. 3 left panel), and (11), the surface signature of the soliton is simulated and then filtered with a Gaussian filter (Fig. 4 middle right panel). Finally, it is integrated to get IP (lower right panel of Fig. 4). The $a$ coefficients are then computed through the pseudo-CARMA method [see (19)]. Based on the condition (20b), the $a$ coefficients are iteratively computed while adjusting the width of the Gaussian filter. The resulting coefficients are gathered in Table IV.

3) *ISM/CSM Pycnocline Depth Estimations:* For the ISM, several solutions are obtained for the pycnocline depth; only positive and real solutions are collected in Table V. The pycnocline depth from (23a) lays at $h_1 = 64, 7$ m (that is 8.6 m or 15% too deep), whereas (23b) leads to too shallow a pycnocline:
<table><caption>TABLE V
$h_1$ ESTIMATE FROM $a$ FOR SIMULATED NRCS AND ISM</caption><thead><tr><th>$h_1$ estimation</th><th></th></tr></thead><tbody><tr><td>from $a_{\frac{3}{2}}$</td><td>64.7 m</td></tr><tr><td>from $a_{\frac{3}{2}} \cdot \frac{3}{2}$</td><td>24.2 m</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Coefficients $\frac{\alpha_1}{\beta_1 c \Phi'(0)}$ for upper panel and $\frac{\alpha_2}{\beta_1 c^2 \Phi'(0)^2}$ for the lower panel for a wide range of the pycnocline depths $z_p$ and thicknesses $t_p$. Black point indicates the initial value set for the soliton experiment. Black (white) dashed lines indicate the CSM solutions $\frac{\alpha_1}{\beta_1 c \Phi'(0)} = \frac{a_{\frac{3}{2}}}{\Delta x^2}$ (for upper panel) and $\frac{\alpha_2}{\beta_1 c^2 \Phi'(0)^2} = \frac{a_{\frac{3}{2}}}{\Delta x^2}$ (for lower panel) following the solutions of couples $(z_p, t_p)$ which corresponds to the surface SAR signature of the soliton experiment (to the surface SAR signature extracted through the five transects).</figcaption></figure>

$h_1 = 24.2$ m. For the CSM, Fig. 6 shows the values of parameters $\frac{\alpha_1}{\phi'(0)\beta_1 c}$ (upper panel) and $\frac{\alpha_2}{\phi'(0)^2\beta_1 c^2}$ (lower panel) computed for a range of $z_p$ and $t_p$ and for $R = 70$. Black lines depict the values of $(z_p, t_p)$ for which the coefficients $a_{\frac{3}{2}}$ and $a_{\frac{3}{2}} \cdot \frac{3}{2}$, are constant [see (31)], while the black point indicates the real initial values ($z_p = 56$ m and $t_p = 4$ m). For coefficients $a_{\frac{3}{2}}$, or $a_{\frac{3}{2}} \cdot \frac{3}{2}$, a single value leads to an infinity of solutions ($z_p, t_p$).

The value of the coefficient $a_{\frac{3}{2}}$ gives that the pycnocline depth lays between $-64 < z_p < -54$ m considering a range of thickness $t_p$ between 3 and 8 m; the value of the coefficient $a_{\frac{3}{2}} \cdot \frac{3}{2}$, implies a pycnocline depth between $-43 < z_p < -37$ m (always considering a thickness range $3 < t_p < 8$). Estimating the pycnocline depth from the only SAR surface signature seems intricate, but it still provides an insight on the subsurface stratification. Moreover, this highlights the error when considering an interfacial two-layer ocean.

For this experiment, the CSM leads to more coherent solutions for $z_p$. For both methods, the coefficients $a_{\frac{3}{2}}$ lead to more coherent estimates than the coefficients $a_{\frac{3}{2}} \cdot \frac{3}{2}$, do.

## B. October 2008 LAIW Event

In this section, we consider the two SAR images presented in the introduction as well as the CTD measurements. A simulation of the tidal currents [39] at this period suggested that the conditions for internal waves generation at Camarinal Sill were met several hours before and that those SAR images could depict a LAIW. Even if the LAIW SAR signature does not always appear as bright and dark successive bands, the arc-of-circle is considered as LAIW surface signature. However, considering the Alboran Sea as a two-layer ocean can seem unrealistic especially in light of the density profiles in Fig. 2. Indeed, for stations acquired at 15:54, 22:22, or 23:24, a strong stratification can be highlighted respectively at 100, 35, and 20 m. For the other stations, the profiles are composed of several stratification maxima. While exaggerating the two-layer model, they can be seen as profiles with a thick pycnocline encompassing all these stratification maxima.

1) *Transects Extraction:* Five transects (denoted from “T1” to “T5”) have been extracted from the 1 October 2008 SAR image (we chose this image due to a better resolution), perpendicularly to the arc-of-circle shaped wave fronts. In this area, the bathymetry is almost constant (Fig. 1). The resolution of the extracted transects $\delta_{x_p}$ is given by the SAR image resolution divided by $\cos(\phi)$ where $\phi$ is the angle between the SAR range axis and the LAIW propagation direction (see Section II-A). The resolution of the transects ranges are: 34.2 m (T1); 27.7 m (T2); 29 m (T3); 30.5 m (T4), and 26.7 m (T5). These transects are marked in Fig. 7 (in gray thin lines) with minima corresponding to the darkest bands and maxima corresponding to the brightest bands. Averaged wind speeds are estimated as 4–5 m.s$^{-1}$ along the transects, using a model [40], [41], [42] based on the CMOD5 model [43]. The incidence angle $\theta_i$ (satellite characteristics) of the SAR is around 24°. The averaged bathymetry (from GEBCO) $\overline{H}$ under the transects is $\overline{H} = 792 \pm 46$ m with low averaged slopes $\partial H/\partial x_p$ ranging from -0.64% to -1.32%. These values satisfy the wind, incidence angle, and weak bathymetry variation conditions mentioned above. The parameter $\tau_r$ is expected to be 10–100 Bragg wave periods. In 1984, Alpers and Hennings [27] calculated $\tau_r = 4.7 - 47$ s and $\overline{R}/\overline{V} = 130$ s for SEASAT (in Band L). For ENVISAT (Band C), the Bragg wavelength is $\lambda_B = 7.2$ cm and the Bragg period is $T_B = 0.21$ s, and $\tau_r$ is evaluated as 2.1 – 21 s. As the altitude of the satellite
<figure><img src="image_8.png" /><figcaption>Fig. 7. Transects extracted from 1 October 2008 SAR image. From top to down are gathered the transects from northern one (at top) to southern one (at bottom) as represented in Fig. 1. IP (black line referenced on right y-axis) and $\sigma_0/Q$ (gray thin line referenced on left y-axis) are illustrated for each transect.</figcaption></figure>

is not significantly different, the *R/V* value from [27] is also considered for ENVISAT.

For all the transects, one packet of solitons can be highlighted (toward the tail of the transect) following a higher isopycnal excursion (at the right part of the transect). The noise removed to the signal (by the Gaussian filtering) is considered zero-averaged (actually lower than $10^{-14}$).

2) Coefficients *a* and the Pycnocline Depth *z*<sub>p</sub> Estimations: Our method is then applied to the IP profiles computed through (12) (bold dark lines in Fig. 7). The resulting *a* coefficients are gathered in Table VI while Table VII summarizes the solutions for the pycnocline depth considering the ISM. In this case also, several solutions exist for the pycnocline depth, and only the positive and real solutions were collected. Only transect T2 does not lead to any solution for the ISM, suggesting the limitation of this model as previously noticed. All the $a_{3/2}$ coefficients calculated from transects were out of realistic range chosen for CSM except the transect T5. Only transects T5, T4, and T3

<caption>TABLE VI</caption>

*a* COEFFICIENTS COMPUTED FROM THE SAR TRANSECTS

<table><thead><tr><th>Transect</th><th>a<sub>3/2</sub></th><th>a<sub>3/2, 3/2</sub></th></tr></thead><tbody><tr><td>T1</td><td>-8.22e-06</td><td>5.39e-06</td></tr><tr><td>T2</td><td>-2.44e-05</td><td>-3.2e-05</td></tr><tr><td>T3</td><td>-1.42e-05</td><td>1.09e-05</td></tr><tr><td>T4</td><td>-2.12e-05</td><td>-3.68e-05</td></tr><tr><td>T5</td><td>-5.3e-05</td><td>-5.1e-05</td></tr></tbody></table>

$a_{3/2}$ and $a_{3/2, 3/2}$ values are computed as averaging values between the $(a_{3/2}, 1, a_{3/2}, 2)$ values and $(a_{3/2}, 1, a_{3/2}, 2)$ values. Deviation from the mean is the same magnitude as $1e-05$.

provide $a_{3/2, 3/2}$ values encompassing the realistic range chosen for the CSM.

The pycnocline depth solutions calculated through ISM are heterogeneous, depending on whether they are deduced from $a_{3/2}$ or $a_{3/2, 3/2}$. All pycnocline depths deduced from $a_{3/2}$ are very deep. However, the pycnocline depths deduced from $a_{3/2, 3/2}$ range
<table><caption>TABLE VII
$h_1$ ESTIMATE FROM a COEFFICIENTS FROM THE SAR TRANSECTS AND ISM</caption><thead><tr><th>Transect</th><th>$h_1$ estimation<br>from $a_{\frac{3}{2}}$</th><th>$h_1$ estimation<br>from $a_{\frac{3}{2}, \frac{3}{2}}$</th></tr></thead><tbody><tr><td>T1</td><td>323 m</td><td>74.5 m</td></tr><tr><td>T2</td><td>163 m</td><td></td></tr><tr><td>T3</td><td>222 m</td><td>58.9 m</td></tr><tr><td>T4</td><td>173 m</td><td>36 m</td></tr><tr><td>T5</td><td>97 m</td><td>31 m</td></tr></tbody></table>

from -31 (at the southern transect) to -74.5 m (at the northern transect). These values are coherent with profiles in Fig. 2 where a strong stratification can be highlighted at these depths. In particular, the CTD stations acquired on 30 September 2008 at 9:43 and 12:22, which are the closest stations from the SAR transects, indicate a pycnocline depth lying, respectively, near 30 and 65 m. They were acquired close to the moment of the propagation of an LAIW [39]. Estimating a pycnocline depth between 30 and 75 m a few hours later seems coherent. In Fig. 2, the CTD profiles from the lower panel were acquired after the LAIW propagation, and, thus, the effects of LAIW isopycnal excursion are less important (close to the rest situation). These CTD profiles (lower panel) were acquired northward and exhibit a pycnocline depth (first maximum of $N^2$) deeper and deeper especially for the stations from 18:29 to 22:22. This indicates a pycnocline depth deeper north than south, in coherence with the calculated $h_1$ (last column in Table VII). To estimate the pycnocline depth through the CSM, the $a_{\frac{3}{2}}$ and $a_{\frac{3}{2}, \frac{3}{2}}$ values (Table VI) are reported in Fig. 6 as dashed white lines. They are very small and constrain the solution into a deep or thick pycnocline layer solutions. In particular, the $a_{\frac{3}{2}}$ values calculated from transects are out of the range chosen for Fig. 6, except for T5 (the pycnocline lying below 80 m and being thicker than 30 m). This is coherent with the values of the pycnocline depth $h_1$ computed through ISM. The $a_{\frac{3}{2}, \frac{3}{2}}$ values imply thin pycnocline lying below 90 m for T5 and below 65 m for T4, whereas it implies a pycnocline lying above 50 m considering a thick pycnocline ($t_p > 30$ m). Here again, the more southern the extracted transect was, the shallower the estimated pycnocline is.

Gathering the ISM and CSM solutions leads to a spectrum of different solutions from a thin and very deep pycnocline (ISM: -323 m < $h_1$ < -97 m and CSM: $z_p$ < -70 m with $t_p$ < 100 m) to a thick and shallower pycnocline (ISM: -74 m < $h_1$ < -31 m and CSM: $z_p$ > -50 m with $t_p$ > 30 m).

## IV. DISCUSSION

Only coefficients $a_{\frac{3}{2}}$, lead to satisfying solutions for the pycnocline depth $z_p$ in the soliton experiment, whereas only $a_{\frac{3}{2}, \frac{3}{2}}$ coefficients lead to realistic solutions according to the CTD data. Both the $\tilde{a}_{\frac{3}{2}, \frac{3}{2}}$ and $\tilde{a}_{\frac{3}{2}}$ coefficients should theoretically lead to the same results. The error can be linked to the linear approximation $(\frac{\eta}{h})^2 \ll 1$. This approximation has consequences: 1) on the $\alpha_2$ value (which is part of $a_{\frac{3}{2}, \frac{3}{2}}$, coefficients and which actually depends on $F(z)$ as used in [44]) and 2) on (12) and thus in all the definitions (20c) and (20d).

On the other hand, the departure between the solutions provided by $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$, can be explained by the nonvalidity of the “two-layer” ocean assumption. The more realistic this approximation is, the more relevant is the choice to use only the first vertical mode (assumed in this study). However, the realistic profiles in Fig. 2 show often several stratification maxima. In the suggested two-layer ocean model, these profiles are modeled through a very thick pycnocline gathering all these maxima. When the real stratification profiles differ from this assumption, the higher vertical modes (second, third, etc.) become less negligible. Yet, this assumption is still used in many studies since it is the simplest ocean model configuration to estimate the pycnocline depth from the LAIW SAR surface signature.

Considering the entire transect provides more information, but it also leads to more errors. Since the error can be amplified through nonlinearity in backscattering, considering the head soliton first and then the tail packet solitons separately could be an alternative for a future work. Moreover, considering the whole transect assumes that the pycnocline depth and thickness are constant, which is not necessarily exact. In Fig. 7, the profiles have a shape close to the dnoidal function. A dnoidal experiment could be an alternative experiment to the LAIW soliton experiment described in Section III-A. In the same way, the theoretical dnoidal function coefficients $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$, could be calculated for the same range of pycnocline depth and thickness and then their values could be reported in Fig. 6 and compared with the LAIW soliton experiment.

One of the most limiting constraints of the CSM method is the computing cost to expand the $\Phi(z)$ to a $R$ order quite high when the pycnocline depth is close to the surface. This study was conducted on a local computer. Using supercomputers could make the study more precise especially in order to evaluate a thicker and shallower pycnocline. Besides, the size of the range for ($z_p, t_p$) values can considerably increase the computing cost. Rather than calculating the coefficients $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$, for the whole ($z_p, t_p$) values, the ($z_p, t_p$) optimal values can be obtained through minimizing the difference between 1) the values of $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$, calculated from the surface signature on the one hand, and 2) the values of $a_{\frac{3}{2}}$, and $a_{\frac{3}{2}, \frac{3}{2}}$, calculated from the stratification model on the other hand.

Solutions for the pycnocline depth were found through ISM or CSM, whereas neither ISM nor CSM is really better to estimate the subsurface stratification. On the opposite, they must be jointly used as they provide complementary information.

Finally, all along this study, the aim was to evaluate two models for pycnocline depth estimation, considering the pycnocline as a limit. Even if the image of a boundary is convenient to understand the dynamics and the exchanges between ocean and atmosphere, the pycnocline is rather a layer where stratification is strong. Defining stratification thresholds, the pycnocline depth, and thickness that could limit the pycnocline extension, is here again a convenient image, the position of the pycnocline being qualitative.

This study was focused on the region of interest of the Alboran Sea especially because of the simultaneity of the LAIW surface signature captured on the SAR image and measurements campaign. However, it would be interesting to transpose the study
on other similar regions: where LAIWs are generated and where the ocean stratification can be described through the stratification model explained above.

## V. CONCLUSION

In this article, we have presented a parametric method to estimate the eKdV parameters, and then to invert these parameters in order to estimate the pycnocline depth. Our approach is derived from the Yule-Walker equations extended to nonlinear autoregressive models. The eKdV parameters have been derived from a simple two-layer interfacial model and a continuously stratified model. Our calculation has shown a departure of the eKdV coefficients in particular for the dissipative effect coefficient (and thus on the SAR signature of the internal waves). Our inversion method has been tested on a simulated soliton shape experiment as well as on SAR images acquired simultaneously to CTD measurements. Some pycnocline depth estimates show a fairly good agreement with the in-situ measurements, while some other parameters are not consistent. In particular, we retrieve the pycnocline depth variation from south to north of the Alboran Sea. A possible explanation of the estimate discrepancy is that the two-layer model is questionable as observed in the in-situ measurements.

## ACKNOWLEDGMENT

The authors would like to thank the reviewers for all useful and helpful comments on this manuscript. The authors would also like to thank the research team from the Departamento de Física Aplicada from the Universidad de Cadiz for providing in-situ CTD density profiles. Finally, M. Dessert would like to thank the research and technical teams from EXWEXs for their help and implication during this study.

## REFERENCES

[1] P. Brandt, W. Alpers, and J. O. Backhaus, "Study of the generation and propagation of internal waves in the Strait of Gibraltar using a numerical model and synthetic aperture radar images of the European ERS 1 satellite," *J. Geophysical Res.: Oceans*, vol. 101, no. C6, pp. 14237-14252, 1996.

[2] C. Jackson, "Internal wave detection using the Moderate Resolution Imaging Spectroradiometer (MODIS)," *J. Geophysical Research: Oceans*, vol. 112, no. C11, pp. 1-10, 2007. [Online]. Available: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2007JC004220

[3] J. C. B. da Silva, A. L. New, M. A. Srokosz, and T. J. Smyth, "On the observability of internal tidal waves in remotely-sensed ocean colour data," *Geophysical Res. Lett.*, vol. 29, no. 12, pp. 10-1-10-4, 2002. [Online]. Available: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2001GL013888

[4] X. Pan, G. T. F. Wong, F.-K. Shiah, and T.-Y. Ho, "Enhancement of biological productivity by internal waves: Observations in the summertime in the northern south China sea," *J. Oceanogr.*, vol. 68, no. 3, pp. 427-437, Jun. 2012.

[5] H. Kim, Y. B. Son, and Y.-H. Jo, "Hourly observed internal waves by geostationary ocean color imagery in the East/Japan sea," *J. Atmospheric Ocean. Technol.*, vol. 35, no. 3, pp. 609-617, Mar. 2018.

[6] A. K. Liu, Y. S. Chang, M.-K. Hsu, and N. K. Liang, "Evolution of nonlinear internal waves in the East and South China Seas," *J. Geophysical Res.: Oceans*, vol. 103, no. C4, pp. 7995-8008, 1998.

[7] D. L. Porter and D. R. Thompson, "Continental shelf parameters inferred from SAR internal wave observations," *J. Atmospheric Ocean. Technol.*, vol. 16, no. 4, pp. 475-487, Apr. 1999.

[8] X. Li, P. Clemente-Colón, and K. S. Friedman, "Estimating oceanic mixed-layer depth from internal wave evolution observed from Radarsat-1 SAR," *Johns Hopkins APL Tech. Dig.*, vol. 21, no. 1, pp. 130-135, 2000.

[9] Q. Zheng, Y. Yuan, V. Klemas, and X.-H. Yan, "Theoretical expression for an ocean internal soliton synthetic aperture radar image and determination of the soliton characteristic half width," *J. Geophysical Res.: Oceans*, vol. 106, no. C12, pp. 31415-31423, Dec. 2001.

[10] Z. Zhao, V. Klemas, Q. Zheng, X. Li, and X.-H. Yan, "Estimating parameters of a two-layer stratified ocean from polarity conversion of internal solitary waves observed in satellite SAR images," *Remote Sens. Environ.*, vol. 92, no. 2, pp. 276-287, Aug. 2004.

[11] J.-M. Le Caillec, "Study of the SAR signature of internal waves by nonlinear parametric autoregressive models," *IEEE Trans. Geosci. Remote Sens.*, vol. 44, no. 1, pp. 148-158, Jan. 2006.

[12] H. Shen, W. Perrie, and C. L. Johnson, "Predicting internal solitary waves in the gulf of Maine," *J. Geophysical Res.: Oceans*, vol. 125, no. 3, 2020, Art. no. e2019JC015941. [Online]. Available: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2019JC015941

[13] J. Utiz, D. Stramski, B. Gentili, F. D'Ortenzio, and H. Claustre, "Estimates of phytoplankton class-specific and total primary production in the Mediterranean Sea from satellite ocean color observations," *Glob. Biogeochemical Cycles*, vol. 26, no. 2, 2012, Art. no. GB2024. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2011GB004055

[14] S. Pierini, "A model for the Alboran Sea internal solitary waves," *J. Phys. Oceanogr.*, vol. 19, no. 6, pp. 755-772, Jun. 1989.

[15] H. Lacombe and C. Richez, "The Regime of the Strait of Gibraltar," in *Elsevier Oceanography Series, Ser. Hydrodynamics of Semi-Enclosed Seas*, vol. 34, J. C. J. Nihoul, Ed. Amsterdam, Netherlands: Elsevier, Jan. 1982, pp. 13-73.

[16] L. Armi and D. Farmer, "The flow of Mediterranean water through the strait of Gibraltar. The flow of Atlantic water through the Strait of Gibraltar," *Prog. Oceanogr.*, vol. 21, no. 1, pp. 1-105, 1988.

[17] D. M. Farmer and L. Armi, "The flow of Atlantic water through the Strait of Gibraltar," *Prog. Oceanogr.*, vol. 21, no. 1, pp. 1-103, Jan. 1988.

[18] C. R. Jackson, J. da Silva, G. Jeans, W. Alpers, and M. J. Carusso, "Non-linear internal waves in synthetic aperture radar imagery," *Oceanography*, vol. 26, no. 2, pp. 68-79, 2013.

[19] W. J. Plant, "Bragg scattering of electromagnetic waves from the air/sea interface," in *Surface Waves and Fluxes: Volume II – Remote Sensing, ser. Environmental Fluid Mechanics*, G. L. Geernaert and W. L. Plant, Eds. Dordrecht: Springer Netherlands, 1990, pp. 41-108.

[20] T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, "A unified directional spectrum for long and short wind-driven waves," *J. Geophysical Res.: Oceans*, vol. 102, no. C7, pp. 15781-15796, 1997.

[21] V. N. Kudryavtsev, V. K. Makin, and B. Chapron, "Coupled sea surface-atmosphere model: 2. Spectrum of short wind waves," *J. Geophysical Res.: Oceans*, vol. 104, no. C4, pp. 7625-7639, 1999.

[22] A. Bringer, B. Chapron, A. Mouche, and C.-A. Guérin, "Revisiting the short-wave spectrum of the sea surface in the light of the weighted curvature approximation," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 1, pp. 679-689, Jan. 2014.

[23] R. Romeiser, A. Schmidt, and W. Alpers, "A three-scale composite surface model for the ocean wave-radar modulation transfer function," *J. Geophysical Res.: Oceans*, vol. 99, no. C5, pp. 9785-9801, 1994.

[24] R. Romeiser and W. Alpers, "An improved composite surface model for the radar backscattering cross section of the ocean surface: 2. Model response to surface roughness variations and the radar imaging of underwater bottom topography," *J. Geophysical Res.: Oceans*, vol. 102, no. C11, pp. 25251-25267, 1997.

[25] X. Zhang, Z.-S. Wu, and X. Su, "Electromagnetic scattering from deterministic Sea surface with oceanic internal waves via the variable-coefficient gardener model," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 2, pp. 355-366, Feb. 2018.

[26] T. Song, C. Tong, and L. Cong, "Numerical analysis of influences from internal waves on electromagnetic scattering from sea surface," *IEEE Access*, vol. 6, pp. 42014-42021, 2018.

[27] W. Alpers and I. Hennings, "A theory of the imaging mechanism of underwater bottom topography by real and synthetic aperture radar," *J. Geophysical Res.: Oceans*, vol. 89, no. C6, pp. 10529-10546, 1984.

[28] D. D. J. Korteweg and D. G. d. de Vries, "On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves," *London, Edinburgh, Dublin Philos. Mag. J. Sci.*, vol. 39, no. 240, pp. 422-443, May 1895.

[29] T. B. Benjamin, "Internal waves of permanent form in fluids of great depth," *J. Fluid Mechanics*, vol. 29, no. 3, pp. 559-592, Sep. 1967.

[30] H. Ono, "Algebraic solitary waves in stratified fluids," *J. Phys. Soc. Jpn.*, vol. 39, no. 4, pp. 1082-1091, Oct. 1975.

[31] R. I. Joseph, "Solitary waves in a finite depth fluid," *J. Phys. A: Math. Gen.*, vol. 10, no. 12, pp. L225-L227, Dec. 1977.
[32] T. Kubota D. R. S. Ko, and L. D. Dobbs, “Weakly nonlinear internal waves in a stratified fluid of finite depth,” *J. Hydrodyn.*, vol. 12, no. 4, pp. 157–165, 1978.

[33] E. N. Pelinovski, O. Poloukhina, A. V. Slyunyaev, and T. G. Talipova, “Internal solitary waves,” in *Solitary Waves in Fluids*, ser. WIT Transactions on State of the Art in Science and Engineering, vol. 9. Ashurst, England: WIT Press, 2007, pp. 85–110.

[34] O. Kurkina, T. Talipova, E. Pelinovsky, and T. Somere, “Mapping the internal wave field in the Baltic Sea in the context of sediment transport in shallow water,” *J. Coastal Res.*, vol. SI 64, pp. 2042–2047, 2011.

[35] U. Libal and K. H. Johansson, “Yule-walker equations using higher order statistics for nonlinear autoregressive model,” in *Proc. Signal Process. Symp.*, 2019, pp. 227–231.

[36] K. R. Helfrich and W. K. Melville, “Long nonlinear internal waves,” *Annu. Rev. Fluid Mechanics*, vol. 38, no. 1, pp. 395–425, 2006.

[37] V. I. Vlasenko, “Multimodal soliton of internal waves,” *Izv. Akad. Nauk SSSR. Fiz. Atmos. i Okeana*, vol. 30, no. 2, pp. 173–181, 1994.

[38] S. R. Massel, “Waves in a continuously stratified sea of constant water depth,” in *Internal Gravity Waves in the Shallow Seas* (GeoPlanet: Earth and Planetary Sciences Series), S. R. Massel, Ed. Cham, Switzerland: Springer, 2015, pp. 45–86.

[39] Á. Vázquez, S. Flecha, M. Bruno, D. Macías, and G. Navarro, “Internal waves and short-scale distribution patterns of chlorophyll in the Strait of Gibraltar and Alboran Sea,” *Geophysical Res. Lett.*, vol. 36, no. 23, 2006, Art. no. L23601. [Online]. Available: http://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2009GL040959

[40] T. V. La, C. Messager, M. Honnorat, and C. Channelliere, “Detection of convective systems through surface wind gust estimation based on Sentinel-1 images: A new approach,” *Atmospheric Sci. Lett.*, vol. 19, no. 12, 2018, Art. no. e863. [Online]. Available: https://onlinelibrary.wiley.com/doi/pdf/10.1002/asl.863

[41] T. V. La and C. Messager, “Convective system sea surface wind pattern detection and variability observation from a combination of Sentinel-1 and Radarsat-2 images,” *Remote Sens. Lett.*, vol. 11, no. 5, pp. 446–454, Feb. 2020.

[42] T. V. La et al., “Use of Sentinel-1 C-band SAR images for convective system surface wind pattern detection,” *J. Appl. Meteorol. Climatol.*, vol. 59, no. 8, pp. 1321–1332, Aug. 2020.

[43] H. Hersbach, A. Stoffelen and S. de Haan, “An improved C-band scatterometer ocean geophysical model function: CMOD5,” *J. Geophys. Res.*, Vol. 112, 2007, Art. no. C03006.

[44] O. Kurkina, E. Rouvinskaya, A. Kurkin, A. Giniyatullin, and E. Pelinovsky, “Vertical structure of the velocity field induced by mode-I and mode-II solitary waves in a stratified fluid,” *Eur. Phys. J. E*, vol. 41, no. 3, pp. 1–8, Mar. 2018.

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

**Morgane Dessert** is currently working toward the Ph.D. degree with EXWEXs and Laboratory for Ocean Physics and Satellite Remote Sensing (LOPS), University of Brest, Brest, France.

She worked as an Engineer with the European Institute for Marine Studies (IUEM), Plouzané, France, from 2012 to 2017, and has been with EXWEXs since 2018.

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

**Marc Honnorat** received the M.S. degree in digital image and signal processing from Cranfield University, Cranfield, U.K., in 2002, and the Ph.D. degree in applied mathematics from Grenoble INP, Grenoble-Alpes University, Grenoble, France, in 2007.

From 2007 to 2014, he worked as a Research Engineer with the French National Institute for Research in Digital Science and Technology (Inria), Le Chesnay-Rocquencourt, France. Since 2016, he has been a Research Engineer with EXWEXs, Brest, France. He has experience in applied mathematics,

numerical ocean and weather modeling, data assimilation, and machine learning.

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

**Jean-Marc Le Caillec** (Senior Member, IEEE) received the Engineering degree in telecommunications from Telecom Bretagne (IMT Atlantique), Brest, France, in 1992, and the Ph.D. degree in mathematics and signal processing from the University of Rennes I, Rennes, France, in 1992.

He is a Professor with IMT Atlantique, Nantes, France, in the 2IP (Information and Image Processing) department. From 1997 to 1999, he worked with Thomson AirSys (now Thales AirSys), La Défense, France. He joined Telecom Bretagne, Brest, France,

as an Associate Professor in 1999 and became a Professor in 2007. From 2014 to 2020, he was in charge of the division CID (Knowledge, Information, Decision) of the Lab-STICC (Laboratory of sciences and technics for information, communication, and knowledge, UMR 6285). His main research interests include statistics, nonlinear system modeling, mathematics, and signal processing for applications in remote sensing and finance.

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

**Christophe Messager** received the master's degree in meteorology and oceanography from the Université de Toulon et du Var, La Valette-du-Var, France, in 1994, the M.Sc. degree in physical oceanography from the Centre d'Océanologie de Marseille, Marseille, France, in 1995, the Ph.D. degree in atmospheric science from the Université Joseph Fourier, Saint-Martin-d'Héres, France, in 2005, and then the “Habilitation à Diriger de Recherche-Accreditation to Supervise Research,” Université de Bretagne Occidentale, Brest, France, in 2015.

He was the former Founder and the Head of the ICEMASA International French-South African Marine Laboratory, then the Head of the Laboratoire de Physique des Océans-Ocean Physics Laboratory until 2015. From spring 2012, he was the Founder of Extreme Weather Expertises company. In 2020, he founded the B-Space company for which he is also the leader. His research interests include atmospheric science, meteorology and climate, and especially convective systems and ocean-atmosphere interactions, as well as the applications of remote sensing in satellite meteorology and oceanography.

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

**Xavier Carton** received the Ph.D. degree in physical oceanography from the University of Paris, Paris, France, in 1988, and the Sc.D. degree in geostrophic vortex stability and interactions from the University of Brest, Brest, France, in 1999.

He is the Professor of fluid dynamics and applied mathematics. He worked for the French Navy and for IFREMER. Since 2004, he has been a Professor with the University of Brest. He is specialized in flow instability, vortex dynamics, and outflows from marginal seas.