1 Introduction

ITs are internal gravity waves at tidal frequencies. They are generated when barotropic tides scatter on variable bottom topographies, such as seamounts, mid-oceanic ridges, and continental shelves (Bell, 1975; Munk and Wunsch, 1998). The vertical structure of ITs is modeled by normal modes (Gill, 1982). High-mode ITs dissipate locally near generation sites, whereas low-mode ITs can propagate several thousands of kilometers before they break in regions that are far away from the generation sites (Zhao et al., 2010, 2016; Pickering and Alford, 2012). The breaking of ITs provides an important energy source for diapycnal mixing, which supplies approximately 1 TW of energy for the maintenance of meridional overturning circulation (Garrett, 2003; Whalen et al., 2020). In situ observations indicate that the IT generation sites are hotspots for turbulent mixing. For instance, the maximum value of diapycnal diffusivity exceeds 10⁻³ m² s⁻¹ near the Hawaii Ridge (Klymak et al., 2006). Moreover, the diapycnal diffusivity reaches 3 × 10⁻² m² s⁻¹ in Luzon Strait (Yang et al., 2016) and 0.9 × 10⁻⁴ m² s⁻¹ near the Mariana Trench in the western Pacific Ocean (Shang et al., 2017). Moreover, the mixing induced by ITs considerably affects the redistribution of heat (Jithin and Francis, 2020), nutrients (Sharples et al., 2007), and underwater acoustics (Powell et al., 2013). Therefore, mapping the dissipation of ITs and IT-driven mixing plays an essential role in understanding oceanic energy budget and mass transport. However, mapping the global IT-driven mixing by widespread observations is difficult owing to the small-scale turbulence. Thus, parameterization is an effective way of estimating the IT-driven dissipation and turbulent mixing.

Considering the local dissipation of ITs near generation sites, St. Laurent et al. (2002) proposed parameterization for IT-driven mixing (hereafter referred to as LSJ02). They pointed out that one-third of IT energy dissipates locally near generation sites, and the dissipation rate decays exponentially upward from the seafloor. The LSJ02 parameterization has been widely used in several ocean circulation models, thus improving the simulation of meridional overturning circulation (Jayne, 2009; Melet et al., 2013; Schmittner and Egbert, 2014). In addition, this parameterization has been adopted for estimating tidal mixing in the South China Sea (Wang et al., 2016) and Andaman Sea (Jithin et al., 2020). Nevertheless, several studies have indicated the nonuniform spatial distribution of IT energy's local dissipation efficiency. For instance, the fraction of IT energy dissipating locally can reach 0.36 – 0.60 in generation sites, such as the Luzon Strait (Jan and Chen, 2009; Alford et al., 2011; Kerry et al., 2013), and 1 in semiclosed basins, such as Indonesian Seas (Koch-Larrouy et al., 2007). The vertical distribution of energy dissipation is also crucial. Koch-Larrouy et al. (2007) (hereafter referred to as KL07) found a good agreement between depth-integrated IT energy and buoyancy frequency. Accordingly, they improved the vertical structure function in the LSJ02 parameterization by setting it as a function of stratification. Nevertheless, in the LSJ02 and KL07 parameterizations, only the local dissipations of high-mode ITs are considered, whereas the radiation of IT energy is neglected. Considering the remote dissipation of low-mode ITs, de Lavergne et al. (2020) proposed a new parameterization (hereafter referred to as dL20) based on the energy conservation equation. In particular, they replaced background diffusivities in the LSJ02 parameterization with far-field tidal mixing.

The applicability of these three parameterizations to the northern Pacific Ocean has not been thoroughly evaluated. Furthermore, the map of the bathymetry of the northern Pacific reveals significant spatial variations (Fig. 1). The Luzon Strait, Ryukyu Island Chain, Izu-Bonin Ridge, Mariana Island Arc (Kunze, 2017), and Hawaiian Islands (Zhao et al., 2016) located in this area are considered the strong generation sites of ITs. Meanwhile, the Philippine Sea and open Pacific Basin are the long-range propagation areas of the ITs (Niwa and Hibiya, 2004; Zhao et al., 2010; Kerry et al., 2013; Wang et al., 2018). Therefore, both the local dissipation of high modes and the remote dissipation of low modes are expected in the northern Pacific.

