Estimation of Propagation Speed and Direction of Nonlinear Internal Waves from Underway and Moored Measurements

: Propagation speed and direction of nonlinear internal waves (NLIWs) are important pa ‐ rameters for understanding the generation and propagation of waves, and ultimately clarifying re ‐ gional ocean circulation. However, these parameters cannot be directly measured from in ‐ situ in ‐ struments , but can only be estimated from post ‐ processing in situ data. Herein, we present two methods and an optimal approach to estimate the propagation speed and direction of waves using underway and moored observations. The Doppler shift method estimates these parameters from ap ‐ parent observations concerning a moving ship using the Doppler shift induced by the changing relative distance of the NLIWs from the moving ship. The time lag method estimates the parameters using the distance between two locations of the NLIW observed at different times and the time lag. To optimize the speed and direction of NLIWs, the difference in the propagation direction inde ‐ pendently estimated by the two methods needs to be minimized concerning the optimal propaga ‐ tion speed to yield the optimal propagation direction. The methods were applied to two cases ob ‐ served in the northern East China Sea in May 2015 and August 2018. This study has practical sig ‐ nificance for better estimating the propagation speed and direction of NILWs particularly over a broad continental shelf.


Introduction
Nonlinear internal waves (NLIWs) are ubiquitous in stratified seas and are accompanied by isopycnal fluctuations with a sharp vertical density gradient. They play an important role in underwater acoustics, regional circulation, local biogeochemistry, and energetics, mostly via vertical mixing in the stages of generation, propagation, evolution, and dissipation. NLIWs affect the transportation of momentum, heat, and energy via turbulent dissipation and mixing [1][2][3][4][5]. Marine ecosystems are significantly influenced by vertical nutrient supply, chlorophyll bloom, and biological redistribution, which can be modulated by NLIWs [6][7][8][9]. The NLIWs drive sediment resuspension and transportation; thus, they affect marine geophysics and underwater acoustics [10][11][12][13][14][15]. Vertical isopycnal displacements, which allow the wave amplitude to be defined, and propagation speed and direction, are fundamental parameters of NLIWs that are useful, but cannot be directly measured from in situ sampling, for a clear understanding of their generation, propagation, evolution, and dissipation. Estimating the propagation speed and direction can be important for assessing regional ocean circulation, biogeochemical cycles, energetics, underwater acoustics, and the dynamics of NLIWs.
Methods to estimate the propagation speed and direction have been suggested but are mostly limited by sampling strategies that have not yet been validated. The most common method using multiple moorings aligned in the propagation direction of NLIWs aims to divide the distance between the mooring locations by the arrival time differences [16][17][18][19][20][21][22]. However, it is not practical to deploy many moorings along the ray of NLIWs, particularly where the continental shelf is wide and multiple NLIWs are generated from multiple sources with different unknown propagation directions. Another method is to use the principal direction of the wave-induced horizontal velocity and the temporal difference of enhanced echo intensity from acoustic Doppler current profiler (ADCP) measurements [23][24][25][26][27]. This method is useful, but not very practical, as extracting the propagation speed and direction is not straightforward. Using remote sensors, such as synthetic aperture radar (SAR) and spectroradiometer, the propagation speed and direction can be estimated from the horizontal curvature of satellite images [28][29][30][31][32]; however, the limited spatiotemporal satellite sampling from polar orbits does not allow NLIWs to be easily detected. Therefore, it is necessary to develop a method to estimate the propagation speed and direction of NLIWs from widely used ship-based in situ measurements.
In the northern East China Sea (ECS), NLIWs are mainly formed by strong tidal forces that interact with bathymetric features. NLIWs in this region have been observed in association with strong semidiurnal internal tides over slope areas in the southern and southeastern parts of the ECS [33] and local lee-wave generation by small islands and seamounts near Jeju Island and the Ieodo Ocean Research Station (IORS) in the northern ECS [34][35][36] (Figure 1). Unlike the typical setting where dominant first-mode NLIWs in a twolayered condition propagate from the shelf break towards the coast, high NLIW modes propagating in multiple directions from multiple sources have been identified in the northern ECS [34,37]. Herein, we present a new method for estimating the propagation speed and direction of NLIWs using both moored and underway measurements, and the results of applying the method to two cases of NLIWs observed in the northern ECS in May 2015 and August 2018 ( Figure 1).

