Estimating the Pycnocline Depth From the SAR Signature of Internal Waves in the Alboran Sea

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 coefficients were 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.


I. INTRODUCTION
O CEANS 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 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).
Digital Object Identifier 10.1109/JSTARS.2022.3214298 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. 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-temperaturedensity (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 wellknown 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. 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.

A. Relation Between NRCS and the LAIW-Induced Surface Currents
SAR images provide the normalised radar cross section (NRCS) denoted σ 0 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 −1 ) and lower than roughly 10 ms −1 [18]. Besides, the SAR incidence angle θ i 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, σ 0 can be expressed as in [19] where k B is the Bragg wavevector (in what follows, the vectors are denoted as k while their corresponding norm as k). The Bragg wavenumber k B , pulsation ω B , and period T B can be computed as 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 and ε r is the relative dielectric constant of seawater. Ψ( 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: where 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 σ 0,LAIW and surface current modulation as whereσ 0,LAIW (x r , x a ) = (σ 0 − σ 0 )/σ 0 is the LAIW-induced NRCS anomaly where σ 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 ≈ 130 s for SEASAT satellite [27]). Finally, τ r is the relaxation time conveying the source functions of (4). Physically, τ r is the response time of the wave system to the current variation and was set to τ r ≈ 30 − 40 s [27]. Then, as in [27], (5) is rewritten in the ( x P , x T ) referential where x P is the direction along the LAIW propagation direction, whereas x T is the direction across the LAIW propagation direction. Thus, the rotation and the derivative introduce cos 2 (φ) and cos(φ) sin(φ) factors with φ the angle between the satellite range direction x r and the propagation direction x P of the LAIW. Following the same approach, ∂U 0 ∂x P is assumed larger than ∂U 0 ∂x T , implying that the σ 0,LAIW variations are negligible along x T . Equation (5) becomes From (6), the surface current variation ∂U 0,LAIW ∂x P can be deduced from the NRCS σ 0 . The next sections describe the relation between the LAIW-induced current field U 0,LAIW (x P ) and the density profile ρ(z) in order to connect the NRCS and the density variation due to LAIWs.