The present study aims to evaluate the performance of the LSJ02, KL07, and dL20 parameterizations in the northern Pacific Ocean. The paper is organized into sections. The data and three parameterizations are introduced in Section 2. The calculation results are shown in Section 3, along with comparisons with oceanic observations. Finally, the conclusions are drawn in Section 4.

2 Data and Methods

The IT-driven dissipation and diapycnal diffusion in the northern Pacific are estimated using the LSJ02, KL07, and dL20 parameterizations, in which the rough bottom topography is key to IT generation. The conversion rate from barotropic tides to baroclinic tides depends on the roughness of topography (Decloedt and Luther, 2010), and the accuracy of topographic roughness is limited by the resolution of bottom topography. Thus, higher resolution bathymetric data are required to calculate the topographic roughness. The 1 arc-minute ETOPO1 is a high resolution global bathymetric dataset that is used currently and can distinguish topographic features with horizontal wavelengths of 10 – 20 km. Herein, the topography is taken from the bathymetric data of ETOPO1. The amplitude and phase of barotropic tides are extracted from the regional solution for the Pacific obtained from the Oregon State University Tidal Inversion Software (OTIS, Egbert and Erofeeva, 2002) with 1/12˚ spatial resolution. Fig. 2 shows the co-tidal charts. Furthermore, the climatological density and buoyancy frequency, showing a spatial resolution of 1/4˚, are obtained from the World Ocean Atlas 2013 (WOA13).

2.1 LSJ02 and KL07 Parameterizations

In LSJ02 parameterization, the energy conversion from barotropic to baroclinic tides is calculated as

$ E(x, y) = \frac{1}{2}{\rho _0}{N_b}k{h^2}\left\langle {u_{bt}^2} \right\rangle, $

(1)

where ρ₀ is the background density, N_b is the buoyancy frequency at the seafloor, and k and h, denote the topographic wavenumber and amplitude, respectively. The wavenumber is set as k = 2π/(10 km) in Jayne and St. Laurent (2001). The topographic roughness h² is defined as the variance of the bathymetry over a 1/3˚ × 1/3˚ domain square (Li and Xu, 2014). In addition, $\left\langle {u_{bt}^2} \right\rangle $ is the tidal period mean of squared barotropic tidal velocity, and u_bt is derived by the OTIS model with 1/12˚ spatial resolution.

The energy dissipation rate ε and diapycnal diffusivity κ_v in LSJ02 parameterization are respectively calculated as

$ \left\{ \begin{array}{l} \varepsilon = \left({q/\rho } \right)E(x, y)F(z) \\ {\kappa _v} = {\mathit{\Gamma}} \varepsilon /{N^2} + {\kappa _0} \\ \end{array} \right., $

(2)

where

$ F(z) = \frac{{\exp [ - (H + z)/\zeta ]}}{{\zeta [1 - \exp (- H/\zeta)]}} . $

(3)

The mixing efficiency Г and local dissipation efficiency q are set to 0.2 and 0.3, respectively. Here κ₀ = 10⁻⁵ m² s⁻¹ is the background diffusivity. In Eq. (3), F(z) is the vertical structure function depending on the bathymetry H, z is the depth from the surface, and the e-folding height is set as ζ = 500 m.

Based on the LSJ02 parameterization, the KL07 parameterization indicates the relationships between the magnitude of the depth-integrated energy E_n and buoyancy frequency N as E_n ~ N (when the vertical gradient is negative) and E_n ~ N² (when the vertical gradient is positive). Accordingly, the KL07 parameterization modifies the vertical structure function, and the corresponding expressions are given below.

$ F(z) = \left\{ \begin{array}{l} \frac{N}{{\int {{N^2}{\text{d}}z} }}, \frac{{{\text{d}}N}}{{{\text{d}}z}} < 0 \\ \frac{{{N^2}}}{{\int {{N^2}{\text{d}}z} }}, \frac{{{\text{d}}N}}{{{\text{d}}z}} > 0 \\ \end{array} \right.. $

(4)

In this study, the direction of the positive z-axis is upward.

2.2 dL20 Parameterization

The dL20 parameterization decomposes the energy dissipation of ITs into four components: attenuation by wave-wave interactions ε_wwi, direct breaking of ITs through shoaling ε_sho, dissipation at critical topographic slopes ε_cri, and scattering by abyssal hills ε_hil. The formulas of these four components are respectively given below.

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _{{\text{wwi}}}} = \frac{{{E_{{\text{wwi}}}}}}{\rho }\frac{{{N^2}}}{{\int {{N^2}{\text{d}}z} }}, {\text{ }}{\varepsilon _{{\text{sho}}}} = \frac{{{E_{{\text{sho}}}}}}{\rho }\frac{N}{{\int {N{\text{d}}z} }}, {\text{ }}{\varepsilon _{{\text{cri}}}} = \frac{{{E_{{\text{cri}}}}}}{\rho }\frac{{\exp (- {h_{{\text{ab}}}}/{H_{{\text{cri}}}})}}{{{H_{{\text{cri}}}}(1 - \exp (- H/{H_{{\text{cri}}}}))}}} \\ {{\varepsilon _{{\text{hil}}}} = \frac{{{E_{{\text{hil}}}}}}{\rho }\left[ {{r_{{\text{bot}}}}\frac{1}{{{{(1 + {h_{{\text{ab}}}}/{H_{{\text{bot}}}})}^2}}}\left({\frac{1}{H} + \frac{1}{{{H_{{\text{bot}}}}}}} \right) + (1 - {r_{{\text{bot}}}})\frac{{{N^2}}}{{\int {{N^2}{\text{d}}z} }}} \right]} \end{array}} \right. . $

