Assessment of fine-scale parameterizations at low latitudes of the North Pacific

Fine-scale parameterizations based on shear and stratification are widely used to study the intensity and spatial distribution of turbulent diapycnal mixing in the ocean. Two well-known fine-scale parameterizations, Gregg–Henyey–Polzin (GHP) parameterization and MacKinnon–Gregg (MG) parameterization, are assessed with the full-depth microstructure data obtained in the North Pacific. The GHP parameterization commonly used in the open ocean succeeds in reproducing the dissipation rates over smooth topography but fails to predict the turbulence over rough topography. Failure of GHP parameterization over rough topography is attributed to the deviation of internal wave spectrum from the Garrett–Munk (GM) spectrum. The internal wave field over rough topography is characterized by energetic intermediate-scale and small-scale internal waves that are not described well by the GM model. The MG parameterization that is widely used in coastal environments is found to be successful in reproducing the dissipation rates over both smooth and rough topographies. The efficacy of GHP and MG parameterizations in evaluating the dissipation rates has been assessed. The result indicates that MG parameterization predicts the magnitude and variability of the dissipation rates better than the GHP parameterization.

features where internal wave spectra deviate from GM spectra 16 . Thus using GHP parameterization near topographic features should be cautious and an appropriate fine-scale parameterization for regions where internal wave spectra deviate from GM spectra is needed to better understand the intensity and spatial distribution of turbulent diapycnal mixing.
Recently, MacKinnon and Gregg [2003a] proposed an analytical model, known as the MacKinnon-Gregg (MG) parameterization, by incorporating the internal wave properties over the continental shelf. It is found to be appropriate in coastal environments where the internal wave field might not be described well by the GM model due to the topography interactions, wind stress, and internal solitary waves [17][18][19][20][21][22][23] . For example, MacKinnon and Gregg [2003a;2005b] verified that MG parameterization succeeds in predicting the dissipation on the New England continental shelf where the wave evolution is controlled by wind stress and bottom drag 24,25 . Van der Lee and Umlauf [2011] reported that MG parameterization is applicable in the Bornholm Basin of the southern Baltic Sea where the internal wave field is governed by low-mode near-inertial wave motions. Xie et al. 23 found that MG parameterization can be applied in the central Bay of Biscay where the internal wave field is characterized by large-amplitude internal tides and internal solitary waves. Shang et al. 22 found that MG parameterization can be applied in the upper ocean of the South China Sea.
Though the MG parameterization was proved to be applicable in coastal environments, few studies have been made to assess this parameterization in the open ocean. Internal wave fields over rough topography in the open ocean are strongly influenced by the topography interactions. Features of these internal wave fields are likely to be consistent with those in coastal environments. Thus it is possible that MG parameterization can overcome the defects of GHP parameterization and work well over rough topography in the open ocean. To confirm this possibility, we assess the GHP and MG parameterizations with the full-depth microstructure data obtained in the North Pacific. The assessment would provide a useful reference for researchers on choosing fine-scale parameterizations to explore the intensity and spatial distribution of turbulent diapycnal mixing in the open ocean.