Data and Processing
Shallow-water Acoustic Variability EXperiment 2015 (SAVEX15) was conducted on 14-28 May 2015, focusing on a relatively small area (water depth: ~100 m, area shown in Figure 1a,b) in the northern ECS [38][39][40][41][42]. During the experiment, two moorings (vertical line arrays [VLAs]) and underway conductivity-temperature-depth (UCTD) instruments were used to collect moored temperature (with no conductivity) time series at multiple depths and vertical profiles [43]. The two moorings (VLA1 and VLA2) deployed at water depths of 101 and 102 m, were horizontally separated by ~5.5 km; 25 temperature loggers and 5 Star-Oddi temperature-depth-tilt sensors were attached at nominal depths of 2-80 m with an interval of 1-5 m. The sampling time interval of the moored temperature sensors was 30 s.
Ieodo Ocean Research Station 2018 (IORS18) was conducted in the northern ECS in the vicinity of IORS (32°7.4 N, 125°10.9 E, constructed at a water depth of 41 m) on 28 August-1 September 2018 (Figure 1a,c). During the IORS18 experiment, UCTD was used to collect ship-based vertical profiles of temperature and salinity and the IORS-based time series of temperature and pressure data observed at nominal depths of 2, 5, 11, 16, 22, 32, and 37 m, with a typical sampling time interval of 60 s [44].
The depths of the moored temperature sensors attached to the VLAs were corrected using the tilt and pressure data recorded by five Star-Oddi temperature-depth-tilt sensors. After the removal of outliers, the moored temperature data were vertically interpolated using the Akima spline method [45] at 1 m intervals. The UCTD data were processed following the method described by Ullman and Herbert [46], except for the alignment process to correct the mismatch due to different time delays of the conductivity and temperature sensors. The raw temperature and conductivity measured by the UCTD were filtered with a cut-off period of four scans (0.25 s). The raw pressure measured using the UCTD was filtered with a cut-off period of 32 scans (2 s). The time delay between the conductivity and temperature sensors was corrected using lagged correlation. Then, spikes in the salinity data were removed by aligning the data of the temperature and conductivity sensors. After the alignment processing, the data were corrected for the effect of viscous heating and finally vertically averaged over 1 bin. To discuss the theoretical propagation speed of NLIWs in the northern ECS in the context of long-term and interannual variations, vertical profiles of temperature and salinity routinely observed every other month from 1990 to 2019, at two hydrographic stations of the National Institute of Fisheries Science (NIFS), Republic of Korea, were used in this study (red open circles in Figure 1a). To ensure the quality of the temperature and salinity data, vertical profiles containing unreasonable values (both global and local) were removed. Quality control procedures, such as the spike and gradient tests, were applied to extract reliable salinity and temperature profiles. Because the profiles are only available at standard depths (i.e., surface, 10,20,30,50,75,100,125,150,200, and 250 m), linear interpolation was conducted to determine data at 1 m vertical intervals [47].
Moderate resolution imaging spectroradiometer (MODIS) sensors onboard the National Aeronautics and Space Administration (NASA) satellites Terra and Aqua, provided true-color images from calibrated, corrected, and geo-located radiance (Level-1 B) data, with a spatial resolution of 250 m. As NLIWs induce the divergence and convergence of sea surface currents as they propagate, thereby modifying the sea surface roughness, they are visible in MODIS true-color images if they are in a sun-glint area [48]. In this study, two images obtained by MODIS Terra on 2 August 2015, and MODIS Aqua on 30 July 2018, were used to estimate the propagation direction of NLIWs from sea surface manifestations.

Two-Layered KdV (Korteweg-de Vries) Theories
In classical KdV theory [49], a leading-order weak non-linearity and dispersion are competing but comparable to each other. For the two-layered KdV theory, the thicknesses (ℎ , ℎ ) and densities ( , ) of the upper and lower layers can be used to estimate the parameters of mode-1 NLIWs, including linear phase speed . , theoretical propagation speed . , characteristic width 2∆ . , nonlinear parameter , and dispersion parameter , yielding the wave equation as follows [14,50]: where , , and are the vertical displacement of the isopycnals (or isotherms), time, and horizontal coordinates, respectively. The . , , and can be estimated using the density stratification parameters ( , , ℎ , and ℎ ) in a two-layered system as follows: . , . ℎ ℎ 6 Here, is the gravity acceleration set to 9.80 m s . The thicknesses of the upper and lower layers (ℎ and ℎ ) were determined based on the depth of the maximum density gradient from the density profiles obtained from the UCTD. The densities at the upper and lower layers ( and , respectively) were determined as the minimum density within the upper layer and the maximum density within the lower layer, respectively. The solution of Equation (1) for the displacement , yields the nonlinear soliton as follows: Here, the theoretical propagation speed . and characteristic width 2∆ . were calculated from the linear phase speed . and the amplitude ( ) of the vertical displacement of are as follows: . .
is a cubic nonlinear parameter in the two-layer system. The theoretical propagation speed . and characteristic width 2∆ . based on the eKdV theory in the two-layered system are as follows:

