## Abstract

Internal waves can transfer energy from large-scale to microscale processes; however, the spectra of these waves remain poorly known. A method that combines modal harmonic decomposition and maximum-likelihood method is proposed in this study to estimate four-dimensional internal wave spectrum using limited mooring observations. Using this method, a four-dimensional internal wave spectrum was obtained for the first time based on the mooring measurements collected during the South China Sea (SCS) Internal Wave Experiment in July 2014. The spectrum was then validated by comparing with the spectrum based on Fourier analysis and with the modified Garrett–Munk internal wave spectrum, respectively. The power of the internal wave spectrum decreased obviously with increasing frequency and wavenumber, with a falloff rate of *ω*^{−2} beyond tidal frequencies, and with falloff rates of $kh\u22122$ and $kz\u22122.5$ for horizontal and vertical wavenumbers, respectively. In addition, at a fixed frequency and vertical wavenumber, the propagation direction and phase speed of internal waves can be obtained through the four-dimensional spectrum. In summary, we verified the feasibility of estimating four-dimensional internal wave spectrum using limited mooring observations in this study, and the method we proposed should be applicable to other regions where such mooring observations are available.

## 1. Introduction

Internal wave spectrum reflects energy distribution of internal wave field as a function of wave frequency and wavenumber. Among processes with different scales, the wavenumber range of internal waves lies between large scale and turbulent scale, which serves as a continuum for energy transfer and eventually furnishes turbulent mixing in the ocean interior (Henyey et al. 1986). Therefore, the parameterization of small-scale mixing is highly dependent on the performance of internal wave spectrum in different regions of the ocean (Polzin et al. 1995, 2014).

There are several “universal” internal wave spectral models. The first one is the Garrett–Munk (GM) spectrum based on the assumption that internal wave field is horizontally isotropic and vertically symmetric. The GM model has evolved over time, and now it has several incarnations. GM72 spectrum was initially introduced by Garrett and Munk (1972) and is considered as a milestone on studies of internal wave spectrum. However, it was built on the wintertime measurements at site “D” on the continental rise north of the Gulf Stream. GM72 spectrum cuts off sharply beyond some specified wavenumber as a top-hat model, which ignores much of small scales. GM75 spectrum (Garrett and Munk 1975) replaces the top hat used in GM72 spectrum with a tapered cutoff to conform to the even higher wavenumber data from Millard (1972). GM75 spectrum is much closer to the real ocean spectrum than GM72 spectrum, as suggested by Garrett and Munk (1975). Cairns and Williams (1976) provided a modification of GM75 spectrum, which is referred to as GM76 spectrum; it has a more rapid roll-off from the low-wavenumber plateau region to the high-wavenumber asymptote than GM75 spectrum. For simplicity, GM79 spectrum (Munk 1981) was introduced by replacing the wavenumber used in GM76 with vertical mode number. These GM spectra have certain universality for internal wave field in the world oceans, particularly for deep open oceans, and one of their important roles is to be used in mixing parameterization (e.g., Henyey et al. 1986; Polzin et al. 1995).

With the existence of these GM spectra, one may ask how well these internal wave spectral models perform in the marginal seas. Levine (2002) presented a modification of GM79 spectrum, which mainly replaces the *e*-folding scale of buoyancy frequency with the Wentzel–Kramers–Brillouin (WKB) scaled thickness of the vertical waveguide; it compared well with the observations from the continental shelf. Another representative internal wave spectral model is the Internal Wave Experiment (IWEX) spectrum (Müller et al. 1978) based on a 42-day internal wave experiment in the Sargasso Sea, which considered contamination of internal waves by fine structures of temperature and currents and by noise in current measurements. Although this model improved adaptability, it is too complicated to use in practice when compared with the GM spectra, since it requires 12 parameters, which include anisotropic parameters, at each frequency. The background fields vary significantly from one marginal sea to another, and hence the parameters employed in internal wave spectrum should be adjusted accordingly.