(5)

As shown above, E_wwi, E_sho, E_cri, and E_hil are the respective energies of the four components. In addition, h_ab denotes the height above the seafloor, H_cri is the e-folding scale, r_bot is the fraction of E_hil dissipating near the bottom, and H_bot is the height scale that depends on the seafloor (de Lavergne et al., 2019, 2020). The total dissipation rate ε and diapycnal diffusivity κ_v are given below.

$ \left\{ \begin{array}{l} \varepsilon = {\varepsilon _{{\text{wwi}}}} + {\varepsilon _{{\text{sho}}}} + {\varepsilon _{{\text{cri}}}} + {\varepsilon _{{\text{hil}}}} \\ {\kappa _v} = \Gamma \varepsilon /{N^2} \\ \end{array} \right.. $

(6)

In this study, the dissipation rate provided by de Lavergne (2020, https://doi.org/10.17882/73082) is used directly, and the diapycnal diffusivity is computed with Eq. (6).

To evaluate the applicability of the three parameterizations to a real ocean, this study compares the dissipation rates and diapycnal diffusivities obtained by the three parameterizations with the observations from the work of Kunze (2017), which were estimated using fine-scale parameterizations based on hydrographic data (ftp.nwra.com/outgoing/kunze/iwturb.).

3 Results and Analysis 3.1 Spatial Patterns of Dissipation and Diapycnal Diffusion

As shown in Fig. 3, IT energy is horizontally inhomogeneous with spatial variations of 5 – 6 orders of magnitude. The Luzon Strait, Ryukyu Island Chain, Mariana Island Arc, and Hawaiian Islands are the main generation sites of ITs. Energy conversion reaches 10⁻¹ – 10⁰ W m⁻² in regions with topographic roughness higher than 5 × 10⁵ m² and 10⁻⁶ – 10⁻⁵ W m⁻² in regions with smooth topography.

Fig. 4 shows the full-depth averaged dissipation rate and diapycnal diffusivity maps from parameterizations and observations. The results estimated using the LSJ02 and KL07 parameterizations show spatial patterns similar to those of baroclinic energy. Strong dissipation occurs in generation sites. For instance, the dissipation rate exceeds 10⁻⁸ W kg⁻¹ in the Luzon Strait and Ryukyu Island Chain, as well as reaches 10⁻⁸ W kg⁻¹ in the Mariana and Hawaiian Islands. In addition, the maximum diffusivity estimated using the LSJ02 parameterization exceeds 10⁻³ m² s⁻¹. The value of κ_v estimated by the KL07 parameterization is generally smaller than that of the LSJ02 due to the discrepancy of the vertical structure functions. Fig. 4c shows the dissipation rate estimated using the dL20 parameterization. In the generation sites, the results obtained by the three parameterizations are close to the observations (Figs. 4d, h), while in the far-field regions, the patterns of results from the dL20 parameterization match better with the observations. Apart from local dissipation, this scheme estimates the dissipation of low-mode ITs in far-field regions. For example, the dissipation rate and diffusivity reach 10⁻⁸ W kg⁻¹ and 10⁻⁴ m² s⁻¹, respectively, in the Philippine Sea. Kerry et al. (2013) found that 4.78 and 3.82 GW of M₂ IT energy radiate into the Philippine Sea from the Luzon Strait and Mariana Island Arc, respectively, and these can provide sufficient energy for mixing in the Philippine Sea.

At the same time, the dL20 parameterization estimates tidal energy attenuation through wave-wave and wave-topography interactions (Fig. 5), amounting to the respecttive 82% and 18% of the total tidal energy dissipation rates in the northern Pacific. The wave-wave interactions play a major role in low and middle latitudes (5˚ – 30˚N). In addition, the dissipation of low-mode ITs caused by interactions with critical and shoaling slopes mainly occurs in ridges and continental shelves. Furthermore, the decay of low-mode IT energy due to scattering on abyssal hills, particularly in the Philippine Basin, is non-negligible.

In terms of magnitude and horizontal distribution, the results estimated by the three parameterizations vary. To assess precisely the three-dimensional (3D) structure of energy dissipation and diapycnal diffusion, the current study compares the calculated results with the observations (Kunze, 2017) along the 24˚N section. This section crosses the Philippine Sea, Mariana Island Arc, Central Pacific, and Hawaiian seamounts, thus covering the generation sites and far-field regions. As shown in Fig. 6, significant discrepancies exist in the vertical structures of the dissipation rates obtained by the three parameterizations. The vertical structure in the LSJ02 parameterization is a function of the height above the seafloor. The strongest dissipation occurs at the bottom boundary and exponentially decays upward. Nevertheless, the vertical distribution from the KL07 parameterization is similar to N² under the thermocline; that is, the strongest dissipation shifts to the thermocline and gradually decreases with depth (Alberty et al., 2017; Cuypers et al., 2017). The results from the KL07 parameterization match better the observations in the upper layer compared with that from the LSJ02 parameterization. At the same time, the vertical structure in the dL20 parameterization is similar to that in the LSJ02 parameterization near the seafloor and in the KL07 parameterization in the upper ocean.

According to Fig. 7, the dissipation caused by the interaction with critical slopes and abyssal hills is significant near the bottom boundary, with the value being approximately O(10⁻⁹) W kg⁻¹. Additionally, ε_wwi indicates the energy dissipation rate of low modes in the far-field and the high modes near the generation sites through the wave-wave interactions. The dissipation caused by wave-wave interactions is dominant at depths of 0 – 2.5 km, and the maximum value reaches O(10⁻⁸) W kg⁻¹. ε_sho contains the energy dissipation rates from the direct breaking of modes 1 – 10 through shoaling, which mostly acts at the shelf break and shoreward. According to Eqs. (4) and (5), both the vertical structure functions of ε_wwi and ε_sho depend on the buoyancy frequency, which is similar to those of the KL07 parameterization. Therefore, in the interior ocean, the dissipation of high modes estimated by the KL07 parameterization could be enhanced by wave-wave interactions.

In addition, ε_cri implies the dissipation rates from modes 1 – 10 ITs at the critical slopes, which occur primarily at continental slopes. Here, ε_hill is the energy dissipation rate of ITs on the abyssal hill. The dissipation caused by the interaction with critical slopes and abyssal hills is significant near the bottom boundary, with the value being approximately O(10⁻⁹) W kg⁻¹. According to Eqs. (3) and (5), ε_cri decays exponentially upward from the bottom, which shows a similar phenomenon as the result of the LSJ02 parameterization. The vertical structure function of ε_hill consists of two parts: one decaying as a negative power function from the seafloor upwards, and the other is proportional to N². Consequently, the bottom-intensified dissipation estimated by the LSJ02 parameterization could relate to the dissipation driven by the interaction of high modes with critical topography and abyssal hills.

In summary, among the three parameterizations, the pattern estimated by dL20 parameterization exhibits the highest level of consistency with the observations (Fig. 6d), particularly at depths of 0 – 2 km. However, due to the absence of observation, the performance of the three parameterizations in the deep layer is difficult to evaluate. Notably, the dissipation rates estimated by all three parameterizations are greater than the observations at the top of the Mariana and Hawaii ridges.

Fig. 8 shows the diapycnal diffusivity along the 24˚N section. As can be seen, the dissipation rates obtained by the LSJ02 and dL20 parameterizations significantly decay upward from the seafloor. The maximum diapycnal diffusivity estimated by the LSJ02 parameterization reaches 2 × 10⁻² m² s⁻¹ near the rough topography and is greater than that estimated by dL20 parameterization. The diffusivity from the KL07 parameterization generally remains constant below the thermocline, and strong mixing only occurs at the peak of the submarine ridges. In addition, the background diffusivities in the LSJ02 and KL07 parameterizations are set to be constant (10⁻⁵ m² s⁻¹); therefore, the diffusivities are uniform at a distance above 1 km from the seafloor. However, this condition results in overestimations in the far-field regions compared with the observations (Fig. 8d). In the dL20 parameterization, the background diffusivity is replaced by the far-field diffusion driven by low-mode ITs and matches the observations in the open ocean (160˚ – 180˚E in Fig. 8c).

3.2 Applicability of Three Parameterizations

In Section 3.1, the spatial distribution characteristics of IT-driven dissipation and diapycnal diffusion in the Northern Pacific were analyzed. In addition, a preliminary comparison of the results from the parameterizations and observations is made. Next, the calculated results in the four generation sites (areas A1 – A4 in Fig. 3: Luzon Strait, Ryukyu Islands, Ogasawara-Mariana Arc, and Hawaii Islands) and the two far-field regions (areas A5 and A6 in Fig. 3: the Philippine Sea and Kuroshio Extension) are compared with the observations quantitatively in order to accurately assess the applicability of the three parameterizations to a real ocean. Here, a coefficient α is defined as

$ \alpha = \left| {{{\log }_{10}}\left({{\varepsilon _{\text{p}}}/{\varepsilon _{\text{o}}}} \right)} \right|, $

(7)

where ε_p and ε_o are the dissipation rates from the parameterizations and observations, respectively (Chen et al., 2021). The percentage of α < 1 is also calculated, which indicates that the error of ε is lower than one order of magnitude. As the radiation of IT energy is neglected in the LSJ02 and KL07 parameterizations, only the corresponding results in the generation sites are discussed. The results estimated by the dL20 parameterization are compared with the respective observations recorded in the generation sites and far-field regions. The vertical profiles of ε and κ_v obtained by the three parameterizations and observations in the six areas are presented in Fig. 9, while Fig. 10 shows the comparison results of all the stations marked in Fig. 3.

As shown in Figs. 9(a – d, g – j), ε and κ_v estimated by the LSJ02 parameterization are 2 – 3 orders of magnitude larger than the observations within 1 – 2 km above the seafloor. The results with α < 1 only account for 43% of the outcomes from the four generation sites (Fig. 10a). For the KL07 parameterization, most of the estimation results satisfy α < 1 (Fig. 10g). Moreover, 80.7% of the errors between the KL07 parameterization and the observations are lower than one order of magnitude in the four generation sites. The estimated ε and κ_v are highly consistent with the observations in the generation sites, except the Luzon Strait, for which the calculated results are higher than the observations. Furthermore, the percentage of α < 1 is the largest in the Ogasawara-Mariana Arc, where 91.4% of the data are located in the gray band (Fig. 10j). As shown in Fig. 11, the percentages of α < 1 calculated from KL07 parameterization are close in the four depth ranges, amounting to approximately 80%.

Furthermore, ε and κ_v calculated by the dL20 parameterization show good consistency with the corresponding observations in the upper and middle layers. This is because 77.1% of the errors are less than one order of magnitude in the four generation sites (Fig. 10m), and the proportion decreases gradually with depth (Fig. 11). The results from the dL20 parameterization are 1 – 2 orders of magnitude higher than the observations at depths below 3 km, where the percentage of α < 1 only accounts for 33%. In addition, the percentages of α < 1 for the dL20 parameterization are 83.2%, 78.4%, and 75.8% in the Ryukyu Islands, Ogasawara-Mariana Arc, and Hawaii Islands, respectively. However, similar to the KL07 parameterization, the dL20 parameterization shows low performance in the Luzon Strait. Overall, in the generation sites, the dissipation rates calculated by the KL07 parameterization show the highest level of consistency with the observations among the three parameterizations. The model results in East Timor from the work of Cuypers et al. (2017) indicated that the KL07 parameterization likely underestimates the near-bottom dissipation of high-mode ITs. However, evaluating the performance of the three parameterizations entirely near the seafloor is difficult due to the lack of observations.

For the far-field regions, the dissipation rates estimated by the dL20 parameterization agree well with the observations. In particular, 87.5% and 89.9% of the errors for the Philippine Sea (Fig. 10l) and Kuroshio Extension area (Fig. 10r), respectively, are lower than one order of magnitude. At the sea bottom of the two regions, the energy dissipation is slightly higher than the observed values (Fig. 10f). Such a result may be attributed to the dissipation of low-mode tides by scattering on abyssal hills in dL20 parameterization. In addition, ε and κ_v are less intense than the observational results in the Kuroshio Extension region at depths of 0 – 1.5 km (Figs. 9f and 9l). The ITs also show no significant effects on this area, and mixing therein may be driven by other processes, such as mesoscale eddies and near-inertial internal waves (Li and Xu, 2014; Whalen et al., 2015).

4 Summary and Discussion

In this study, the IT-driven dissipation and diapycnal diffusion in the northern Pacific Ocean are estimated by the LSJ02, KL07, and dL20 parameterizations. The respective performances of the three parameterizations in the different regions are evaluated by comparing the calculation results with the observations. In particular, the calculation results from the three parameterizations show significant differences in terms of magnitude and vertical structure. The main conclusions are listed as follows.

Both the LSJ02 and KL07 parameterizations ignore the radiation of IT energy and only estimate the local dissipation and IT-driven mixing. However, both parameterizations focus on different depths, and the dissipation mechanisms for ITs are obviously different. In particular, the LSJ02 parameterization shows a bottom-intensified feature and a bottom dissipation that attenuates upward from the seafloor. This parameterization is only applicable to submarine ridges. Moreover, only 43% of the calculation results across the four generation sites have errors of less than one order of magnitude. The ITs can dissipate through direct shear instability (St. Laurent and Garrett, 2002) and hydraulic jumps (Legg and Huijts, 2006; Nikurashin and Legg, 2011) near the rough bottom topography. Nevertheless, the KL07 parameterization focuses on the upper and interior ocean, while the strongest dissipation shifts to the thermocline. Intense dissipation also occurs in the thermocline and decreases downward toward the sea bottom. Furthermore, the dissipation rate across all generation sites shows the highest level of consistency with the observations, as 80.7% of the results have errors of less than one order of magnitude. Additionally, the local dissipation of high-mode ITs can be enhanced by the nonlinear wave-wave interactions (MacKinnon et al., 2013; Yi et al., 2017).

The dL20 scheme considers the dissipation of low-mode ITs in the far-field regions through wave-wave and wave-topography interactions, thereby significantly enhancing the accuracy of remote tidal dissipation. For the dL20 parameterization, strong dissipation occurs in the upper and bottom layers, whereas the dissipation is relatively weak in the middle layer near the source regions. The errors lower than one order of magnitude account for 77.1% of the generation sites. Hence, the estimation results are in great agreement with the observations in the upper and middle layers, while the dissipation rate is significantly stronger than the observations at depths below 3 km. Moreover, 88.7% of the results are in accordance with the observations within one order of magnitude in the far-field regions.

In this study, the performance of the three tidal mixing parameterizations is preliminarily evaluated. However, the performance of the parameterizations near the bottom is not evaluated due to the limited observations in the abyssal ocean. Thus, other near-bottom observational studies should be performed for the further validation and modification of parameterizations. In addition, the vertical structure functions of the KL07 parameterization and the dissipation caused by the wave-wave interactions of the dL20 parameterization both depend on the buoyancy frequency, whose variation in the upper ocean is not only influenced by ITs but also by wind-driven near-inertial waves, mesoscale eddies, and so on. However, the physical processes are complicated, and further work is needed to understand the relevant mechanisms.

Acknowledgements

This study is supported by the National Key Research and Development Program of China (No. 2017YFA0604103), the National Natural Science Foundation of China (No. 41876015), the Fundamental Research Funds for the Central Universities (No. 202061001), and the Open Innovative Fund of Marine Environment Guarantee (No. HHB003).