Doppler Shift Method
To estimate the propagation direction of NLIWs using the Doppler shift caused by propagating NLIWs observed from a moving ship, the theoretical propagation speed . and ship speed were assumed to be constant during the measurement period, and the propagation direction was assumed to be orthogonal to the constant phase lines parallel to the wavefront lines (Figure 2a,b). As the estimated propagation direction is in reference to the ship course , the apparent propagation speed can be represented as the difference between . and the ship speed in direction as follows: . cos .
Here, | | is the angular difference between the ship course and the propagation direction of the NLIWs. Because the Doppler-shifted apparent frequency or the inverse of the apparent period can be represented by and is the wavelength of the NLIWs and the Doppler equation [52], the following equation can be used to estimate the : Here, is determined from measurements (Table 1), while . is determined by the Cnoidal model [50] as where is a complete elliptic integral of the first kind and parameter is set to 0.5. Equation (11) can then be rewritten using Equation (12) as follows: Further, the is obtained as follows: The propagation direction of the NLIWs estimated using the method described above has an angular ambiguity caused by the sign of the arccosine part. Thus, a physically reasonable direction between the two was selected. The and are angles in degrees measured counter-clockwise from the east (for example, 180° and 270° correspond to the westward and southward directions, respectively).

Time Lag Method
Independent of the method described in Section 2.2.2, the propagation direction of NLIWs was estimated using the distance between two locations of the NLIW front observed at different times, with the assumption that the NLIWs propagate across the two measurement locations with an angle orthogonal to the constant phase lines at a constant speed (Figure 2a,c). The observed propagation speed was estimated by dividing the distance between the two locations by the arrival time lag , for example, . Then, . was calculated from and angular difference | | between (direction from the first measurement location to the second measurement location) and the propagation direction of NLIWs as Finally, the was obtained from Equation (16) as The propagation direction estimated using the method described above (time lag method) also has an angular ambiguity caused by the sign of the arccosine part. Thus, a physically reasonable direction is selected. The and are angles in degrees measured counter-clockwise from the east.

Estimation of Propagation Speed and Direction
Two propagation directions of NLIWs estimated from the two methods were used to estimate the optimal propagation direction and successive propagation speed. First, a consistent direction between the two directions derived from each method was selected to minimize the ambiguity where the two methods yield angular difference | | typically less than 30°. For example, each one (bold green and blue colors in Figure 2b,c) among two and two are selected as more consistent between the two methods and physically reasonable, while inconsistent and (deemed green and blue colors in Figure 2b,c) among the two and two were not selected. To optimize the propagation speed and direction, the difference in consistent propagation directions from the two methods | | was minimized by iteratively changing . at intervals of 0.01 m • s instead of using the constant propagation speed derived from the twolayered KdV theory. The updated propagation speed was determined from the iterations, and the final propagation direction of NLIWs was determined by averaging the two directions /2 when | | reached its minimum value for the updated propagation speed. Iterations were performed for a range of 30% deviation from .
(typically requiring 38 iterations to reach the minimum), which is comparable to the range of the interannual variation of propagation speed reported in a previous study [47], as discussed in Section 4.

Estimation of Propagation Direction Using Satellite Images
The propagation directions from the MODIS images ( Figure 3) were estimated from the horizontal curvature of the leading fronts from the sea surface manifestations [28]. The propagation direction was calculated by the direction of the center of the straight line connecting the endpoints of the manifestations to the center of the arc for the manifestations. For example, a manifestation of NLIW is shown in the blue box in Figure 3a,c. The orange curvature line from both endpoints (Points A and B) is the leading front of the NLIW. Point C is the center of line AB. Point D is the center of the arc AB. In this case, the angle between the vectors from C to D, measured counter-clockwise from the east, is the propagation direction.