The South China Sea (SCS) is the largest marginal sea of the Pacific, where performance of internal wave spectrum has also been investigated. There are complex topography and dynamic processes with different scales in the SCS, such as large-scale currents (e.g., Wang et al. 2014, 2016), mesoscale eddies (e.g., Chen et al. 2011; Wang et al. 2012), and so on, which may influence internal wave field and result in a unique internal wave spectrum. Although it is much easier to obtain internal wave frequency spectrum than wavenumber spectrum, studies on internal wave frequency spectrum in the SCS are still rare because of the requirement for high-temporal-resolution sampling. In the northern SCS, internal waves have geographical features. The ratio of total internal wave energy to the GM spectral energy ranges from 2 to 13 for different topographic regions, including the region west of the Luzon Strait, near the Dongsha Island, at the shelf break, and on the continental shelf (Lien et al. 2005). Both Shang et al. (2009) and Xie et al. (2010) pointed out that the internal wave frequency spectrum in the northern SCS is affected by the level of nonlinear interactions among internal waves.

Most studies on internal wave wavenumber spectrum in the northern SCS are based on images of seismic reflection and synthetic aperture radar (SAR). Using seismic measurements with a horizontal resolution of 6.5 m, Dong et al. (2009) presented a 1D horizontal wavenumber spectrum in the northeastern SCS, which was generally consistent with GM76 spectrum. Using SAR images, Liu et al. (2005) provided 2D horizontal wavenumber spectra and examined wave propagation, wavelength, and amplitude. Images of seismic reflection and SAR provide an alternative technique to consider horizontal wavenumber spectrum of internal wave simply; however, the information in frequency and vertical wavenumber domains are missing.

To an extent, these studies have undoubtedly improved our understanding of internal wave spectrum in global ocean, and specifically some marginal seas. However, they mostly provided 1D or 2D spectrum, such as 1D frequency spectrum, 1D horizontal wavenumber spectrum, or 2D horizontal wavenumber spectrum, and failed to provide a 4D spectrum and hence a detailed picture of internal wave field. In this paper, we provide a method to estimate 4D internal wave spectrum. Using this method, a 4D internal wave spectrum was obtained for the first time based on the mooring measurements collected during the SCS Internal Wave Experiment in July 2014. The paper is organized as follows. The field experiment and the method used to derive the 4D internal wave spectrum are given in section 2. Validation and analysis of the 4D internal wave spectrum are presented in section 3, followed by a summary in section 4.

## 2. Data and method

### a. SCS Internal Wave Experiment

To investigate the internal wave feature in the northern SCS, the SCS Internal Wave Experiment was conducted near Dongsha Island during 24–28 July 2014. Four bottom-anchored moorings were deployed around (21°N, 118°E), with three of them (A1, A2, and A4) forming a triangle and the remaining one (A3) at the center of the triangle (Fig. 1). The horizontal spacing between two adjacent moorings ranged from 2.9 to 8.2 km. Water depths of moorings A1–A4 were 984, 943, 999, and 959 m, respectively. Each mooring was equipped with a full-depth thermistor chain consisting of six CTDs and 53–54 thermometers mounted at different depths (Table 1). The depth interval was 10 or 20 m, and the sampling interval was 6 s, which was adequate to capture high-frequency internal waves, since the maximum buoyancy frequency around the thermocline was 2 × 10^{−2} s^{−1}. More information on these moorings is given in Table 1.

The temperature information at the four moorings is shown in Fig. 2, from which some internal solitary waves can be identified clearly in the upper layer. The A waves (consisting of rank-ordered packets) and B waves (usually having a single large wave) both arrived diurnally, which was attributable to variability in the barotropic tidal currents in the Luzon Strait (e.g., Alford et al. 2010; Ramp et al. 2010; Huang et al. 2014). Although the soliton structure in time is likely to appear in frequency spectra as a series of harmonics extending out to high frequencies, the discrete arrivals of the coherent solitons are not well captured by the spectra (because they throw out phase information). Thus, the spectral results to be discussed below did not include the internal solitary waves. Generally, internal solitary waves can be effectively analyzed by tracking the arrival times of specific peaks at different moorings. For example, on 27 July 2014, an internal solitary wave was observed in sequence by mooring A4 at 0448 LST, by mooring A3 at 0504 LST, and by mooring A1 at 0534 LST, which implies a propagation path from A4 to A3 to A1, namely, from southeast to northwest. This result agrees with previous reports (e.g., Liu et al. 2004; Chang et al. 2006; Zhao and Alford 2006; Li et al. 2008; Huang et al. 2017).