Results
Our Field observations were conducted at seven stations (A1-A7) in the North Pacific on cruises from October to November in 2015 (Fig. 1a). The study region hosts a wide range of topographic features including featureless abyssal basins, flat ridges, and complex trenches. Detail information about the stations, such as observation time, latitude/longitude, water depths, and topographic roughness (Tr), is given in Table 1. Tr represents the variance of bottom height 11 . The ETOPO2v2 (2-Minute Gridded Global Relief Data) were used to evaluate Tr in a 1/2° × 1/2° box which is approximately a 55 km × 55 km domain. Here the smooth topography is defined as regions with Tr less than 10 5 m 2 , and the rough topography as that with value larger than 10 5 m 2 . Stations A1-A6 were located over smooth topography where values of Tr range from 10 3 to 10 5 m 2 , and station A7 was located over rough topography where the value of Tr is 1.96 × 10 6 m 2 (Fig. 1b).
Profiles of dissipation rate estimated from microstructure data (ε OB ), GHP parameterization (ε GHP ), and MG parameterization (ε MG ) are shown in Fig. 2. To be consistent with the resolution of ε GHP , ε OB and ε MG profiles have been broken into half-overlapping 300-m-long segments (starting from the bottom) and averaged onto the 150-m grid. ε OB decease with increasing depth from O (10 −9 W kg −1 ) in the upper layer to O (10 −11 W kg −1 ) in the deep layer. Both GHP and MG parameterizations succeed in predicting the decreasing trend of the dissipation rates at stations A1-A6 ( Fig. 2a-f). However, profiles of ε GHP and ε MG from station A7 display different patterns (Fig. 2g). ε MG show a decreasing trend with increasing depth as ε OB while ε GHP trend to slightly increase with increasing depth below 500 m where ε GHP deviate from the observed values. To evaluate the success of GHP and MG parameterizations in predicting the magnitude of the dissipation rates, we consider relation between parameterized and observed dissipation rates (Fig. 3a,b). If parameterizations succeed in predicting the magnitude of the dissipation rates, the relation between parameterized and observed dissipation rates will be close to the one-to-one relations of log 10 (ε OB ) = log 10 (ε GHP ) and log 10 (ε OB ) = log 10 (ε MG ). ε GHP from stations A1-A6 (Fig. 3a) show a good relation to ε OB , like the one-to-one relation. However, ε GHP from station A7 show a bad relation to ε OB with most of ε GHP deviating largely from the one-to-one relation. These observations suggest that GHP parameterization predicts the magnitude of the dissipation rates at stations A1-A6 better than that at station A7. For MG parameterization (Fig. 3b), ε MG from all of the stations (A1-A7) show a good relation to ε OB with all of the values gathering around the one-to-one relation. Comparing Fig. 3a,b, one can see that MG parameterization predicts the magnitude of the dissipation rates better than the GHP parameterization.
To assess the efficacy of GHP and MG parameterizations in reproducing the variability of the dissipation rates, we calculate the coefficient of determination (R 2 ). R 2 is the ratio of the difference between the variance of the observed values and the variance of the residuals from the parameterization to the variance of the observed values 26 . It can be interpreted as the proportion of the variability of the observed values that can be explained by the parameterized values. R 2 for each station are shown in Fig. 3c. At stations A1-A6, values of R MG 2 and R GHP 2 are positive and larger than 0.5, which indicates that both GHP and MG parameterizations can effectively predict the variability of the dissipation rates at stations A1-A6. At station A7, R MG 2 is still positive and larger than 0.5 while R GHP 2 is negative. A negative R GHP 2 suggests that GHP parameterization fails to predict the variability of dissipation rates at station A7. For all of the stations, the values of R MG 2 are larger than that of R GHP 2 , which implies that MG parameterization predicts the variability of the dissipation rates better than the GHP parameterization.
The above analysis indicates that MG parameterization succeeds in predicting the dissipations at all of the stations while GHP parameterization only succeeds at stations A1-A6. The success of fine-scale parameterizations greatly depends on the internal wave fields. To explore the feature of the internal wave field at each station, baroclinic velocity was formally decomposed onto a set of orthogonal vertical modes. The vertical structure of each mode is governed by 24,27 where c j is the separation constant (eigenvalue), and H is the water depth. The vertical mode shapes are calculated from numerical solution of equation (1) using stratification profile. Vertical velocity and vertical displacement associated with each mode are proportional to Ψ j , while the horizontal velocity is proportional to dΨ j /dz. Baroclinic velocity data were fit to the baroclinic modes using a least square regression. Figure 4a shows a sample stratification (N 2 ) profile from station A4. Stratification is strong in the upper 300 m (N 2 ~ 10 −4 s −2 ) and become weak below (N 2 ~ 10 −6 -10 −5 s −2 ). Figure 4b-g show the first six baroclinic modes calculated with the stratification profile in Fig. 4a. The modes have been normalized. The velocity fitted from the first six modes is shown in Fig. 4h. The modal fit captures the dominant low-mode signal well, though it does not reproduce the small-scale fluctuations. To look at the distribution of horizontal kinetic energy in modes, we calculate the ratio (Ra) of the energy in low modes (modes 1-6) to the total energy. The ratio of the energy in high modes (modes >6) to the total energy is given as 1-Ra. The result is shown in Fig. 4i. Note that background flows associated with mean current or mesoscale eddies may also make important contribution to the observed horizontal velocity. Here we remove the background flows by removing the depth-mean of each velocity component. This method works well though it might not remove the background flows completely. As one can see from Fig. 5 that large-scale (>1000 m) motions are roughly consistent with the GM spectrum. The internal wave fields are dominated by low-mode internal waves at stations A1-A6 with the energy in the first six modes accounting for more than half of the total energy. The Ratios even exceed 70% at stations A5 and A6. However, a different energy distribution is found at station A7. Low-mode internal waves are no longer the dominant components. Only 22.7% energy are captured by the first six modes, and 77.3% energy are from high modes. These observations indicate that high-mode internal waves are more active at station A7 than other stations.
The above modal analysis indicates that internal wave fields at different stations show different features. The internal wave field at station A7 is characterized by energetic high-mode internal waves. To find out whether the internal wave fields can be described well by the GM model, we compare the observed spectra with the GM where E is the non-dimensional energy parameter; b is the vertical length scale; N 0 is the buoyancy frequency scale; and j is the mode number. The vertical wavenumber is given as The observed and GM dropped spectra are shown in Fig. 5. The observed dropped spectra at stations A1-A6 (Fig. 5a-f) roll off smoothly as − k z 2 as predicted by GM model (equation 2). However, observation at station A7 (Fig. 5g) did not show a smooth decrease with vertical wavenumber as predicted by GM model. Instead, it shows some significant elevated peaks at wavenumber band of 10 −3 -10 −2 cpm (indicated by the arrows). These elevated peaks correspond to the energetic intermediate-scale internal waves which cannot be described well by the GM model. In addition, the  energy level at higher wavenumbers (k z > 10 −2 cpm, small-scale internal waves) is higher than the GM dropped spectrum, which indicates that more energy is transferred from large scales to small scales than that predicted by the GM model.

