Retrieval of Internal Solitary Wave Amplitude in Shallow Water by Tandem Spaceborne SAR

: The accurate estimation of the upper layer thickness in a two-layer ocean is a crucial step in the retrieval of internal solitary wave (ISW) amplitude from synthetic aperture radar (SAR) data. In this paper, we present a method to derive the upper layer thickness and the consequent ISW amplitude by combining two consecutive SAR images with the extended Korteweg-de Vries (eKdV) equation. An ISW case observed twice by the Chinese C-band SAR GaoFen-3 (GF-3) and the German X-band SAR TerraSAR-X (TS-X) with a temporal interval of approximately 11 min in shallow water to the southeast of Hainan Island in the northwestern South China Sea was used to demonstrate the applicability of the method. Using the in situ measurements of temperature and salinity near the observed ISW, the proposed method yielded an ISW amplitude of − 4.52 m, in close proximity to − 5.66 ± 1.24 m derived by applying the classic Korteweg–de Vries (KdV) equation based on the continuously stratiﬁed theory. Moreover, the climatological dataset of the World Ocean Atlas 2013 (WOA13) was also used with the proposed method in the Hainan case, and the results showed that the method can still provide a reasonable estimate of ISW amplitude in shallow water even when in situ oceanic stratiﬁcation measurements are absent. The application of our method to derive the ISW amplitude from consecutive SAR images seems highly promising with the increasing emergence of tandem satellites in orbits.


Introduction
Internal solitary waves (ISWs) are a class of high-frequency, non-linear, and non-hydrostatic gravity waves that are widely observed in the coastal oceans and marginal seas [1][2][3][4]. Internal solitary waves can enhance ocean mixing and are an important link in energy dissipation from the barotropic tide [5,6]. Additionally, ISWs can have a profound effect on offshore drilling operations [7], nutrients, pollutants, and sediments transport [8], acoustic propagation [9], and nearshore ecosystems [10]. For these reasons, the generation, propagation, steepening, and dissipation of ISWs have received much research attention for several decades.
Spaceborne synthetic aperture radar (SAR) is a powerful remote sensing instrument for observing ISWs [11]. Internal solitary waves can induce variations in surface currents and locally modulate the distribution and intensity of Bragg waves, which renders the ISWs visible on SAR images as The accurate estimation of upper layer thickness is crucial in determining the ISW amplitude using a two-layer water fluid model because the retrieved amplitude is sensitive to the upper layer thickness [21,22]. Although a variety of methods have been proposed to estimate the upper layer thickness, such as the maximum buoyancy frequency (N max ) method, eigenfunction method, empirical orthogonal function (EOF) method, KdV method, Benjamin-Ono (BO) method, KdV-I method, and the KdV-II method; the estimations using these methods may have large discrepancies [32] and these methods based on in situ measurements are costly and can provide only point estimations. Therefore, some studies have attempted to infer upper layer thickness from SAR images as they can provide a wide view of ISW evolution with high spatial resolutions around tens of meters. For example, Zhao et al. [33] proposed that the upper layer thickness can be calculated by determining the polarity conversion point of the observed ISWs in SAR images, where the depression ISWs are converted into elevation ISWs and the upper layer thickness is exactly half that of the local water depth. However, this method is limited to the conversion point estimation. Li et al. [34] and Yang et al. [35] also derived the upper layer thickness from SAR images under the assumption that the ISW phase speed is close to the group velocity of ISW packets. Thus, the derived phase speed is an average speed within the whole tide period instead of the "true" phase speed of ISWs in a certain location, which may lead to large uncertainties in ISW amplitude estimations. In such cases, the consecutive remote sensing image pairs within a short time (tandem images) have potential applications in estimating the accurate upper layer thickness and subsequent amplitudes because ISWs can be observed multiple times in the tandem images, and therefore, we can derive a more accurate phase speed [36].
Therefore, one can find that although the two-layer model is more practical in retrieving ISW amplitude from SAR images than the continuously stratified model does, the facing challenge is the accurate estimation of upper layer thickness. In this study, we proposed a method of deriving the optimum upper layer thickness based on the accurate estimation of ISW phase speed using a tandem SAR image pair acquired within a short temporal interval. Consequently, the amplitude of the ISW is derived. An ISW case repeatedly observed by the Chinese C-and SAR GaoFen-3 (GF-3) and the German X-band SAR TerraSAR-X (TS-X) with a short temporal interval in shallow water on southeastern Hainan Island is used to demonstrate the method.