1) Definitions:
The following three main models are used for modeling LAIWs depending on the ratio between the horizontal scale λ and the vertical scale H, which is the sea bottom depth: 1) the Korteweg-deVriès model when λ > H [28]; 2) the Benjamin-Ono model when λ ≈ H [29], [30]; 3) the Joseph-Kubota model when λ < 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 λ = 2 km, which is in accordance with typical horizontal length of the Alboran Sea [14]. Under these conditions (λ > 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 where α conveys the nonlinear effects, β 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 ≈ 20 − 50 m and the magnitude of the LAIW-induced pycnocline elevation was A ≈ 1 − 10 m. When α 2 ≈ β (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 α ≈ β. In the eKdV framework, the space-time variation of the LAIW isopycnal elevation field ξ(x, z, t) is given by where c is the phase velocity under linear approximation, α 1 and α 2 are nonlinearity coefficients at first and second orders, β 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 ξ 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 ξ can be expressed as a series as in [33] ξ where η(x, t), the isopycnal elevation, still satisfies the eKdV (8) and where N (z) is the Brunt-Väisälä pulsation and g is the gravity acceleration. The function Φ(z) is the solution of Taylor-Goldstein equation [see (10a)], with F (z) being the first-order nonlinear correction of the Φ(z) function. From the definition of the isopycnal deformation [see (9)], the LAIW-induced surface current can be expressed as [34] Assuming that ( η h ) 2 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 where IP * stands for (nonnormalized) integrated profile, 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 λ and a vertical length H such as δ x P λ 1 (where δ x p is the spatial sampling step along the LAIW propagation direction as discussed again in Section III-B1) and since the LAIW is considered to be propagating as a mature wave, i.e., ∂η ∂t V being the LAIW velocity. Now, the eKdV can be expressed only through spatial derivatives (13) leads to a derived eKdV equation linking the NRCS σ 0 and the eKdV coefficients as 1 being the normalized integrated profile. Up to the first order, the derivatives can be approximated as 1 is the normalized sampling spatial step (along the LAIW propagation direction) and, thus, (14) becomes a weighed summation of linear and nonlinear (quadratic and cubic) products of IP(s) at several lags of s where a n are the weighting coefficients and M and N are the summation upper bounds and d (s) = o(δ s ) 1 is the independent error function having a null statistical mean (i.e., d = 0). In what follows, the bar stands for statistical mean along the LAIW propagation direction. In order to reduce the error 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 IP(s) [right-hand term in (17)] requires N > M.
Then, expression (17) is multiplied successively by The noninteger lags (as s − (N − 1)δ s /2) are obtained by averaging two adjacent samples in this matrix.
Finally, the average operator · is applied to each equation leading to a linear system of N + 2 M equations. For writing convenience, the function h(i, j, . . . , k) is introduced: (18) As a remark, the indices are permutative: h i,j,...,k = h k,...,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 On the one hand, coefficients a can be calculated from 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 Φ (z), c, α 1 , and β 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 ρ up is the upper layer density, while ρ bot is the lower layer density 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 Φ (0) → − 1 h 1 , thus leading to coefficients where h 2 = H − h 1 is the thickness of the bottom layer, is the density jump, related to the stiffness of the pycnocline and ρ = ρ up + ρ bot 2 is the average density over the water column (see [36]). Combining (22b), (22d), and (20c) leads to a relation between the a 3 2 ,· and a 3 2 , 3 2 ,· coefficients and the pycnocline depth h 1 with α 1 , α 2 , β 1 , and c nonlinearly depending on h 1 . Combining (23) and (19) constitutes the interfacial stratification model (ISM) linking the NRCS σ 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, Φ(z) in (10a) must be computed by solving an eigenvalue Sturm-Liouvillle problem with boundary conditions Φ(0) = Φ(H) = 0. This system has an infinite number of solutions defined by the pairs of eigenvalues ϕ (k) and eigenfunctions (also called vertical modes) W (k) (z) as where K is the order up to which the function Φ(z) is described (i.e., the number of vertical modes). Each vertical mode W (k) (z) is described as where R (k) is the order up to which the kth mode is described.
The kth vertical mode has a propagation velocity c (k) that can be computed from the kth eigenvalue ϕ (k) through 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 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 and when the quantities 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 Φ and the depth of the middle of the pycnocline are relatively close [see (27)] and when a 3 2 ,· and a 3 2 , 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 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 γ r ), then a continuous density profile, clearly separating two layers, is introduced as in [ where t p expresses the pycnocline thickness. This approach involves introducing N (z) into (10a) and retrieving γ r and then calculating the pair (ϕ, 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] where α 1 , α 2 , Φ (0), and c can then be calculated from (z p , t p ); but as seen in Section III-A3, the two pycnocline parameters cannot be uniquely estimated from these coefficients. Combining (30a), (30b), (30c), and Φ(z) and c calculated from (25) and (26) leads to The a 3 2 ,· and a 3 2 , 3 2 ,· coefficients are expressed through α 1 , α 2 , β, c, and Φ (0) which all depend on the stratification parameters (z p , t p ). Combining (31) and (19) constitutes the CSM binding the NRCS σ 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 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 referencial 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 ≈ 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 Gardnertype soliton propagating in a two-layer ocean.
1) Stratification Simulation: Fig. 4 illustrated at left panel the density profile ρ(z) [defined as in (29)]. At upper right panel, where η 0 is the maximum amplitude. The maximum amplitude is set to η 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 Φ(z) is expanded up to the optimal order R opt . In Fig. 5,  (27) and (28)]. Moreover, the thicker the pycnocline is, the higher the R must be to satisfy 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:  In particular, β 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: ,.
Δx 2 (for upper panel) Δ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).  (z p , t p ).
The value of the coefficient a 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 3

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 δ x p is given by the SAR image resolution divided by cos(φ) where φ 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 θ i (satellite characteristics) of the SAR is around 24 • . The averaged bathymetry (from GEBCO) H under the transects is H = 792 ± 46 m with low averaged slopes ∂H/∂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 τ r is expected to be 10-100 Bragg wave periods. In 1984, Alpers and Hennings [27]   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 p 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   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 3 2 and a 3 2 , 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 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 3 2 , 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

IV. DISCUSSION
Only coefficients a 3 2 ,· lead to satisfying solutions for the pycnocline depth z p in the soliton experiment, whereas only a 3 2 , 3 2 ,· coefficients lead to realistic solutions according to the CTD data. Both theã 3 2 , 3 2 ,· andã 3 2 ,· coefficients should theoretically lead to the same results. The error can be linked to the linear approximation η h 2 1. This approximation has consequences: 1) on the α 2 value (which is part of a 3 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 3 2 ,· and a 3 2 , 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 3 2 ,· and a 3 2 , 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 Φ(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 3 2 ,· and a 3 2 , 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 3 2 ,· and a 3 2 , 3 2 ,· calculated from the surface signature on the one hand, and 2) the values of a 3 2 ,· and a 3 2 , 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.