SAVEX15
On 23 May 2015 during the SAVEX15, the existence of NLIWs was confirmed from a series of UCTD profiling measurements, for example, at 15:56 UTC (vertical red dashed line in Figure 4b), and the time-depth pattern of temperature variations observed at the northern mooring VLA1 at 15:02 UTC (vertical grey dashed line in Figure 4a). The NLIW was observed at VLA1 3239 s earlier than that observed in the UCTD measurements (Figures 4 and 5a). The 13 °C isotherm displacements, derived from the UCTD and VLA1 observations, commonly range from 27 to 35 m (Figure 5a). The NLIW has an amplitude ( ) of 6.1 and 5.9 m estimated from UCTD and VLA1 observations, respectively ( Figures  4 and 5a). The two-layered theoretical parameters of the NLIWs are listed in Table 1 Table 1).   When the ship moved at a speed of 0.47 m • s and a direction of 148° (northwestward), the apparent propagation direction of NILWs had an angular difference of 60° with ship course (Figure 6a), derived from the Doppler shift method using Equation (15), resulting in 208° (southwestward) or 88° (northward). From the distance ( = 2233 m) between the two measurement locations (VLA1 and ship) and the time lag of the NLIW arrivals ( = 3239 s), the observed propagation direction of NILWs was estimated to have an angular difference of 19° with (Figure 6b), derived from the time lag method using Equation (17), resulting in 211° (southwestward) or 249° (slightly more southwestward). Thus, more consistent propagation directions of 208° and 211° were selected to optimize the propagation speed and direction. By minimizing | |, the optimal propagation speed ( ) of 0.64 m • s was derived from the iterative calculations, yielding 208° and 211° with | | 3°, and the resultant propagation direction ( ) of 210° (southwestward).

IORS18
On 30 August 2018 during the IORS18, the existence of NLIWs was confirmed from the isotherm displacements observed by a series of vertical profiling UCTD measurements, particularly at 11:11 UTC (vertical red dashed line in Figure 7b), and the time-depth pattern of temperature measurements at the IORS, at 10:32 UTC (vertical grey dashed line in Figure 7a). The NLIWs observed at IORS were 2370 s earlier than those observed at the UCTD (Figure 8a). The 18 °C isotherm displacements observed from the UCTD and IORS were comparable, ranging from 23 to 31 m (Figure 8a). The NLIWs (the leading NLIW among a set observed by the UCTD) had an amplitude ( ) of 6.8 and 7.0 derived from the IORS and UCTD observations, respectively (Figures 7 and 8a). The densities at the upper and lower layers ( and ) were 1018.18 and 1024.99 kg • m , and thicknesses of the upper and lower layers (ℎ and ℎ ) were 24.0 and 28.0 m, respectively ( Figure 8b and Table 1).  When the ship moved at a speed of 1.38 m • s and a direction of 334° (southeastward), the apparent propagation direction of NLIWs had the angular difference of 126° with the ship course (Figure 9a), derived from the Doppler shift method using Equation (15), resulting in 208° (southwestward) or 100° (northward). From the distance ( = 3045 m) between the two measurement locations (IORS and ship) and the time lag of the NLIW arrivals ( = 2370 s), the observed propagation direction of NLIW was estimated to have an angular difference 45° with (Figure 9b), derived from the time lag method using Equation (17), resulting in 198° (southwestward) or 288° (southeastward). Thus, more consistent propagation directions of 208° and 198° were selected to optimize the propagation speed and direction. By minimizing | |, optimal propagation speed ( ) of 1.06 m • s was derived from the iterative calculations, yielding 205° and 205° with | | 0°, and the resulting propagation direction ( ) of 205° was obtained.