Discussion
With the full-depth microstructure data obtained in the North Pacific, we have assessed the GHP and MG parameterizations. GHP parameterization succeeds in predicting the dissipation over smooth topography (stations A1-A6) but fails over rough topography (station A7). One possible explanation for the failure of GHP parameterization over rough topography is that the internal wave spectra over rough topography diverge from the GM model. GHP parameterization was developed based on the assumption that the waves are statistically stationary with respect to wave-wave interactions by which energy is transferred from large to small scales 7 . It is typically evaluated for the internal wave field with GM spectral shape, but inappropriate for the regions where the internal wave field deviates from the GM model 6,15,16 . Station A7 is characterized by energetic intermediate-scale and small-scale internal waves that are not described well by the GM model (Fig. 5). These energetic internal waves might be caused by topography interactions rather than the wave-wave interactions as predicted by GM model. Studies [29][30][31][32] have indicated that topography interactions could push energy directly into high wavenumbers and account for the departure of the observed spectrum from the wave-wave interactions underlying the GM model. Station A7 is located at the intersection of the Yap Trench and Mariana Trench (Fig. 1) where the topography is rather rough (Tr = 1.96 × 10 6 m 2 ). Energetic high-mode internal waves (Fig. 4i), elevated peaks at wavenumber band of 10 −3 -10 −2 cpm, and high energy level at large wavenumbers (>10 −2 cpm) in the spectrum of station A7 (Fig. 5g) might result from topography interactions.
In the classical lee-wave problem 33 , a constant flow (e.g., geostrophic flow) in a stratified fluid over a variable bottom topography generates upward radiating internal waves. The vertical wave number of the generated waves is related to the horizontal wave number (k h ) of topography 32 ,  (Fig. 5g). In addition to constant flow, strong tidal currents (e.g., barotropic diurnal and semidiurnal tides) over topographic features with widths less than a tidal excursion can also cause internal waves, which can propagate away from the topography as the forcing frequency or high-frequency (superharmonics) internal waves 31,34,35 . The dispersion relationship for internal waves is where ω is the frequency of internal waves. The observed area is characterized by strong diurnal and semidiurnal tidal currents 36,37 . Thus for generated internal waves with diurnal (ω ≈ 1.16 × 10 −5 s −1 ) and semidiurnal (ω ≈ 2.31 × 10 −5 s −1 ) frequencies, k z is 10 −3 -10 −2 cpm, which is consistent with the observed internal waves at wavenumber band of 10 −3 -10 −2 cpm in the spectrum of station A7 (Fig. 5g). Other possible reasons for the failure of GHP parameterization over rough topography are wave-mean interactions. If the wave-mean interactions dominate nonlinear transports, significant biases are possible 38 . In order to gauge the extent to which wave-mean interactions may be significant, it is convenient to evaluate the ratio of time scales characterizing wave-wave and wave-mean interactions. Invoking the ray tracing equations and assuming waves are randomly aligned with the mean shear, the ratio can be shown to be 38 nlin wm z c z where S represents the mean shear, k z c represents the cutoff wave number. Equation (5) states that, if the waves are randomly aligned with the mean shear, wave-mean interactions dominate the spectral energy transport in vertical wave number space at k z c for mean shears in excess of N/8. Figure 6 shows the ratio for vertical wave numbers equal to k z c . Although many of the values from stations A1-A6 are greater than one, GHP parameterization is effective at these stations. It might be understood that the finescale wavefield tends to be aligned normal to the mean shears. These phenomena also occurred in other observations 6,39 , in which the ratio exceeded one and finescale parameterizations were still applicable. Here, it is necessary to compare the ratio between stations A1-A6 and station A7. Ratio of station A7 is O(1) at most of depth and there is no significant discrepancy among the stations, which suggests that wave-mean interactions are not the major cause of the failure of GHP parameterization at station A7.
Although GHP parameterization fails over rough topography where the internal wave field deviates from the GM model, our observation indicates that the MG parameterization succeeds in reproducing the dissipation not only over smooth topography but also over rough topography. In addition, our statistical results indicated that the MG parameterization predicts the magnitude and variability of the dissipation rates much better than the GHP parameterization. MG parameterization was developed based on the assumption that there is no statistical relationship between shear in low-and high-mode waves, and the strength of low-mode shear is decoupled from properties of high-mode waves 18,19 . MG parameterization was first proposed to predict the dissipation rates over the New England shelf 18 . Later it is found to be suitable for the coastal environments [19][20][21] where the internal wave fields are more complex than that predicted by the GM model. Here we further verify that MG parameterization is also applicable in the deep ocean where the internal wave field deviates from the GM model.
SCIenTIfIC RepoRts | (2018) 8:10281 | DOI:10.1038/s41598-018-28554-z GHP parameterization was most commonly used to explore the intensity and spatial distribution of turbulent diapycnal mixing in the open ocean. With the GHP parameterization, elevated turbulent diapycnal mixing (κ ≥ 10 −4 m 2 s −1 ) have been found over rough topography [8][9][10][11][12][13][14] . This study implies that dissipation over rough topography predicted by GHP parameterization might be inaccurate, and the MG parameterization can overcome the defects of GHP parameterization and allow for more accurate intensity and spatial distribution of turbulent diapycnal mixing in the open ocean. One challenge of the application of these fine-scale parameterizations is the necessary local adjustment of ε 0 . In spite of this, our result shows that the range of ε 0 for MG parameterization is smaller than that for GHP parameterization (0.75 × 10 −9 ε ≤ ≤ 10.47 × 10 −10 W kg −1 ). A small range of ε 0 for MG parameterization implies that using the MG parameterization would lead to a smaller distortion in the broad mapping of dissipation rates than using the GHP parameterization if an averaged ε 0 were adopted in the parameterizations. Yet our observation only covers a small region of the Pacific. In order to complete the assessment of fine-scale parameterizations, more extensive simultaneous measurements of fine-structure and microstructure are needed in the near future.