Brief Description of the Experiment
An ISW observational experiment was conducted on the southeast of Hainan Island in the northwestern South China Sea from 10 to 12 June 2017. During the experiment, a GF-3 and a TS-X image with a short temporal interval of approximately 11 min were acquired, whose spatial coverage are illuminated in Figure 1. The two SAR data across the shelf break cover most sea regions between the southeast of Hainan Island and the Xisha Island, and they both clearly present signals of ISWs which are observed to the west of the shelf break where the water depth is less than −200 m. More details about the SAR data are given in the next subsection. Two conductivity-temperature-pressure (CTD) casts of S1 (109.8190 • E, 18.1596 • N) and S2 (109.9268 • E, 18.1167 • N) located near the observed ISWs were performed during the experiment and their locations are indicated by red dots in Figure 1. The two casts were used to collect the background stratification (density) profiles of the study area, which are shown later in the paper. German X-band SAR TerraSAR-X (TS-X) with a short temporal interval in shallow water on southeastern Hainan Island is used to demonstrate the method.  The two locations S1 and S2 where conductivity-temperature-pressure (CTD) measurements were taken are marked by red dots. Table 1 lists the technical specifications of the TS-X and GF-3 images acquired in the experiment. The TS-X SAR image was acquired in the scanSAR mode with a single-polarization of vertical-vertical (VV). The GF-3 image was acquired in the standard stripmap mode with a dual-polarization of VV and vertical-horizontal (VH). Portions of the GF-3 (VV polarization) and TS-X image in the same geographical area are shown in Figure 2a,b, respectively. The manifestations of ISWs patterns are clearer in the GF-3 image than that in the TS-X image, as the former had better spatial resolution of 8 m. On the other hand, the incidence angle of the GF-3 sub-image was approximately 15.23-21.70 degrees, while that of the TS-X sub-image was 42.50-48.72 degrees. Both images were processed by radiometric calibration, speckle filtering, and geolocation correction. In the two images, one can clearly observe three ISW packets distributed southeast of Hainan Island. The spatial-temporally collocated ERA5 reanalysis wind model suggests that the sea surface wind speed was no more than 5 m/s at the acquisition time of the SAR data, which favors a higher probability for observing these ISWs on SAR images during summer [17]. These observed ISWs in Figure 2 are interpreted as mode-1 depression ISWs because their SAR signatures appear as bright-dark stripes. The convex curvatures imply that these ISWs were propagating northwestwards onto the Hainan shelf. The most prominent packet in these observed ISWs is labeled P1 in Figure 2a. A transect (AA', the red dashed line in Figure 2a) parallel to the wave propagation direction was selected to analyze the variation in the normalized radar cross-section (NRCS) of the ISWs. To suppress the noise, the NRCS Remote Sens. 2019, 11, 1706 5 of 18 along AA' was averaged in the direction perpendicular to the wave propagation direction, and then, the averaged NRCS was smoothed by a moving average filter and plotted in Figure 3. The horizontal distribution curve of pixel backscattering intensity values (averaged and smoothed NRCS) clearly reveals the leading three distinct solitons of packet P1. In this paper, we chose the first leading soliton, which is named L1, to derive the amplitude. From Figure 3, we obtain the wavelength of soliton L1 as approximately 440 m, which is defined as the distance between the crest of L1 and its following soliton. The local nautical chart shows that the water depth near L1 was approximately −74 m. Therefore, the wavelength of L1 in this case was greater than the water depth, which suggests that the ISW satisfies the shallow water internal wave criteria.

Processing of the Spaceborne SAR Data and Extraction of ISW Characteristics
patterns are clearer in the GF-3 image than that in the TS-X image, as the former had better spatial resolution of 8 m. On the other hand, the incidence angle of the GF-3 sub-image was approximately 15.23-21.70 degrees, while that of the TS-X sub-image was 42.50-48.72 degrees.  Both images were processed by radiometric calibration, speckle filtering, and geolocation correction. In the two images, one can clearly observe three ISW packets distributed southeast of Hainan Island. The spatial-temporally collocated ERA5 reanalysis wind model suggests that the sea surface wind speed was no more than 5 m/s at the acquisition time of the SAR data, which favors a higher probability for observing these ISWs on SAR images during summer [17]. These observed ISWs in Figure 2 are interpreted as mode-1 depression ISWs because their SAR signatures appear as bright-dark stripes. The convex curvatures imply that these ISWs were propagating northwestwards onto the Hainan shelf. The most prominent packet in these observed ISWs is labeled 1 in Figure  2(a). A transect (AA', the red dashed line in Figure 2(a)) parallel to the wave propagation direction was selected to analyze the variation in the normalized radar cross-section (NRCS) of the ISWs. To suppress the noise, the NRCS along AA' was averaged in the direction perpendicular to the wave propagation direction, and then, the averaged NRCS was smoothed by a moving average filter and plotted in Figure 3. The horizontal distribution curve of pixel backscattering intensity values (averaged and smoothed NRCS) clearly reveals the leading three distinct solitons of packet 1. In this paper, we chose the first leading soliton, which is named 1, to derive the amplitude. From Figure 3, we obtain the wavelength of soliton 1 as approximately 440 m, which is defined as the distance between the crest of 1 and its following soliton. The local nautical chart shows that the water depth near 1 was approximately -74 m. Therefore, the wavelength of 1 in this case was greater than the water depth, which suggests that the ISW satisfies the shallow water internal wave criteria.

CTD Data and WOA13 V2 Data
The mean water density profiles shown in Figure 4a were sampled by the CTD casts of S1 and S2 on 12 June 2017. The soliton 1 was located approximately 6 km from the both CTD casts and was observed on SAR images about two days before the CTD casts. Thus, the in situ water density profile in Figure 4a is treated as a mean stratification profile of the ocean background for the time in which the ISW case took place in the study. Figure 4b presents the corresponding background buoyancy frequency profile. A pycnocline is the cline or layer where the vertical water density

CTD Data and WOA13 V2 Data
The mean water density profiles shown in Figure 4a were sampled by the CTD casts of S1 and S2 on 12 June 2017. The soliton L1 was located approximately 6 km from the both CTD casts and was observed on SAR images about two days before the CTD casts. Thus, the in situ water density profile in Figure 4a is treated as a mean stratification profile of the ocean background for the time in which the ISW case took place in the study. Figure 4b presents the corresponding background buoyancy frequency profile. A pycnocline is the cline or layer where the vertical water density changes dramatically. Thus, as shown in Figure 4, the study area features a strong and thick seasonal pycnocline located approximately between −9 m and −33 m. The thickness of the seasonal pycnocline accounts for approximately 32% of the total water depth, making the determination of the upper layer thickness in a two-layer model rather ambiguous and difficult.  The stratification profile is needed to be determined as an input to the proposed ISW amplitude estimation method. The conventional approach to underline the density stratification is to use in situ or climatological data. So apart from the in situ CTD data, we also used the background stratification profile from the monthly objectively analyzed mean climatological dataset of WOA13 V2 as an input to the proposed method, in order to examine the practicability of the proposed method in cases where in situ CTD measurements were not available. The WOA13 V2 dataset has a horizontal grid resolution of 0.25° and a 5 m depth vertical interval for water depths of less than 100 m. In this dataset, the point of (109.875° E, 18.125° N) is the nearest grid point to the location of soliton 1. Therefore, the temperature and salinity information from WOA13 V2 at this point in June was chosen to analyze the background stratification. The temperature and salinity profiles were interpolated to 1 m depth intervals and then substituted into the Thermodynamic Equation of Sea Water 2010 (TEOS-10) to calculate the density profile and its corresponding background buoyancy frequency profile, which are shown in Figure 5a,b, respectively. The profiles in Figure 5b are similar to the in situ measurements shown in Figure 4b, and it also shows a strong and thick seasonal pycnocline, which is located approximately between −12 m and −47 m. The buoyancy frequency value from the WOA13 V2 dataset shown in Figure 5b is about 0.15 s −1 smaller than that of our in situ measurements. The difference between WOA13 V2 climatological data and in situ measurements may be due to the impacts of weekly synoptic weather events and seasonally/yearly climate signals on the water column condition in shallow water regions. Given the free availability of WOA13 V2 data and the similarity between the WOA13 V2 and the in situ measured buoyancy frequency profile, we also applied the proposed method to derive amplitude by using the WOA13 V2 data in the study area, which is then compared with the derived amplitude by using the in situ stratification data. The stratification profile is needed to be determined as an input to the proposed ISW amplitude estimation method. The conventional approach to underline the density stratification is to use in situ or climatological data. So apart from the in situ CTD data, we also used the background stratification profile from the monthly objectively analyzed mean climatological dataset of WOA13 V2 as an input to the proposed method, in order to examine the practicability of the proposed method in cases where in situ CTD measurements were not available. The WOA13 V2 dataset has a horizontal grid resolution of 0.25 • and a 5 m depth vertical interval for water depths of less than 100 m. In this dataset, the point of (109.875 • E, 18.125 • N) is the nearest grid point to the location of soliton L1. Therefore, the temperature and salinity information from WOA13 V2 at this point in June was chosen to analyze the background stratification. The temperature and salinity profiles were interpolated to 1 m depth intervals and then substituted into the Thermodynamic Equation of Sea Water 2010 (TEOS-10) to calculate the density profile and its corresponding background buoyancy frequency profile, which are shown in Figure 5a,b, respectively. The profiles in Figure 5b are similar to the in situ measurements shown in Figure 4b, and it also shows a strong and thick seasonal pycnocline, which is located approximately between −12 m and −47 m. The buoyancy frequency value from the WOA13 V2 dataset shown in Figure 5b is about 0.15 s −1 smaller than that of our in situ measurements. The difference between WOA13 V2 climatological data and in situ measurements may be due to the impacts of weekly synoptic weather events and seasonally/yearly climate signals on the water column condition in shallow water regions. Given the free availability of WOA13 V2 data and the similarity between the WOA13 V2 and the in situ measured buoyancy frequency profile, we also applied the proposed method to derive amplitude by using the WOA13 V2 data in the study area, which is then compared with the derived amplitude by using the in situ stratification data. Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 18 (a) (b) Figure 5. Background profiles of (a) density and (b) buoyancy frequency in June obtained from the monthly WOA13 V2 dataset.

The Conventional Method of ISW Amplitude Estimation Based on the Classic KdV Equation in a Continuously Stratified Ocean Model
The conventional ISW amplitude estimation method applying the classic KdV equation in a continuously stratified ocean was developed by Small et al. [18]. Here, we briefly describe the method. The following Equations (1)-(7) are from Reference [18] and are presented here for understanding the method.
The propagation of an ISW along the -direction is described by the KdV equation: where ( , ) is the maximum vertical displacement of isopycnal surfaces, is the spatial variable in the direction of wave propagation, is the time, is the linear wave speed of ISWs, and and are the quadratic non-linear and dispersion parameters, respectively. For a continuously stratified system, these parameters are obtained as follows: where is the water depth, ( ) is the vertical structure function for the first mode, which is normalized by its maximum absolute value. ( ) and can be obtained by solving the following Sturm-Liouville equation for a Boussinesq fluid in two dimensions and with no current shear or rotation: where ( ) is the buoyancy frequency, which can be obtained by in situ CTD measurements or climatological data. An analytical solution of Equation (1) is as follows: where is the characteristic half width, is the ISW amplitude, and is the non-linear phase speed. These parameters are obtained as follows:

The Conventional Method of ISW Amplitude Estimation Based on the Classic KdV Equation in a Continuously Stratified Ocean Model
The conventional ISW amplitude estimation method applying the classic KdV equation in a continuously stratified ocean was developed by Small et al. [18]. Here, we briefly describe the method. The following Equations (1)-(7) are from Reference [18] and are presented here for understanding the method.
The propagation of an ISW along the x-direction is described by the KdV equation: where η(x, t) is the maximum vertical displacement of isopycnal surfaces, x is the spatial variable in the direction of wave propagation, t is the time, c 0 is the linear wave speed of ISWs, and α and β are the quadratic non-linear and dispersion parameters, respectively. For a continuously stratified system, these parameters are obtained as follows: where H is the water depth, φ(z) is the vertical structure function for the first mode, which is normalized by its maximum absolute value. φ(z) and c 0 can be obtained by solving the following Sturm-Liouville equation for a Boussinesq fluid in two dimensions and with no current shear or rotation: where N(z) is the buoyancy frequency, which can be obtained by in situ CTD measurements or climatological data. An analytical solution of Equation (1) is as follows: where l is the characteristic half width, η 0 is the ISW amplitude, and c is the non-linear phase speed. These parameters are obtained as follows: Therefore, if we have the ISW characteristic half width l and the background stratification profile, the ISW amplitude can be calculated using Equation (6). The background stratification profile can be obtained from in situ measurements or climatological data. The characteristic half width of ISWs can be derived by the curve fitting method proposed in Reference [20]. However, since the expression for ISW signature on SAR images used in the curve fitting method of Reference [20] is based on a two-layer theory, we need to deduce the expression for ISW signature on SAR images in a continuously stratified ocean.
For a given radar wave number k 0 and incidence angle θ, the NRCS (denoted σ 0 ) depends on the surface capillary-gravity wave spectrum density [37]. Under the equilibrium and steady condition, the variations of surface capillary-gravity wave spectra density are mainly modulated by three terms: wind input, viscosity dissipation, and ISW-induced current [19,20]. On the horizontal scale of an internal wave packet, the first two terms of wind input and viscosity dissipation are uniform, and thereafter only contribute to the background of an internal wave SAR image. Thus, the modulation effects caused by the internal waves on σ 0 are what we focus on in this study. Following the ideas presented in Reference [19], the theoretical expression for the soliton-induced radar backscatter cross-section per unit area (σ ISW 0 ) is given in Equation (8).
where θ and k 0 represent the incidence angle and wave number of radar waves, respectively, k sw = 2k 0 sinθ indicates the wave numbers of the resonant surface Bragg waves, χ represents the SAR beam looking angle relative to the wave direction, ω is the angular frequency of the resonant surface Bragg waves, and m 3 is the non-dimensional constant. The indices i and j denote the polarizations of the incident and backscattered microwaves, respectively, and g ij (θ) represents the first-order scattering coefficients. u is the ISW-induced horizontal velocity, which can be obtained as follows: Substituting Equation (5) into Equation (9) and then substituting the result into Equation (8), we have the following equation: where M = 16 By taking a normal transect through an ISW pattern on the SAR image as the transect AA' shown in Figure 2a, we can obtain a horizontal distribution curve of pixel backscattering intensity value I. Fitting the curve to Equation (10), we know that I can be expressed in the following form: where A, B, and C are the fitting parameters. From Equation (11), we can see that the pixel values of ISWs on the SAR image are relevant to the four parameters (A, B, l, and C). These parameters are determined by the best fitting curve to the pixel values along a transect across an ISW on the SAR image by Equation (11). Having the ISW characteristic half width l and the background stratification profile, one can calculate the ISW amplitude using Equation (6).

The Proposed Method of ISW Amplitude Estimation Based on the eKdV Equation in a Two-Layer Ocean Model
As mentioned in the Introduction, a simplified two-layer oceanic model is more practical than a continuously stratified model for estimating the ISW amplitude. Thus, the proposed method was applied to a two-layer ocean and employs the eKdV equation because the eKdV equation is a good indicator of the ISW amplitude estimation in the two-layer shallow water [21].
The eKdV equation is written as follows [38,39]: where α 1 is the cubic non-linearity parameter. For a two-layer system with a rigid lid and no background flow, in the Boussinesq approximation, the following Equations (13)-(16) can be obtained: where h 1 is the upper layer thickness and h 2 is the lower layer thickness. For depression wave, h 1 < h 2 . ∆ρ/ρ 0 = 2(ρ 2 − ρ 1 )/(ρ 2 + ρ 1 ) is the relative layer density difference, where ρ 1 (ρ 2 ) is the uniform density of the upper (lower) layer, which is defined as follows [14]: here, ρ(z) is the background density profile. From Equation (14), we know that α has a negative sign if considering the case in which the wave amplitude η 0 < 0 (depression wave). According to Equation (15), the dispersion parameter β is positive. As the term inside the square bracket of Equation (16)

can be equally transformed into
, α 1 is negative. One of the ISW solutions to the eKdV equation is as follows: where here, b is a parameter that satisfies b < 1. When b > 1, solution (19) would contain meaningless points of discontinuity, which make the denominator of solution (19) equal to 0. As the upper limit ( b → 1 ) is approached, the waves begin to broaden until a limiting table-topped wave amplitude is reached.
To obtain the eKdV type σ ISW 0 , we substitute solution (19) into Equation (9) and then substitute the result into Equation (8). Therefore, σ ISW 0 is written as follows: dφ dz |z = 0 . Similar to Equation (11), Equation (23) is modified into the following form: where I still represents the SAR image pixel backscattering intensity value, and A, B, and C are, again, the fitting parameters of I to Equation (24). From Equation (24), the pixel backscattering intensity value I is found to be related to the five parameters A, B, b, γ, and C. Figure 6 shows an illustrated plot of I versus ξ for setting A = −1, B = 0, C = 0, γ = 1, and b = 0.5. In the I − ξ plot, there is one maximum point (ξ max , I max ) and one minimum point (ξ min, I min ), marked by red dots. According to Equation (24), we obtain that ξ max and ξ min are as follows: Therefore, ξ max and ξ min satisfy the following equation: I max and I min are obtained as follows: From Equations (28) and (29), we obtain the following equation: Based on Equations (27) and (30), parameters B and C can be directly obtained after locating the maximum and minimum points in the pixel backscattering intensity profile of an ISW on the SAR image. From Equations (20)- (22), (28), and (29), the other three parameters γ, b, and A are obtained: According to Equation (32), the amplitude η 0 is dependent on the upper layer thickness h 1 , total water depth H, background density profile ρ(z), and phase speed c. In this study, the water depth was determined from a local nautical chart, the background density profile was taken from in situ measurements or WOA13 dataset, the phase speed was determined by using two consecutive remote sensing images, and the upper layer thickness was found by a best fitting curve of the pixel values of an ISW using Equation (24). Finally, the amplitude was obtained in terms of Equation (32). Actually, the idea behind using consecutive remote sensing image pairs to find accurate upper layer thickness and estimating consequent ISW amplitude is similar to that of Reference [16] in which Doppler velocities from along-track InSAR data were considered to find a plausible upper layer thickness and internal soliton amplitude within their expected ranges. Moreover, a point that must be addressed is that another wave amplitude solution η 0 = −α + α 2 + 6α 1 (c − c 0 ) /α 1 can also be acquired by solving Equation (20), which is not considered here because in this case, the parameter b = α − α 2 + 6α 1 (c − c 0 ) / α + α 2 + 6α 1 (c − c 0 ) was greater than 1 for α 1 < 0 and α < 0. According to Equation (32), the amplitude is dependent on the upper layer thickness ℎ , total water depth , background density profile ( ), and phase speed . In this study, the water depth was determined from a local nautical chart, the background density profile was taken from in situ measurements or WOA13 dataset, the phase speed was determined by using two consecutive remote sensing images, and the upper layer thickness was found by a best fitting curve of the pixel values of an ISW using Equation (24). Finally, the amplitude was obtained in terms of Equation (32). Actually, the idea behind using consecutive remote sensing image pairs to find accurate upper layer thickness and estimating consequent ISW amplitude is similar to that of Reference [16] in which Doppler velocities from along-track InSAR data were considered to find a plausible upper layer thickness and internal soliton amplitude within their expected ranges. Moreover, a point that must be addressed is that another wave amplitude solution = (− + + 6 ( − ))/ can also be acquired by solving Equation (20), which is not considered here because in this case, the parameter = ( − + 6 ( − ))/( + + 6 ( − )) was greater than 1 for < 0 and < 0.

Results and Analysis
This section presents the results of the ISW amplitude estimation in the Hainan case obtained using the two methods described above.

Results and Analysis
This section presents the results of the ISW amplitude estimation in the Hainan case obtained using the two methods described above. Figure 7 shows the best fitting curve between the pixel backscattering intensity value of the soliton L1 observed on the GF-3 SAR image and the theoretical model result given in Equation (11). The pixel backscattering intensity values of the soliton L1 observed on the SAR images (the red points in Figure 7) are written as I 01 , I 02 , I 03 , I 04 , . . . . . . I 0n , respectively. The optimal fitting coefficients (A m , B m , l m , C m ) obtained by the curve fitting among the observation values (I 0i , i = 1, 2, 3, . . . . . . n) and the theoretical model (11) are −10.75, 1.35, 144.93 m, and −11.66, respectively. Here, the curve fitting method rather than the Cramér-Rao Bound method [31] was applied for estimating the characteristic half width. It is because the amplitude of the ISW backscattering intensity defined in Reference [31] was small (about 0.85) in our case, from which we could deduce that the curve fitting and Cramér-Rao Bound methods may have little difference in the characteristic half width estimation (refer to Figure 4b in Reference [31]). Thus, the curve fitting method was adopted here for its relative simplicity. Determined by the background density profile shown in Figure 4, the quadratic parameter α and dispersion parameter β were −0.0158 and 157.06, respectively. Therefore, the amplitude of the soliton L1 using this conventional method was estimated to be −5.66 m.

The Derived ISW Amplitude Using the Conventional Method
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 18 characteristic half width. It is because the amplitude of the ISW backscattering intensity defined in Reference [31] was small (about 0.85) in our case, from which we could deduce that the curve fitting and Cramér-Rao Bound methods may have little difference in the characteristic half width estimation (refer to Figure 4b in Reference [31]). Thus, the curve fitting method was adopted here for its relative simplicity. Determined by the background density profile shown in Figure 4, the quadratic parameter and dispersion parameter were −0.0158 and 157.06, respectively. Therefore, the amplitude of the soliton 1 using this conventional method was estimated to be −5.66 m. Small et al. [18] concluded that the conventional method of using the ISW characteristic half width estimated from SAR images to derive soliton amplitudes can provide basic estimates of amplitude, but the results contain some uncertainties and may be affected by some factors, such as the characteristic half width, water depth, and background stratification. Therefore, this estimation of −5.66 m cannot be directly treated as a real value for comparison with the amplitude derived using the proposed method. To clarify the reliability of this amplitude estimation, we conducted the following uncertainty analysis. Here, the analysis method was similar to that adopted in Reference [19].
The theoretical model (11) in the case of the optimal curve fitting estimation can be rewritten as follows: The total root mean square (RMS) deviation of the pixel backscattering intensity values between the observations on SAR images and theoretical model (34) is as follows: Small et al. [18] concluded that the conventional method of using the ISW characteristic half width estimated from SAR images to derive soliton amplitudes can provide basic estimates of amplitude, but the results contain some uncertainties and may be affected by some factors, such as the characteristic half width, water depth, and background stratification. Therefore, this estimation of −5.66 m cannot be directly treated as a real value for comparison with the amplitude derived using the proposed method. To clarify the reliability of this amplitude estimation, we conducted the following uncertainty analysis. Here, the analysis method was similar to that adopted in Reference [19].
The theoretical model (11) in the case of the optimal curve fitting estimation can be rewritten as follows: The total root mean square (RMS) deviation Dev of the pixel backscattering intensity values between the observations on SAR images and theoretical model (34) is as follows: From model (34), we know that the Dev of pixel backscattering intensity values is caused by the uncertainties in the four parameters A m , B m , l m , and C m . Thus, the Dev can also be written as follows: where ∆l m is the uncertainty in l m and As the amplitude η 0 is a function of α, β, and l (refer to Equation (6)), and α and β, which are derived from the field observation profile, are close to the real value, the uncertainty in the amplitude estimation (∆η 0m ) is caused by only ∆l m and can be written as follows: According to Equation (38), ∆η 0m is directly proportional to ∆l m . To give a maximum estimation of the uncertainty in the amplitude retrieval, we assume that Dev = ∂I m ∂l |∆l m | by neglecting the contributions of the three parameters A m , B m , and C m to Dev in Equation (36). Therefore, the uncertainty in parameter l m can be derived by the following equation: Substituting the optimal curve fitting parameters A m , B m , l m , and C m into Equations (35) and (37)-(39), respectively, we obtain Dev, ∆l m , and ∆η 0m as 0.32, 15.87 m, and 1.24 m, respectively. Thus, the amplitude of the soliton L1 should be −5.66 m with an uncertainty of ± 1.24 m.

The Derived ISW Amplitude Using the Proposed Method
From the profile of the pixel backscattering intensity value of the soliton L1 (Figure 7), we can derive the maximum point (ξ max , I max ) and minimum point (ξ min, I min ) located at approximately (1.26, −7.50) and (1.44, −15.62), respectively. Thus, B and C are estimated to be 1.35 and −11.56, respectively. As the internal soliton usually features a wide stripe on SAR images, we define the point where the positive peak of the ISW signature is on the SAR image as its location. Therefore, the propagation distance of an ISW within a short time is defined as the distance between two positive peaks, whose geographic locations can be determined from SAR images with the help of Sentinel Toolbox software (http://step.esa.int/main/download/). In the present study, the propagation distance of the soliton L1 determined from the two consecutive GF-3 and TS-X SAR images within 11 min was 435.6 m (thus, the phase speed was approximately 0.66 m/s). From the background density profiles shown in Figure 4, it is difficult to exactly determine the upper layer thickness, which seems to be in a large range of 9 m to 33 m. Using a simple iterative approximation scheme, the optimal upper layer thickness of 23 m was eventually determined by making the best fitting curve between Equation (24) and the SAR image pixel values along the transect AA' of the soliton L1. Figure 8 shows the best fitting curve between the pixel backscattering intensity value of the soliton L1 observed on the GF-3 SAR image and the theoretical model given in Equation (24). Therefore, the amplitude of soliton L1 was subsequently calculated to be approximately −4.52 m, which agrees with the amplitude estimation of −5.66 ± 1.24 m derived by the conventional KdV equation in a continuously stratified ocean.
A great challenge of using SAR images to retrieve the ISW amplitude is the need for synchronous in situ background stratification data. Here, we showed that our method could partially overcome this challenge using the climatological data WOA13 V2. The density profile ( Figure 5) obtained from the WOA13 V2 climatological datasets was used to derive the amplitude of soliton L1 in the Hainan case. Our method yields an amplitude of −7.35 m, which is comparable to the reference value of −5.66 ± 1.24 m. However, if the density profile of WOA13 V2 was used for the conventional method based on the KdV equation in a continuously stratified ocean, the estimated amplitude would be −16.30 m, which is almost three times the reference value. Clearly, the proposed method yields more credible results than the conventional method in the case where in situ measurements are unavailable. This result may be because the 5 m depth vertical interval of the background stratification profile from the WOA13 V2 dataset is relatively coarse for the amplitude estimation in shallow water using the conventional method. In contrast, the coarse oceanic density profile obtained from the WOA13 V2 dataset has fewer effects on the amplitude estimation using the proposed method, which suggests that the proposed method using consecutive SAR images and the eKdV equation in a two-layer ocean is a reliable way to estimate the ISW amplitude independent of in situ oceanic density profile measurements.
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 18 WOA13 V2 dataset is relatively coarse for the amplitude estimation in shallow water using the conventional method. In contrast, the coarse oceanic density profile obtained from the WOA13 V2 dataset has fewer effects on the amplitude estimation using the proposed method, which suggests that the proposed method using consecutive SAR images and the eKdV equation in a two-layer ocean is a reliable way to estimate the ISW amplitude independent of in situ oceanic density profile measurements.

Discussion
Determining the upper layer thickness is a crucial step in estimating ISW amplitude from SAR images. In some reported cases, the sea water features a thin and strongly stratified pycnocline, sandwiched between two homogeneous layers, and therefore the upper layer thickness can be accurately estimated. Two common approaches have been proposed to determine the upper layer thickness from the background density profile. One approach is the eigenfunction method, which defines that the depth where the first mode eigenfunction has its maximum is the upper layer thickness. The other one is the method, which determines the depth of the maximum buoyancy frequency as the upper layer thickness. If the eigenfunction method is applied to the Hainan case, based on the two-layer eKdV theory, the upper layer thickness and amplitude of the soliton 1 are estimated to be approximately 29 m and −5.37 m, respectively. If the method is applied, the upper layer thickness is approximately 13 m, and the amplitude of the soliton 1 is approximately −12.01 m. Therefore, one can find that the uncertainties of upper layer thickness estimation can lead to large discrepancy of amplitude estimation. The stratification profile in the Hainan case suggests that the uncertainties of upper layer thickness estimation arises from the typical stratification of a thick pycnocline, which makes the determination of the upper layer thickness in a two-layer model rather ambiguous and difficult. Our proposed method solves this problem by using tandem SAR images to determine the accurate ISW phase speed, which then is used to estimate the optimum upper layer thickness and the consequent ISW amplitude.

Discussion
Determining the upper layer thickness is a crucial step in estimating ISW amplitude from SAR images. In some reported cases, the sea water features a thin and strongly stratified pycnocline, sandwiched between two homogeneous layers, and therefore the upper layer thickness can be accurately estimated. Two common approaches have been proposed to determine the upper layer thickness from the background density profile. One approach is the eigenfunction method, which defines that the depth where the first mode eigenfunction has its maximum is the upper layer thickness. The other one is the N max method, which determines the depth of the maximum buoyancy frequency as the upper layer thickness. If the eigenfunction method is applied to the Hainan case, based on the two-layer eKdV theory, the upper layer thickness and amplitude of the soliton L1 are estimated to be approximately 29 m and −5.37 m, respectively. If the N max method is applied, the upper layer thickness is approximately 13 m, and the amplitude of the soliton L1 is approximately −12.01 m. Therefore, one can find that the uncertainties of upper layer thickness estimation can lead to large discrepancy of amplitude estimation. The stratification profile in the Hainan case suggests that the uncertainties of upper layer thickness estimation arises from the typical stratification of a thick pycnocline, which makes the determination of the upper layer thickness in a two-layer model rather ambiguous and difficult. Our proposed method solves this problem by using tandem SAR images to determine the accurate ISW phase speed, which then is used to estimate the optimum upper layer thickness and the consequent ISW amplitude.
While a case study shows the effectiveness of the proposed method, there are some aspects needed to be discussed.
The key to the proposed method is to seek two consecutive images to estimate the accurate phase speed. As the derived ISW phase speed in the proposed method is used to fit the theoretical model of radar backscatter intensity to the SAR observation for determining the optimum upper layer thickness, the image pair applied in the proposed method should include at least one SAR image. In this study, we use two consecutive SAR images (tandem SAR images) to drive the method. A development trend of modern spaceborne SAR is constellation configurations, such as Cosmo-SkyMed, TerraSAR-X and its twin Tandem-X, and the recently launched Radarsat constellation mission; therefore, acquisitions of tandem SAR images become increasingly feasible. In addition to the combination of tandem SAR images, a combination of SAR image and an optical image is also feasible to construct consecutive image pairs to apply the proposed method, which shall be tested in our further work. However, for constructing a consecutive image pair used for the proposed method, two aspects remain further studies. One is the geolocation accuracy of the two consecutive images. The geolocation accuracy directly determines the retrieval accuracy of the phase speed and the consequent amplitude. The other aspect is the temporal interval between the two consecutive images. The appropriate temporal interval for estimation of phase speed of ISWs should depend on the propagation and evolution characteristics of ISWs. In principle, the larger the variation of the ISW phase speed between two observations is, the shorter the temporal interval between two remote sensing images should be.
There is not yet a method for calculating the exact equivalent two-layer ocean model of a continuously stratified ocean [33]. A reasonable approach to judge whether a two-layer fluid model is or not a good representation of the actual continuous density profile is to compare the coefficients of the KdV-type equation (including the linear wave speed c 0 , quadratic non-linear parameter α, cubic parameter α 1 and dispersion parameter β). In this study, the derived four coefficients of c 0 , α, α 1 , and β in the two-layer eKdV equation using the derived optimum upper layer thickness of 23 m are 0.63, −0.0226, −0.0018, and 123.54, respectively, which agree with the four coefficients estimated based on the continuous stratified eKdV equation [40] without considering current shear or rotation with values of 0.60, −0.0158, −0.0003, and 157.06, respectively. This should be the reason why the ISW amplitude derived by the conventional method and the proposed method is close to each other in the Hainan case.

Summary and Conclusions
Amplitude is one of the most important parameters of ISWs, while the retrieval of their amplitudes from satellite remote sensing images, which generally reflects sea surface variations, remains difficult, although observations of ISWs from space have been available for a few decades. Combining satellite remote sensing images, e.g., SAR images, and the KdV equation for application in a continuous stratification can generally yield reasonable ISW amplitude estimations. However, this method requires simultaneous in situ measurements of background stratification (density) profiles, which are often not available. To overcome this difficulty, the simplified two-layer ocean model was introduced, which improves the practicability of deriving ISW amplitudes from SAR images. However, uncertainties in estimating upper layer thickness in the two-layer ocean may lead to large discrepancies in amplitude estimations. The Hainan case shown in this study is a typical one with such an existing such problem, where the thick pycnocline relative to water depth makes it difficult to determine the upper layer thickness from the oceanic density profile. Therefore, we proposed a new method to solve this problem in the retrieval of ISW amplitudes from SAR images.
In the proposed method, we firstly developed a theoretical model describing the radar backscatter cross-section of ISWs based on the eKdV theory in a two-layer ocean. By fitting the theoretical model estimation of radar backscatter to the SAR observations, we exploited consecutive SAR images to determine the optimum upper layer thickness for amplitude estimation. The retrieved ISW amplitude in the Hainan case agrees well with the result estimated by the classic KdV equation in the measured stratification. A key parameter needed in this method is the ISW phase speed. Different from the previous methods (assuming phase speed is equivalent to the group velocity), our estimation of phase speed was based on two consecutive SAR images acquired with a short temporal interval on a scale of dozens of minutes, which ensures derivation of an instantaneous phase speed and further leads to accurate estimations of ISW amplitudes. Besides, using the representative climatological dataset of WOA13 V2 instead of in situ measurements can also yield a more reasonable estimation of ISW amplitudes than the results using the conventional method based on the KdV theory in a continuously stratified ocean in this shallow water case. This result suggests that the proposed method can provide basic yet reliable information on ISW amplitudes even when in situ measurements are absent.
In conclusion, the present study demonstrates a method of retrieving ISW amplitude from SAR images, which are characterized by a consecutive image pair within a short time interval. Since the method has a potential capability to get reliable amplitude estimations when the in situ measurements are absent, it could have a wide and attractive application in retrieving ISW amplitude under the trend that modern spaceborne SAR tends to be in constellation configurations.