Discussion
Herein, we discussed whether the propagation speeds of NLIWs estimated using the proposed method are reasonable based on the KdV theory and previous observations. Interannual variations of the theoretical propagation speed ( . , NIFS-SAVEX15) in May from 1994 to 2019 derived from the NIFS historical hydrographic data near the SAVEX15 area range from 0.36 to 0.71 m • s , with a temporal mean and standard deviation of 0.50 and 0.09 m • s , respectively (red line in Figure 10a). A long-term decreasing trend was observed in May  Figure 10a).
Similarly, interannual variations in the theoretical propagation speed ( . , NIFS-IORS18) in August from 1994 to 2019 derived from the NIFS historical hydrographic data near the IORS ranged from 0.42 to 0.86 m • s , with a temporal mean and standard deviation of 0.66 and 0.10 m • s , respectively (blue line in Figure 10a). However, the longterm trend in August .
(NIFS-IORS18) was positive (at a rate of 0.003 m • s • yr ; blue dotted line in Figure 10a) because of the increasing density stratification between the layers, that is, decreasing and increasing with no significant change in ℎ and ℎ in August (blue lines in Figure 10b  To determine whether the propagation directions of NLIWs estimated using the proposed method are physically reasonable, we compared the values with those derived from satellite images and previous observations. The values estimated during the two experiments (SAVEX15 and IORS18) using the proposed method were consistent with those derived from MODIS images. Despite the fact that the surface manifestations of NLIWs observed in the two MODIS images were distant (56-123 km from the SAVEX15 area and 42-190 km from the IORS) from the locations of NLIWs observed during the two experiments, south-westward-propagating NLIWs (propagation direction of 186-209°) were consistently found in the two images ( Figure 3). Previous observations based on satellite SAR and optical images taken between 1993 and 2004 [35] and SAR images taken between 2014 and 2015 [36] in the northern ECS also support the southwestward-propagating NILWs (propagation direction ranging from 212° to 245°), which may be dominant among the various NLIWs propagating in multiple directions from multiple sources in the northern ECS (Figure 1a).
The two-layered classical (ordinary) KdV theory used in this study has clear limitations. The classical KdV theory is very simplified and assumes weak non-linearity and weak dispersiveness. In fact, NLIWs observed in many areas have been better explained by the eKdV theory than by the KdV theory. However, propagation speeds, in cases of SAVEX15 and IORS18, derived based on the eKdV theory ( . ) including the cubic non-linearity are not significantly different from those based on the KdV theory ( . ), yielding the difference less than 0.02 m • s due to relatively small (Table 1). Furthermore, as in the case of NLIWs in the South China Sea, finite-depth theory may be theoretically more appropriate than shallow-water theory where the KdV and eKdV theories are based on [53]. The theoretical propagation speeds in the forms of Equations (6) and (9) are limited to the case of no background pedestal condition that could not be considered in this study and might affect the speed significantly. The rigid lid assumption of the KdV and eKdV theories at the top boundary is not fully realistic, although reasonable in many cases, because the resonant interaction between the surface and internal waves supports the possible need for the presence of a free surface at the top boundary to yield more realistic theoretical estimates [54]. In addition, the results presented in this study are limited to only mode-1 NLIWs by applying the two-layered system, yet the vertical profiles of mean density observed during the two experiments (Figures 5b and 8b) support normal mode decompositions (J. Klinck's Matlab program dynmodes.m, available online at http://github.com/sea-mat/dynmodes; accessed on 5 October 2021) for the first three modes corresponding to 49%, 18%, and 14% for SAVEX15 and 50%, 21%, and 13% for IORS18, respectively. Multi-mode NLIWs beyond the mode-1 NLIWs in the region, not investigated in this study, yet explaining about half of NLIWs, need to be examined in the future.

Concluding Remarks
We present a novel method to estimate the propagation speed and direction of NLIWs using widely collected underway and moored observations, and the results of applying the method to two cases of NLIWs observed in May 2015 (SAVEX15) and August 2018 (IORS18). Two-layered KdV theory and satellite images were used to discuss the results of the proposed method. The propagation direction of NLIWs was estimated with respect to a moving ship using the Doppler shift relationship (1. Doppler shift method) and independently using the time lag between the NLIWs observed at two different locations (2. time lag method). Then, the propagation speed and direction were optimized to minimize the difference in propagation directions derived from the two methods by iterating the propagation speed in the range of ±30% at a resolution of 0.01 m • s . The results derived from the proposed method are robust, as the range of iterative propagation speeds is comparable to the interannual variation of theoretical propagation speeds estimated using historical hydrographic data, yielding an error of less than 15% for the propagation direction. Because in situ observations of NLIWs are still challenging to collect and propagation speed and direction cannot be directly measured from subsurface instruments, our proposed method for estimating the propagation speed and direction of NLIWs using common underway and moored measurements is of practical importance, particularly over a broad shelf, such as the northern ECS, where the multi-directional propagation of multi-mode NLIWs from multiple sources is often observed.