Methods
Data. Microstructure data were obtained with an expendable Vertical Microstructure Profiler (VMP-X), a free-falling instrument ballasted to fall at 0.6-0.8 m s −1 . The instrument was equipped with standard turbulence sensors for measuring turbulent velocity shear (∂u/∂z and ∂v/∂z), temperature, and pressure 40 . The maximum measurement depth is 6000 m, which satisfied the maximum depths of the stations (Fig. 1a). Fine-structure temperature and salinity data were collected using a Seabird 9-11 Plus CTD. The CTD data were processed according to standard procedures as recommended by the manufacturer, and bin averaged to 2-m resolution. Buoyancy frequency was computed from finite differencing with the CTD temperature and salinity data. Fine-structure velocity data was measured with two lowered acoustic Doppler current profilers (an upward 300 kHz LADCP and a downward 300 kHz LADCP) mounted on the CTD package. Vertical bin size of LADCPs was 10 m and there were 20 bins for each LADCP. Velocity profiles were processed on a 10-m depth grid with software developed at the Lamont-Doherty Earth Observatory, Columbia University. Baroclinic velocity was computed by removing the depth-mean of each velocity component. Shear was calculated by first-differencing velocity over 10-m intervals.
OB k k 2 1 2 where ν is the kinematic viscosity and < >denotes the spatial average. The lower integration limit, k 1 , is set to 1 cpm, and the upper limit, k 2 , is the highest wavenumber not contaminated by vibration noise. Shear spectrum is computed from the shear signal with half-overlapping 6-m depth segments, thus the vertical resolution of ε OB is 3 m. The noise level 40 of the VMP-X is 10 −11 W kg −1 .
MG parameterization. MG dissipation rate (ε MG ) can be expressed in terms of fine-scale shear and stratification as 18 MG 0 0 0 where S 0 = N 0 = 3 cph and ε 0 is an adjustable constant that tunes the parameterized dissipation rate to the observational data. ε MG are estimated with the 10-m shear and buoyancy frequency; thus ε MG have a vertical resolution of 10 m. MG parameterization is found to be not sensitive to the shear resolution. It can also reproduce the observed dissipation with low shear resolution, e.g., 20-m shear, 50-m shear (not shown). ε 0 in equation (7) shows variability in different regions [19][20][21]23 , spanning from 10 −10 to 10 −8 W kg −1 . In order to fit the parameterized dissipations best to the observed dissipations, for each station ε 0 is set to the value that minimizes the sum of squared residuals (a residual being the difference between an observed value and the parameterized value). ε 0 are 0.75 × 10 −9 , 3.12 × 10 −9 , 2.92 × 10 −9 , 3.08 × 10 −9 , 2.01 × 10 −9 , 4.19 × 10 −9 , and 1.89 × 10 −9 W kg −1 for stations A1-A7, respectively, with a mean of 2.57 × 10 −9 W kg −1 .
GHP parameterization. GHP dissipation rate (ε GHP ) depends on fine-scale shear and strain variances as 5,6,8 increasingly blue at higher wavenumbers (k z > 0.140 rad m −1 ) because of instrument noise, so the upper bound for the shear variance integration is set to 0.140 rad m −1 , corresponding to vertical wavelength λ z = 45 m. The lower bound for the shear variance integration is the lowest resolved wavenumber (0.021 rad m −1 ), corresponding to vertical wavelength of 300 m. The strain spectra roll off roughly as − k z 3/2 at the high wavenumbers (k z > 0.419 rad m −1 ). The upper bound for the strain variance integration is set to 0.419 rad m −1 , corresponding to λ z = 15 m. The lower bound for the strain variance integration is set to 0.042 rad m −1 (corresponding to λ z = 150 m) because strain variance at the low wavenumbers (k z < 0.042 rad m −1 ) might be influenced by background stratification. With the obtained V z 2 and ξ z 2 , ε GHP are calculated with equation (8), which have a vertical resolution of 150 m. The reference dissipation rate ε 0 in equation (8) shows variability in different studies 8,[10][11][12][13]16 , spanning from 10 −10 to 10 −9 W kg −1 . In order to fit the parameterized dissipations best to the observed dissipations, for each station ε 0 is set to the value that minimizes the sum of squared residuals. ε 0 are 0.76 × 10 −10 , 3.31 × 10 −10 , 4.68 × 10 −10 , 4.43 × 10 −10 , 4.07 × 10 −10 , 10.47 × 10 −10 , and 0.30 × 10 −10 W kg −1 for stations A1-A7, respectively, with a mean of 4.0 × 10 −10 W kg −1 .
Data availability. The research data can be accessed from the corresponding author Xiao-Dong Shang, whose email is xdshang@scsio.ac.cn.