Because of the lack of salinity measurements, the only feasible way is to use isothermal displacement as a proxy for isopycnal displacement in subsequent analysis of internal wave spectrum. This is reliable because the temperature-induced density variation ($\u2212\rho 0\alpha \Delta T$), which was calculated based on a linear equation of state, accounted for most of the whole density variation caused by both temperature and salinity, only with a small root-mean-square (RMS) error of 0.01–0.18 kg m^{−3} (Fig. 3).

The isothermal displacement profiles *ξ*(*z*, *t*) were derived from $\xi \u2061(z,t)=z\u2061(T,t)\u2212\u2329z\u2061(T)\u232at$ following Fjeldstad (1933) and Krauss (1966), where *z*(*T*, *t*) is the depth of isothermal at different moments, and $\u2329z\u2061(T)\u232at$ is the time mean of *z*(*T*, *t*) for each mooring. To remove the influence of variable stratification on the vertical structure of internal waves, the isothermal displacement was WKB-scaled $\xi ^=\xi N\u2061(z)/N\xaf$ and the depth coordinate was WKB-stretched $z^=\u222bz0[N\u2061(z\u2032)/N\xaf]dz\u2032$ (e.g., Leaman and Sanford 1975; Althaus et al. 2003). Here, the caret denotes the WKB-scaled value. Because of the sparse salinity measurements, the buoyancy frequency profile *N*(*z*) is a monthly-mean estimation in July 2014 based on the *World Ocean Atlas* (*WOA*) dataset (https://www.nodc.noaa.gov/OC5/woa13/woa13data.html), and $N\xaf$ is the depth-averaged value of *N*(*z*). This is reliable since the atlas-based buoyancy profile matches well with the CTD-derived buoyancy frequency data (Fig. 7a). The WKB scaling aims to convert the ocean to one of constant stratification in which vertical standing modes are sinusoids. Then each WKB-scaled isothermal displacement profile was linearly interpolated onto uniform 5-m vertical interval in the WKB-stretched depth coordinate. The isothermal displacements at the four moorings show a dominant diurnal cycle below 700 m and a dominant semidiurnal cycle above 700 m with their maximum amplitudes both exceeding 100 m (Fig. 4). Hereinafter, only the synchronous measurements from the four moorings were used to estimate 4D internal wave spectrum according to the requirement of method described next.

### b. Method

Because of the limited sampling locations in the horizontal (only four moorings), a method for estimating frequency–wavenumber spectrum called the maximum-likelihood method (MLM; Capon 1969) was used to estimate 4D internal wave spectrum, rather than the conventional spectrum estimation method based on direct Fourier transform. However, there will be a spatial aliasing if the original measurements are directly used in the MLM, because there are also signals with wavelength smaller than the array spacing. To solve this problem, the isothermal displacement was preprocessed by modal harmonic decomposition (Buijsman et al. 2010), which separated different frequency/vertical mode components from the total field.

#### 1) Modal harmonic decomposition

The modal harmonic decomposition consists of three steps. The first step is harmonic decomposition. WKB-scaled isothermal displacement during the whole period is fitted to $\xi ^\u2061(x,z^,t)=\u2211mAm\u2061(x,z^)\u2009cos\u2061(\omega mt)+Bm\u2061(x,z^)\u2009sin\u2061(\omega mt)$ using Fourier transform. Here, $\u2061(x,z^)$ means the horizontal location and WKB-stretched depth of the isothermal displacement. The *ω*_{m} is the frequency set deriving from the Fourier transform. The corresponding amplitude and phase of isothermal displacement at frequency *ω*_{m} are $Am2+Bm2$ and $arctan\u2061(Bm/Am)$, respectively. Figure 5 shows the amplitudes $Am2+Bm2$, which can easily distinguish the near-inertial waves, diurnal, semidiurnal internal tides and their harmonics. At the frequencies lower than 10 cpd, amplitudes were larger than *O*(1) m; especially, large amplitudes of *O*(10) m were mainly found at diurnal frequency below 700 m and at semidiurnal frequency within nearly the whole water column. Amplitudes at frequencies between 10 and ~80 cpd were *O*(0.1) m, and they were even smaller for frequencies higher than 100 cpd.

The second step is modal decomposition. The internal wave has the character of a standing wave in the vertical, while it propagates in the horizontal, since the ocean’s boundaries in the vertical, bottom and surface naturally confine the propagation of internal waves (Gerkema and Zimmerman 2008). This property is exploited in the method of vertical modal decomposition. To find out whether this method is suitable for our observation data, the frequency–vertical wavenumber spectrum of WKB-scaled vertical displacement was presented (Fig. 6), which indicated that although the component with downward-propagating phase is slightly larger than that with upward-propagating phase, especially for high-frequency waves, their differences are still small, not exceeding 0.08%. That is, the downward-propagating energy and upward-propagating energy are almost equal; namely, the assumption of vertically symmetric internal wave field is valid in our observation area and the modal decomposition method can be applied to our data. A similar conclusion can also be drawn from the four-quadrant frequency–vertical wavenumber spectrum of WKB-scaled horizontal baroclinic velocity (Fig. S1 in the online supplemental material).

Then the vertical modes *ψ*(*z*) can be obtained by solving the hydrostatic Sturm–Liouville equation (Gill 1982):

and

Equation (2) assumes that vertical velocities at both sea surface and bottom are zero. The *N*(*z*) is the monthly-mean buoyancy frequency profile in July 2014 based on the *WOA* dataset (Fig. 7a). The subscript *n* indicates the vertical mode number, corresponding to vertical wavenumber of $kz\u2061(n)=(n\pi /b)\u2061(N\xaf/N0)$ (Munk 1981) in the WKB sense, where *b* = 1300 m, *N*_{0} = 5.2 × 10^{−3} rad s^{−1}, and $N\xaf$ is the depth-averaged value of *N*(*z*). The phase speed of *n*th*-*mode internal waves with a frequency of *ω* can be determined using $cp=\u2061(\omega /\omega 2\u2212f2)cn$. The *ψ*_{n}(*z*) is the vertical structural function of *n*th-mode internal waves, which is sinusoid in the WKB-stretched depth coordinate.

In this study, the vertical structural functions of first 50 modes at each mooring were calculated based on its actual water depth, namely, 984, 943, 999, and 959 m for mooring A1–A4, respectively. Results of the first five modes are shown in a WKB-stretched coordinate (Figs. 7b–e). The resulting vertical wavenumbers of *n*th-mode internal waves at moorings A1, A2, and A4 only have slight differences (within a factor of 1.01–1.06) from that at mooring A3. To make the horizontal wavenumber calculation (described below) simple, the vertical wavenumbers at the four moorings were uniformly set to the vertical wavenumber obtained at mooring A3.

The third step is modal harmonic decomposition. Amplitudes for the internal wave components within the internal wave frequency band $f<\omega m<N\xaf$ are fitted to $Am\u2061(x,z^)=\u2211nAmn\u2061(x)\psi n\u2061(z^)$ and $Bm\u2061(x,z^)=\u2211nBmn\u2061(x)\psi n\u2061(z^)$ using the first 50 modes by least squares method, respectively. The amplitudes $Amn2+Bmn2$ of isothermal displacement with different frequencies and mode numbers are shown in Fig. 8. The signal was dominantly distributed in tidal frequencies for the first 20 modes. As a result, the isothermal displacement of the *n*th mode with frequency *ω*_{m} at location $\u2061(x,z^)$ can be reconstructed as

#### 2) Maximum-likelihood method

Given an observation array consisting of *M* sensors (not in a straight line), the *i*th sensor is situated at the vector location $xi$ about an arbitrary origin. Data series of each sensor are divided into *L* segments, with *L* larger than the number of sensors *M* but not too large to ensure enough frequency resolution; we set *L* = *M* + 1 in this study. The cross-power spectrum between the *i*th and *j*th sensors at frequency *ω* is determined by $Gij\u2061(\omega )=1/L\u2211l=1LSil\u2061(\omega )Sjl*\u2061(\omega )$, where *S*_{il}(*ω*) is Fourier transform of the data from the *i*th sensor and the *l*th segment, and the superscript * denotes complex conjugate. The directional information of the input signal is contained in this cross spectrum. However, this cross spectrum contains noise, and we must find a way to extract useful signal (directional information) from it.

The high-resolution frequency–wavenumber spectrum is defined as follows, according to Capon (1969) and Tokimatsu et al. (1992):

where **k** is the 3D wavenumber vector (in rad m^{−1}), and the weights

$qij\u2061(\omega ,k)$ is the inverse of matrix ${exp[ik\u22c5(xi\u2212xj)]Gij\u2061(\omega )}$. The MLM can be regarded as a filter in wavenumber space, with the wavenumber response being expressed as $B\u2061(\omega ,k,k0)=\u2211j=1MAj\u2061(\omega ,k0)\u2009exp\u2061[i\u2061(k\u2212k0)\u22c5xj]$, whose shape and sidelobe structure are not only determined by array geometry but also different for each wavenumber **k**_{0} at which an estimate is obtained. The MLM can pass traveling wave with a wavenumber of **k**_{0} without being distorted and maximally reject traveling waves with wavenumbers other than **k**_{0}. Therefore, the MLM is useful for the estimation of frequency–wavenumber spectrum when the incoherent noise power is relatively small compared to the power of the traveling waves. In addition, to get reliable direction information, the minimum signal wavelength should be larger than the minimum array spacing; otherwise, a white azimuthal spectrum will appear because of spatial aliasing, and the maximum signal wavelength is about 6–7 times (Satoh et al. 2001) or 8 times (Wu and Huang 2015) the maximum array spacing.

The MLM has a higher wavenumber resolution compared to the conventional spectrum estimation based on Fourier analysis. Generally, it is used to derive spectrum in (*ω*, *k*_{x}, *k*_{y}) domain when location differences of the sensors exist only in the horizontal. For example, it has been widely used in seismology (e.g., Kocaoğlu and Fırtana 2011; Leong and Aung 2012) and oceanography (e.g., Drennan et al. 2003; Young and Babanin 2006; Donelan et al. 2015).

Before being used to derive the 4D internal wave spectrum, the MLM was tested first by examining the spectrum of a given synthetic signal. The signal was generated in the form of $\zeta \u2061(x,y,z,t)=A\u2009sin\u2061(k0zz)\u2009sin\u2061(k0xx+k0yy\u2212\omega 0t)$ plus a random noise, which is consistent with the characteristics of vertical standing waves. Here, signal amplitude *A* was set to 10 m, (*x*_{i}, *y*_{i}) (*i* = 1, 2, 3, 4) represented the same array locations as the four moorings A1–A4, and *z*_{i} (*i* = 1, 2, 3, 4) was set to the same depth value that needs to satisfy $sin\u2061(k0zz)\u22600$. Wavenumbers were set to *k*_{0x} = 2*π* × 10^{−5} rad m^{−1}, *k*_{0y} = −2*π* × 10^{−4} rad m^{−1}, and *k*_{0z} = 2*π* × 10^{−2} rad m^{−1}. As a result, the signal has a horizontal wavelength ~10 km, which is larger than the minimum array spacing (2.9 km). The frequency was calculated using $\omega 0=[\u2061(k0x2+k0y2)N2+k0z2f2]/\u2061(k0x2+k0y2+k0z2)$ with *N* = 4 × 10^{−2} rad s^{−1} and *f* = 5.2 × 10^{−5} rad s^{−1}. Using the MLM, spectrum $P\u2061(\omega 0,kh\u2009cos\u2061\theta ,kh\u2009sin\u2061\theta ,k0z)$ can be estimated, where $kh=k0x2+k0y2$ and *θ* ranges from 0° to 360° with an interval of Δ*θ* = 0.25°. The spectral result showed a significant peak at the given wavenumbers (*k*_{0x}, *k*_{0y}), indicating the signal propagates along −84.25°, with a phase speed of 0.64 m s^{−1} and a wavelength of 9.95 km (Fig. 9a). The cases with different angular resolutions showed that as Δ*θ* value decreases, the spectral integration along *θ* only experiences a slight change (Fig. 9c), but the directional accuracy will be largely improved (Fig. 9b) and the distinction between the real peak and spurious peak in the opposite direction becomes larger when Δ*θ* is smaller than 0.5° (Fig. 9d). However, considering that a too small Δ*θ* value will lead to a substantial increase in computational cost, a value of Δ*θ* = 0.25° is chosen here. Besides, the cases with different random noise levels showed similar patterns except for a slight discrepancy in spectral level (Fig. 9e). Furthermore, by changing values of (*k*_{0x}, *k*_{0y}) to *k*_{0x} = −5.34 × 10^{−5} rad m^{−1}, *k*_{0y} = −3.14 × 10^{−5} rad m^{−1} (Fig. 10), we found the MLM still worked well when the signal had a horizontal wavelength ~100 km, exceeding 8 times the maximum array spacing (8.2 km). These test cases indicate the MLM is capable of producing the valid horizontal wavenumber spectrum for signal in the form of $A\u2009sin\u2061(kzz)\u2009sin\u2061(kxx+kyy\u2212\omega t)$.

In the step of modal harmonic decomposition, the WKB-scaled isothermal displacements were separated into components with different frequencies and vertical wavenumbers. For the component at *ω*_{m} and $kz\u2061(n)$, its reconstructed displacements $\zeta ^mn\u2061(xi,z^i,t)$ (*i* = 1, 2, 3, 4) plus a random noise (~1% of the signal amplitude) was used in the MLM to obtain its horizontal wavenumber spectrum $S\u2061(kx\u2061(mn),ky\u2061(mn),kz\u2061(n),\omega m)$ or azimuthal spectrum $S\u2061(\theta ,kh\u2061(mn),kz\u2061(n),\omega m)$. Here, $kx\u2061(mn)=kh\u2061(mn)cos\u2061(\theta )$, $ky\u2061(mn)=kh\u2061(mn)sin\u2061(\theta )$, with *θ* ranging from 0° to 360° with an interval of 0.25*°* and horizontal wavenumber being set to $kh\u2061(mn)=kz\u2061(n)\u2061(\omega m2\u2212f2)/\u2061(N\xaf2\u2212\omega m2)$ according to the dispersion relationship of internal waves. The $z^i$ (*i* = 1, 2, 3, 4) was set to the same depth, the spectra in all layers were estimated and then depth averaged. Only the components with frequency greater than or equal to 2 cpd were used to obtain their horizontal wavenumber spectra because our data lasting only about four days were not long enough to effectively analyze the signal of diurnal frequency after being divided into five segments (as mentioned above). Combining the horizontal wavenumber spectra for all components, the 4D internal wave spectrum $S\u2061(kx,ky,kz,\omega )$ or $S\u2061(\theta ,kh,kz,\omega )$ can be obtained finally. However, it is notable that for the component with a horizontal wavelength smaller than the minimum array spacing, the MLM may give a white azimuthal spectrum because of the aliasing caused by sparse array spacing. In such a case, only the integral of spectral energy along *θ*, not the distribution of spectral energy along *θ*, makes sense.

## 3. Validation and analysis of results

### a. Validation of 4D internal wave spectrum

Based on the Fourier analysis, the isothermal displacements (without modal harmonic decomposition) can be used to estimate 2D frequency–vertical wavenumber spectrum *S*(*k*_{z}, *ω*)_{FFT}, as well as 1D frequency spectrum *S*(*ω*)_{FFT} and 1D vertical wavenumber spectrum *S*(*k*_{z})_{FFT}. The spectra based on Fourier analysis were referred to as FFT spectra. Taking the FFT spectra as the standards, we modified the GM79 spectra. Then the 4D internal wave spectrum was validated against the FFT spectra and modified GM79 spectra.

#### 1) Modification of GM79 spectrum

Based on the assumption that the internal wave field is horizontally isotropic and vertically symmetric, the spectrum of vertical displacement is defined as $F\zeta \u2061(\omega ,j)=b2N0N\u22121\u2061(\omega 2\u2212f2)\omega \u22122EGM\u2061(\omega ,j)$ according to GM79 spectrum (Munk 1981). The *E*_{GM}(*ω*, *j*) is a dimensionless energy density that is determined by $EGM\u2061(\omega ,j)=E0B\u2061(\omega )H\u2061(j)$, where the modal dependence $H\u2061(j)=\u2061(jt+j*t)\u22121/\u2211j=1\u221e\u2061(jt+j*t)\u22121$ and the frequency dependence $B\u2061(\omega )=(2/\pi )\u2061(f/\omega \omega 2\u2212f2) $$\u2061(f\u2264\omega \u2264N)$ need to satisfy $\u2211j=1\u221eH\u2061(j)=1$ and $\u222bfNB\u2061(\omega )d\omega =1$, respectively. Here, *E*_{0} = 6.3 × 10^{−5} is a nondimensional spectral energy level, *j* is the vertical mode number, *j*_{*} = 3 is a mode scale number, *t* = 2 is the high wavenumber slope, and *b* = 1300 m is the *e*-folding scale of stratification.

Significant differences of energy level and slope between the FFT spectra and GM79 spectra mainly occurred in the vertical wavenumber spectrum (Figs. 11 and 12), based on which the modal dependence of GM79 was modified. Two changes were made. One involved changing the value of *j*_{*} from 3 to 1, based on the higher energy levels of low-mode components of the FFT spectra than those of GM79 spectra (Fig. 11b), with a smaller *j*_{*} indicating lower modes contributing relatively more to the total internal wave energy. Levine (2002) suggested that in the deep ocean *j*_{*} is typically found to be near 3, but in shallow water mode one dominates. Yang and Yoo (1999) found *j*_{*} = 1 could fit the linear internal waves well in shallow water over the continental shelf region of the Mid-Atlantic Bight. And the northern SCS is indeed dominated by energetic internal tides, especially low-mode internal tides. For example, through numerical simulation Jan et al. (2008) indicated the dominance of low vertical modes of baroclinic tides in the northern SCS. The other involved changing the value of high vertical wavenumber slope *t* from 2.0 to 2.5, based on the shape of *S*(*k*_{z})_{FFT} (Fig. 11b). High-wavenumber falloff rate in a specific region is important for both internal wave spectrum itself and its associated processes, such as energy cascade. There were several studies focusing on this issue in different regions, as reviewed by Polzin and Lvov (2011). Here, we provide a value of vertical wavenumber slope as −2.5, which indicates the specific behavior of internal wave field in the northern SCS. Not only the modified GM79 vertical wavenumber spectrum agrees well with the observed one (FFT estimation; Fig. 11b), but also the modified GM79 frequency–vertical wavenumber spectrum is consistent with the observed one, with a significant improvement in reproducing the observed low spectral levels in the regions of high frequency and wavenumber (Fig. 12).

#### 2) Performance of 4D internal wave spectrum

The 2D frequency–vertical wavenumber spectrum *S*(*k*_{z}, *ω*)_{MLM} effectively reproduces the magnitude and shape of observed *S*(*k*_{z}, *ω*)_{FFT} and modified GM79 spectrum *S*(*k*_{z}, *ω*)_{GM79*}, which decreases with increasing frequency and vertical wavenumber (Fig. 12). In addition, both 1D frequency spectrum *S*(*ω*)_{MLM} and vertical wavenumber spectrum *S*(*k*_{z})_{MLM} agree well with the observed ones and modified GM79 spectra, exhibiting the falloff rates of *ω*^{−2} and *k*^{−2.5}, respectively (Figs. 11a,b). Departure occurred at the tidal frequency between these spectra, which is not involved in the GM model. The 1D horizontal wavenumber spectrum *S*(*k*_{h})_{MLM} is also consistent with the modified GM79 spectrum *S*(*k*_{h})_{GM79*}, exhibiting a falloff rate of *k*^{−2} (Fig. 11c). These agreements verify the validity of the 4D internal wave spectrum we obtained.

### b. Analysis of 4D internal wave spectrum

By integrating the 4D internal wave spectrum *S*(*θ*, *k*_{h}, *k*_{z}, *ω*) along *θ*, the 3D internal wave spectrum *S*(*k*_{h}, *k*_{z}, *ω*) can be obtained, which shows how internal wave spectrum varies with frequency and wavenumbers. This 3D spectrum decreases significantly with increasing frequency, horizontal wavenumber, and vertical wavenumber (Fig. 13), with a falloff rate of *ω*^{−2} beyond tidal frequencies, and with falloff rates of $kh\u22122$ and $kz\u22122.5$ for horizontal and vertical wavenumbers, respectively (Fig. 11).

The 4D internal wave spectrum *S*(*k*_{x}, *k*_{y}, *k*_{z}, *ω*) can also be projected to (*k*_{x}, *k*_{y}) domain for a fixed frequency and vertical wavenumber. Such spectrum can provide information about power, direction, and speed of propagating waves. The spectral peak indicates the horizontal wavenumber and the direction of incoming waves. The corresponding phase speed *c* and wavelength *λ* can be obtained by $c=\omega /kx2+ky2$ and $\lambda =2\pi /kx2+ky2$, respectively. For example, spectrum for the first-mode semidiurnal internal waves indicated a propagation direction along 160° (Fig. 14a), which is consistent with the fact that M_{2} internal tide was generated in the Luzon Strait and propagated westward into the northern SCS (e.g., Simmons et al. 2011; Zhao 2014). The phase speed *c* and wavelength *λ* are about 2.3 m s^{−1} and 101 km, respectively, which are close to the values reported by Buijsman et al. (2010). By comparing spectra of the same mode, the result suggests that the energy of first-mode internal wave was mainly distributed in the frequencies lower than 6 cpd and generally propagated in a zonal direction (Fig. 14b). When comparing spectra of the same frequency, the energy of M_{2} internal tide was mainly contained in the first five modes, with the first three modes propagating in an almost zonal direction (Fig. 14c).

## 4. Summary

The contribution of this study is to provide a method for estimating 4D internal wave spectrum with limited mooring measurements, and examining the features of internal wave field. This method is made up of two parts. One is the modal harmonic decomposition, which separates the isothermal displacement into components with different frequencies and vertical modes/wavenumbers. The other is the MLM, via which the horizontal wavenumber spectrum or azimuthal spectrum of each separated component is obtained and finally the 4D internal wave spectrum can be derived.

Although this method has been showed to be effective based on both synthetic signals and mooring measurements in the SCS, it still has several limitations that need attention. First, in order to ensure that the vertical wavenumbers of *n*th-mode internal wave component at all moorings are approximately equal, the water depths of all moorings should be equal or extremely close. Second, the spectral results did not include the internal solitary waves, since the discrete arrivals of the coherent solitons are not well captured by the spectra (because they throw out phase information). Third, the moorings should be deployed at three locations at least and cannot be arranged in a straight line, since the estimation of horizontal wavenumber spectrum needs direction information, which should be sampled by the different moorings. Fourth, the array spacing should be predesigned according to the smallest wavelength of interest, because the azimuthal information of internal waves with a horizontal wavelength shorter than the array spacing cannot be resolved because of spatial aliasing. Last, the sampling period should be long enough to distinguish the low-frequency internal wave components, since the time series of data will be divided into several segments and the number of divided segments should be larger than that of moorings when using the MLM.

This study will attract more attention and raise awareness of estimating 4D internal wave spectrum using limited mooring observations. It provides guidance for applying this method to other regions where such mooring observations are available.

## Acknowledgments

This work is jointly supported by the National Natural Science Foundation of China (Grant 41576009, 91628302, and 91858203), Provincial Natural Science Foundation of Shandong (Grant ZR2019JQ13), Global Change and Air–Sea Interaction Project (Grants GASI-IPOVAI-01-03, GASI-03-01-01-03, and GASI-IPOVAI-01-02), State Key Laboratory of Tropical Oceanography, South China Sea Institute of Oceanology, Chinese Academy of Science (Project LTO1601), National Key Research and Development Program (Grant 2016YFC1401403), the NSFC-Shandong Joint Fund for Marine Science Research Centers (Grant U1406401), the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant 41521091), and the Research Programs of the Chinese Academy of Sciences (ISEE2018PY05).

## REFERENCES

*Atmosphere-Ocean Dynamics.*Academic Press, 662 pp.

*Interne Wellen*. Gebruder Borntraeger, 248 pp.

*Advances in Engineering Mechanics: Reflections and Outlooks*, A. T. Chwang, M. H. Teng, and D. T. Valentine, Eds., World Scientific, 297–313.

*MODE Hot Line News*, No. 18, Woods Hole Oceanographic Institution, Woods Hole, Massachusetts, 1–6.

*Evolution of Physical Oceanography*, B. A. Warren and C. Wunsch, Eds., MIT Press, 264–291.

_{S}profiling

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JTECH-D-18-0046.s1.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).