Impact of Enhanced Wave-Induced Mixing on the Ocean Upper Mixed Layer during Typhoon Nepartak in a Regional Model of the Northwest Paciﬁc Ocean

: To investigate the e ﬀ ect of wave-induced mixing on the upper ocean structure, especially under typhoon conditions, an ocean-wave coupled model is used in this study. Two physical processes, wave-induced turbulence mixing and wave transport ﬂux residue, are introduced. We select tropical cyclone (TC) Nepartak in the Northwest Paciﬁc ocean as a TC example. The results show that during the TC period, the wave-induced turbulence mixing e ﬀ ectively increases the cooling area and cooling amplitude of the sea surface temperature (SST). The wave transport ﬂux residue plays a positive role in reproducing the distribution of the SST cooling area. From the intercomparisons among experiments, it is also found that the wave-induced turbulence mixing has an important e ﬀ ect on the formation of mixed layer depth (MLD). The simulated maximum MLD is increased to 54 m and is only 1 m less than the observed value. The wave transport ﬂux residue shows a dominant role in the mixed layer temperature (MLT) changing. The mean error of the MLT is reduced by 0.19 ◦ C compared with the control experiment without wave mixing e ﬀ ects. The study shows that the e ﬀ ect of wave mixing should be included in the upper ocean structure modeling.


Introduction
Vertical mixing of the ocean is known to be one of the most important processes in regulating the SST and mixed layer structure, which are two key ocean factors that control the exchange of heat and momentum between the ocean and atmosphere. For a model to be able to simulate realistic ocean dynamics, it must include an accurate characterization of the vertical mixing process. However, due to the small-scale turbulence processes involved, vertical mixing usually cannot be explicitly expressed in general ocean circulation models and, so far at least, its implementation must be parameterized.
The vertical mixing in many three-dimensional numerical ocean circulation models is often based on two-equation turbulence closure schemes [1][2][3][4][5][6][7][8]. The turbulent kinetic energy (TKE) equation has a firm physical basis, while the equation for length scale, no matter which form is used, is merely patterned after the TKE equation and consequently, its physical basis is tenuous. Unlike the TKE equation, all terms, like production, destruction, and diffusion term in the length scale equation, Remote Sens. 2020, 12 have to be modeled, which also is a major difficulty in two-equation turbulence models of this type [9]. A common problem of such schemes is an overestimation of SST and underestimation of mixed layer depth (MLD), especially during the summer [1,[10][11][12][13][14][15]. For instance, Martin [15] showed that summertime temperatures at Ocean Weather Station Papa exceeded observations by about 2 • C relative to an annual range of 8 • C. Recently, more attention has been paid to the role of the waves which believed to be the primary reason for the insufficient vertical mixing in the upper ocean. By careful tuning, many studies showed some improvement in the final results [9,13,[16][17][18]. However, more numerical modeling investigation is needed, especially under different real ocean conditions. At the same time, some studies are concentrated on breaking waves which can be applied to ocean model by using surface boundary conditions. However, further studies have shown that the wave breaking generated turbulence is mostly confined to the surface layer with their energy dissipated locally at a depth of the order of the wave amplitude [19,20]. Yet, all these efforts gave us a strong hint that wave could not be neglected in the study of ocean dynamics. Typhoons (or tropical cyclones (TCs)) are probably among the most extreme sea-air interactions. These extremes can drastically change the marine environment and seriously affect human activities. Thus, the ability to simulate changes in the marine environment during TCs is vital. The large waves caused by TCs represent a main feature of the sea surface and exert an important influence on the vertical mixing of the upper ocean. Observations and model experiments show that sea surface cooling caused by the TC can exceed 3 • C over the hundreds of kilometers of the TC center, and cooling of 7-11 • C may occur and could last for 1-3 weeks [21][22][23]. Moreover, a strong TC entrains cold water into the mixed layer (ML) from below, resulting in the cooling of the ML and deepening of the MLD [24]. Wang and Qiao [25] introduced the mixing effect of wave-induced turbulence in the ocean current model, and the results show that ocean waves have a significant influence on the decline of SST and deepening of the MLD during a TC. However, many studies have also shown that even with this added mechanism, the strength of thermocline is still weak [26][27][28][29]. There are many other deficiencies in the simulation results under TC conditions, such as large differences in the range, intensity and spatial structure of sea surface cooling, offset of the strong cooling zone, and abnormal structure of the upper mixed layer and so on. The lack of skill of simulation ability may be due to the insufficient consideration of the physical processes. Based on field observations, Toffoli et al. [30] found that variations of a normally shallow mixed layer depth are observed within a relatively short timescale of approximately 10 h after the influence of TC. The quick response of the mixed layer seems to sustain the role of wave-induced turbulence in the upper ocean mixing. Moreover, using coupled hurricane-ocean-wave modeling, Aijaz et al. [31] found that the effect of wave breaking on SST cooling is minor except the center of the TC track. In addition, the model skill is improved by enhancing the turbulence due to nonbreaking waves.
In this study, a 3D wave-current coupled model consisting of the Princeton Ocean Model (POM) and the Marine Science and Numerical Modeling (MASNUM) wave model is used to investigate to what degree the wave-induced mixing can affect the upper ocean superstructure. Except for the wave-induced turbulence mixing mechanism, another wave-induced mixing mechanism called the wave transport flux residue is also introduced into the ocean circulation model. Preliminary research has shown that the latter can affect the structure of the upper ocean as a driving source. The Northwest Pacific ocean is used as the research area, and TC Nepartak is used as a TC example, and the roles of the two mechanisms in ocean mixing are studied through comparative experiments, especially in the case of TC. A series of comparisons are conducted between the experimental results and the available measured data. This paper is structured as follows: The two physical processes are described in detail in Section 2; the descriptions of material and method used in this study are included in Section 3; the model results are presented in Section 4; the relevant discussions are given in Section 5 and the conclusions are summarized in Section 6.

Wave-Induced Turbulence Mixing
Surface wave breaking and nonbreaking waves can enhance diffusivity, particularly vertical diffusivity, in the ocean. The turbulence generated by surface wave breaking decays rapidly with distance from the sea surface and is therefore considered as a minor contributor to the overall ocean mixing [20,32]. Studies [33,34] have shown that nonbreaking wave-induced turbulence, which is distributed vertically through the water column at a scale comparable to that of the wave length [35], is a predominant mechanism for turbulence generation and enhancement and can significantly affect the physics of ocean mixing.
Xia et al. [28], Qiao et al. [36], Dai et al. [37] and others have shown that the ability of the model to simulate ocean mixed layers and SST can be significantly improved via a rational consideration of wave-induced turbulence mixing (hereafter WITM) during ocean modeling. Qiao et al. [36], Qiao et al. [38] expressed the vertical mixing, Bv, caused by nonbreaking waves as follow: where the parameter α = (1) and is set to 1 following [36,38]; E( → k ) represents the wave number spectrum, which can be calculated from the numerical wave model; ω is the angular frequency of the wave; k is the wave number, and z is the water depth. The nonbreaking WITM is directly incorporated into the ocean model as follows: where K m and K h are the vertical viscosity and diffusion coefficients respectively. K mc and K hc are calculated by the Mellor-Yamada turbulence closure scheme employed in the ocean model.

Wave Transport Flux Residue
Based on the differences in physical characteristics and time-space scales, Yuan et al. [39] divided the ocean motions into four types: turbulence, wave-like motion, eddy-like motion, and circulation. At the same time, a group of Reynolds averages, including first-fold marked by SS , second-fold marked by SM , and third-fold marked by MM , is defined to decompose the ocean motions. For example, by applying the three Reynolds average operations to the velocity, we can decompose it into four parts as shown in Equations (3)-(6), where the U is the total ocean motion, the U is the sums of motions except ocean turbulence, theÛ is the sums of the circulation, eddy-like motion, and the R.H.S. of the Equation (6) represent the circulation, eddy-like motion, wave-like motion, and turbulence, respectively. Moreover, with the first two Reynolds average operations and the ocean motions governing equations, we can get the governing equations for the sum of ocean eddy-like motion and circulation: Remote Sens. 2020, 12, 2808 4 of 26 and the governing equations for the wave-like ocean motion sub-system: Here x i (i = 1, 2, 3) denote the eastward, northward, and vertical of the coordinate, respectively. u SMi (i = 1, 2, 3), T SM , s SM , p SM , ρ SM are the velocity, temperature, salinity, pressure, and density for the wave-like ocean motion, andÛ i (i = 1, 2, 3),T,ŝ,p,ρ for the sum of ocean eddy-like motion and circulation, which also act as the background field components in the Equations (12)- (16). ν o , κ o , D o are the coefficients of molecular viscosity, thermal conductivity, and the diffusion coefficient of salinity. Q and Q SM are the equivalent temperature source terms. It consists of four components which are net downward shortwave, net downward longwave, sensible, and latent fluxes, and is calculated using the Coupled Ocean-Atmosphere Response Experiment (COARE) bulk formula [40]. The latter three components are assumed to occur at the ocean surface through the surface boundary conditions. The downward shortwave flux is specified as an absorption profile which is interpreted by Paulson and Simpson [41]. Each sub-system has two effects in the ocean dynamic cascade system. The governing equation group of each sub-system gives the expression of these two effects. The first is the advection transport and shear instability generations of larger scale motions through which modulate smaller scale motions (the second and third terms on L.H.S. of Equations (13)- (15)). The large-scale variables in the sub-system are solved in its superior system. The second is the ocean mixing induced by smaller scale motions in the form of transport flux residues through which modulate larger scale motions (the last term on R.H.S. of Equation (8), and the second term on R.H.S. of Equations (9) and (10). The small-scale variables in the sub-system are solved in its subordinate system. The transport flux residue expresses the amount left after applying the Reynold's average on the flux of small scale Remote Sens. 2020, 12, 2808 5 of 26 motions [39]. The effects of wave-like ocean motions on circulation, we call wave transport flux residue (hereafter WTFR), include: Furthermore, considering the background field shear, we decomposed the wave-like ocean motions into two components, i.e., where u SMi (i = 1, 2, 3), T SM , s SM are the velocity, temperature, and salinity of the dominant gravitational wave components, andû SMi (i = 1, 2, 3),T SM ,ŝ SM are the perturbation components generated by background field shear. With Equations (17)- (19), the wave transport flux residue can be expanded as: On the R.H.S. of Equation (20), the first term is the radiation stress and was widely used. Here, the ocean wave field is assumed to be locally uniform, i.e., T SM = 0, s SM = 0, so that the first and third terms on the R.H.S. of Equations (21) and (22) (20) is not investigated in this study. The second terms on the R.H.S. of Equations (20)- (22) are the main components of WTFR. With the continuity equation for u SMi (i = 1, 2, 3), they can be expressed as: Since ocean waves are locally uniform, the horizontal changes of statistic parameters for ocean waves within the mixing length are near zero [36]. Therefore, only the vertical terms in WTFR are relevant, and the final expressions are given in Equations (26)- (28). The more detailed derivation processes are summarized in Yang et al. [42].
− u SM3 Remote Sens. 2020, 12, 2808 where H is the local bottom topography. The transport flux residue represents the effect of small-scale motion to large scale motion and is dependent linearly on the background field shear, so this value may sometimes become negative. In positive value condition, the WTRF extract momentum, heat, or salinity from the atmosphere to the ocean. The negative value corresponds to upward flux transfer, i.e., from ocean to atmosphere. In other words, the WTFR may uniform or stratify the ocean horizontally or vertically.
Here we define the vertical WTFR terms in momentum, temperature, and salinity equations: The effects of the WTFR are introduced into the ocean model as follows: The above four equations represent the momentum, temperature, and salinity equations in the ocean model. D is the total water depth, D = H + η where η is the surface elevation, σ is the vertical direction in the sigma coordinate, f is the Coriolis parameter, p is the pressure, and R is the short wave radiation flux. This scheme is applied in the real ocean situation simulation for the first time, especially under TC condition, we consider only the effect of the WTFR in the temperature equation, which is also beneficial to the analysis of the experimental results.

Observation Data
To more objectively analyze the impact of each physical process, we compared the experimental results with the observed data. This subsection lists the observed data we have used and provides a general introduction.

Jason-2
In this study, we used the well-calibrated Geographical Data Record (GDR, https://www.aviso. altimetry.fr) significant wave height (SWH) data measured by the Jason-2 satellite altimeter to test the simulation performance of the wave model. The satellite was launched by NASA and the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) and is used mainly for measuring sea level and SWH with a high degree of accuracy. These data have been widely used in marine and climate research [43,44]. The prelaunch SWH error budget specification for Jason-2 is 10% Remote Sens. 2020, 12, 2808 7 of 26 or 0.4 m [45]. The later validation studies show that Jason-2 measured SWH has a 2-6% difference compared to buoys data [46,47].

MW_IR
To test the performance of the ocean model, we used the optimally interpolated daily SST data product produced by the US Remote Sensing System (MW_IR, http://www.remss.com). These data combine the through-cloud capabilities of microwave data (MW) with the high spatial resolution of infrared (IR) data and have a resolution of nearly 9 km. The accuracy of the data is guaranteed through quality control and optimal interpolation schemes [48]. The measured SSTs are shown to have a mean bias of −0.07 • C and a standard deviation of 0.57 • C when compared to the buoy data [49].

Argo
The Array for Real-time Geostrophic Oceanography (Argo) floats measure the temperature and salinity of oceans (0-2000 m) globally, even under extreme atmospheric and oceanic conditions, such as a TC passage. Therefore, Argo data can well reflect the upper ocean structure and variability under the impact of a TC [50]. During the passage of TC Nepartak, an Argo float numbered WMO 2902666 was found, and it measured twelve profiles. Furthermore, another 69 profiles near the TC track were also used to study the responses of the upper ocean during TC ( Figure 1b). The data were downloaded from the Coriolis Global Data Center in France (ftp://ftp.ifremer.fr).

Buoy
Yang et al. [51] described some very useful observations made a new kind of deep-sea data buoy that is able to survive encounters with TC Nepartak. The observed data are therefore quite valuable and provide a good picture of the rapid changes in the upper ocean temperature field that occur during the direct forcing phase of the TC. In this study, we digitized part of the image data in their paper and used it to compare with the experimental results. The locations of the two buoys, NTU1 and NTU2, are marked in Figure 1b with yellow and blue pentacle stars, respectively.

Ocean Model
The Princeton Ocean Model (POM; [52]) was used to simulate the circulation. POM is a three-dimensional model using Arakawa-C grid staggering, with sigma coordinates in the vertical direction to handle complex topographies and shallow regions. The turbulence closure scheme employed in the POM is the Mellor-Yamada scheme, which is based on the level-2.5 turbulent kinetic energy equations [2].
The model domain covers the geographical area (0 • ∼ 50 • N, 99 • ∼ 160 • E, see Figure 1a), has a horizontal resolution of (1/30) • × (1/30) • that corresponds to 3.7 km in latitude and varies between 3.7 and 2.5 km in longitude. A total of 73 vertical sigma layers are included, and the vertical grid spacing had a higher resolution in the surface and bottom boundary layers and avoided artificial diffusion in large water depths. The highest vertical resolution is D × 3.6 × 10 −4 meters and the lowest resolution is D × 5.5 × 10 −2 meters, where the D is total water depth. The topography of the model is interpolated based on the General Bathymetric Chart of the Oceans (GEBCO_08), and its resolution is (1/120) • × (1/120) • . The minimum depth was set to 5 m, and the maximum depth was set to 5500 m. The mode splitting technique [53,54] is used in this model for the computer economy. The external model is 2-D with a time step of 3 s, and the internal model is 3-D uses a longer time step of 90 s. The boundary conditions are based on the (1/12) • × (1/12) • resolution HYCOM analysis data (http://hycom.org/dataserver/). The boundary forcing variables include sea surface level, velocity, salinity, and temperature. The initial condition was obtained from the model result of Zhao et al. [55] with the initial time 1 January 2013. maximum depth was set to 5500 m. The mode splitting technique [53,54] is used in this model for the computer economy. The external model is 2-D with a time step of 3 s, and the internal model is 3-D uses a longer time step of 90 s. The boundary conditions are based on the (1/12) ∘ × (1/12) ∘ resolution HYCOM analysis data (http://hycom.org/dataserver/). The boundary forcing variables include sea surface level, velocity, salinity, and temperature. The initial condition was obtained from the model result of Zhao et al. [55] with the initial time 1 January 2013. The model is forced by both wind stress and surface heat fluxes. According to the different types of forcing data, the simulation procedure includes three parts. In the first (from the year 2013 to 2015), Global Forecast System (GFS) wind data provided by the National Centers for Environmental Prediction (NCEP) are used, and wind stresses were computed using a bulk formula. The data set has spatial and temporal resolutions of 0.25 ∘ × 0.25 ∘ and 3 hours, respectively. The drag coefficient formulation which can be used for low-to-moderate winds and high wind speed is adopted [56]. The climatological monthly mean data from the Comprehensive Ocean-Atmosphere Data Set (COADS; [57]) are used as the surface heat fluxes. The model is integrated for 3 years reaching a stable state. For the second part (the year 2016), the climatological monthly mean surface heat fluxes are replaced by real time data. This part is used to assess the performance of both the wave and ocean model, and also to evaluate the impact of the two schemes The model is forced by both wind stress and surface heat fluxes. According to the different types of forcing data, the simulation procedure includes three parts. In the first (from the year 2013 to 2015), Global Forecast System (GFS) wind data provided by the National Centers for Environmental Prediction (NCEP) are used, and wind stresses were computed using a bulk formula. The data set has spatial and temporal resolutions of 0.25 • × 0.25 • and 3 h, respectively. The drag coefficient formulation which can be used for low-to-moderate winds and high wind speed is adopted [56]. The climatological monthly mean data from the Comprehensive Ocean-Atmosphere Data Set (COADS; [57]) are used as the surface heat fluxes. The model is integrated for 3 years reaching a stable state. For the second part (the year 2016), the climatological monthly mean surface heat fluxes are replaced by real time data. This part is used to assess the performance of both the wave and ocean model, and also to evaluate the impact of the two schemes on the ocean superstructure under general condition. The GFS meteorological variables are used to calculate the surface heat fluxes using the Coupled Ocean-Atmosphere Response Experiment (COARE) bulk formula [40]. For the last part (during TC Nepartak), to better and emphatic analyze the influence of the two schemes during the TC period, a wind field and wind pressure model is used to construct the wind field and pressure field during the TC Nepartak which will be discussed in Section 3.2.3.

Wave Model
The Marine Science and Numerical Modeling (MASNUM) wave model is a third-generation wave model that was originally designed by Yuan et al. [58,59]. Subsequently, Yang et al. [60] updated the model to increase its suitability for the global scale. The model solves the energy balance equation through the wave number spectrum space. The complicated characteristics inlaid scheme is used for wave energy propagation. Source functions including wind input, dissipation, wave-wave interaction, and wave-current interaction are introduced in this wave model. The wind input source function is based on Phillips' resonance and Miles' shear flow instability mechanisms [61,62]. The wave spectrum of the model is discretized using 24 directions and 25 frequency bins with a minimum frequency of 0.0412Hz. The frequency increases logarithmically, and the increment factor is 1.21. The in situ measurement based formulation of Wu [63], which is applicable even in hurricanes, is adopted as drag parameterization in MASNUM. The model domain and resolution of the wave model are consistent with that of the POM. The boundary conditions of the wave model are from the JONSWAP spectrum [64]. The MASNUM wave model is integrated for 4 years independently, from the year 2013 to 2016, with the same period as the ocean model. The surface wind stress data used by the wave model are also the same as those used by the ocean model.
In the first two parts, the WITM (Bv, see Equation (1)) and the integral part of WTFR (Bsmt, see Equation (27)) are generated as daily output and introduced into the ocean model subsequently. For the third part, the data output and transferring frequency are increased to 3 h equal to the surface forcing data temporal resolutions.

Wind Field and Wind Pressure Model
During TC Nepartak, the wind and pressure fields are generated by a typhoon wind model. The input conditions for the model are from the Joint Typhoon Warning Center (JTWC, https: //metoc.ndbc.noaa.gov/ja/web/guest/jtwc/best_tracks), including the features of the TC, such as warning time, estimation of position, maximum sustained surface wind speed, minimum central pressure, propagation speed and direction, and radius of maximum wind speed. The wind field and wind pressure model gives the pressure P at any radius (r) from the typhoon center as follows [65,66]: where P ∞ is the ambient pressure; P 0 is the central pressure; P(r) is the surface pressure at a distance r from the typhoon center and R is the radius of maximum wind speed. Using Equation (34), the gradient wind speed can be calculated via gradient wind theory: where ∆P is the difference between the ambient pressure and the center pressure; f is the Coriolis parameter; and ρ a is the air density (assumed to be constant at 1.2 kg m 3 ). The wind field derived from the Equation (35) is symmetrical to the typhoon center. Considering the asymmetry of the real wind field, two other factors, the angle of inflow and typhoon forward motion, should be included. In this study, the angle of inflow is set to be constant (θ = 20 • ), and the formula proposed by Ueno [67] is employed to calculate the forward wind field: where V x and V y are the components of the forward speed of the typhoon speed, with V x being positive eastward and V y being positive northward; and → i and → j are the unit vectors in the x and y directions, respectively. The model's wind field comes to: where (x 0 , y 0 ) is the typhoon center position and (x, y) is the calculated position; c 1 and c 2 are the correction coefficients, set to be 1.0 and 0.8, respectively.
The constructed wind field for typhoon is combined with the background wind field as follow [68]: (38) where → V C is the combined wind field and → V Q is the background wind field, and n is a constant, set to be 10. The GFS wind field is employed as the background wind field.

Overview of the Typhoon and Experimental Design
For a better interpretation of the results to be reported soon, it is necessary to provide a short description of the TC that we have considered for the present study. In the next subsection, we describe the experiments we have carried out to test the sensitivity of the model results to the two physical processes described in Section 2.

Overview of the Typhoon
At 1600 UTC 2 July 2016, Nepartak was formed on the surface of the Northwest Pacific Ocean; it then intensified into a strong tropical storm at 1200 UTC 3 July 2016, strengthened at 1800 UTC 4 July 2016 to a typhoon, and then further strengthened at 1800 UTC 4 July 2016 to a super typhoon. At 1200 UTC 6 July 2016, Nepartak increased to its peak intensity of 79.9 m/s and 907 hPa; at 0550 UTC 8 July 2016, it landed at Taitung County, Taiwan Province; at 1345 UTC 9 July 2016, it landed again at Shishi City, Fujian Province; and at 0600 UTC 9 July 2016, it then weakened to a tropical depression ( Figure 1b).
TC Nepartak had a strong intensity and reached a super typhoon level, which had a drastic impact on the upper ocean structure. During the TC period and a long time thereafter, the TC transit area was not affected by other TCs or strong weather events, which was more conducive to the analysis of the experimental results.

Experimental Design
To study the effects of WITM and WTFR on ocean vertical mixing, especially during the TC stage, we designed three experiments, denoted as EX1, EX2, and EX3. EX1 is the control experiment, which is the original ocean current model; EX2 introduces WITM into the ocean current model; EX3 introduces both two physical processes (Table 1).

Measures of Quality for Models Validation
For quantitative comparison, five statistical parameters, including mean error (ME), root mean square error (RMSE), correlation coefficient (r), Nash-Sutcliffe efficiency coefficient (E), and determination coefficient (R 2 ) were used and briefly summarized as below: where M represents the model-generated data, O represents the observational data, the bar over M and O denotes their mean values,M denotes the fitted value, and N is the number of collocated data points. The ME is a statistical quantity that reveals any systematic error between the model and buoy data. The RMSE indicates the degree of discretization of data pairs; lower RMSE values indicate a smaller degree of data discretization between the model and observations. The r value indicates how well the model data fit the observations. E is generally used to verify the model simulation results. The value of E is minus infinity to 1, and E is close to 1, indicating the good quality of the model and high reliability of the model. R 2 is used to judge the quality of data fitting. Its value range is between 0 and 1, and the larger value shows a higher fitting quality.

Results
In this section, first, the comparative results for the entire year of 2016 are used to assess the performance of both the wave and ocean model, the latter also indicates the contribution of the two physical processes to the improvement of the ocean model simulation capability under general condition. Second, a more detailed comparison has been performed using the simulation results of the three experiments and the observation data recorded during the TC Nepartak. The model data are daily averaged during TC period to correspond with the observations which are available as a single data set for each day.

Evaluation Under General Condition: The year 2016
First, we use the SWH data of the Jason-2 satellite to check the simulation results of the wave model. The scatter distribution statistics of the SWH for the year 2016 in the simulated area are shown in Figure 2a. Most of the data points lie near the diagonal line, indicating that the simulation results are consistent with the satellite observations. The error probability distribution (Figure 2b) shows that as the error increases, the probability distribution decreases rapidly, and most of the errors are within an acceptable range of 0.4 m. The r, E, and R 2 value are about 0.95, 0.77, and 0.79, respectively (see Table 2). The statistic comparisons between simulated SWH and the observations show MASNUM's good simulation capabilities which provide a good basis for our study of the mixing effects of waves. results are consistent with the satellite observations. The error probability distribution (Figure 2b) shows that as the error increases, the probability distribution decreases rapidly, and most of the errors are within an acceptable range of 0.4 m. The r, E, and R 2 value are about 0.95, 0.77, and 0.79, respectively (see Table 2). The statistic comparisons between simulated SWH and the observations show MASNUM's good simulation capabilities which provide a good basis for our study of the mixing effects of waves.  The box plot of the spatially averaged SST error in the simulated area is given in Figure 3.   The box plot of the spatially averaged SST error in the simulated area is given in Figure 3.   Figure 3 shows that the average error of the SST simulated by the three experiments is within 1 ℃, and the median error of each month is within 0.5 ℃, which is relatively small. Except for January, the mean SST of each experiment is higher than the observations, with a positive error. However, compared with EX1, the error range of EX2 and EX3 is closer to the zero line, especially in summer (June, July, and August) and early autumn (September). After introducing the mixing effect of the waves, the two experiments (EX2, EX3) show similar simulation results, and the difference between the model simulated SST and the observation is decreased by more than 10% ( Table 2). The results of statistical analysis show that the E and R 2 value is increased from 0.76 to 0.84 and 0.78 to 0.89, which indicate an improvement of the ability of the model to simulate SST.
The above analyses show that both wave and ocean models perform well and the gradual introduction of two physical processes (in EX2 and EX3) has played a positive role in improving the  Figure 3 shows that the average error of the SST simulated by the three experiments is within 1 • C, and the median error of each month is within 0.5 • C, which is relatively small. Except for January, the mean SST of each experiment is higher than the observations, with a positive error. However, compared with EX1, the error range of EX2 and EX3 is closer to the zero line, especially in summer (June, July, and August) and early autumn (September). After introducing the mixing effect of the waves, the two experiments (EX2, EX3) show similar simulation results, and the difference between the model simulated SST and the observation is decreased by more than 10% ( Table 2). The results of statistical analysis show that the E and R 2 value is increased from 0.76 to 0.84 and 0.78 to 0.89, which indicate an improvement of the ability of the model to simulate SST.
The above analyses show that both wave and ocean models perform well and the gradual introduction of two physical processes (in EX2 and EX3) has played a positive role in improving the ocean model simulation capability as a whole. The error of SST has been improved to varying degrees, most notably in summer and early autumn. This is further discussed in Section 5.

Evaluation Under TC Condition: Nepartak
The simulated wave field is evaluated during the TC period using the orbit data of the Jason-2 satellite (Figure 4) ℃, and the median error of each month is within 0.5 ℃, which is relatively small. Except for January, the mean SST of each experiment is higher than the observations, with a positive error. However, compared with EX1, the error range of EX2 and EX3 is closer to the zero line, especially in summer (June, July, and August) and early autumn (September). After introducing the mixing effect of the waves, the two experiments (EX2, EX3) show similar simulation results, and the difference between the model simulated SST and the observation is decreased by more than 10% ( Table 2). The results of statistical analysis show that the E and R 2 value is increased from 0.76 to 0.84 and 0.78 to 0.89, which indicate an improvement of the ability of the model to simulate SST.
The above analyses show that both wave and ocean models perform well and the gradual introduction of two physical processes (in EX2 and EX3) has played a positive role in improving the ocean model simulation capability as a whole. The error of SST has been improved to varying degrees, most notably in summer and early autumn. This is further discussed in Section 5.

Evaluation Under TC Condition: Nepartak
The simulated wave field is evaluated during the TC period using the orbit data of the Jason-2 satellite (Figure 4)  The response of the SST on 11 July 2016 is presented in Figure 5. For this study, similar to Zhang, et al. [69], the response of SST was defined as the SST change relative to its value at pre-typhoon stage (1 July 2016). The model data are interpolated to the observation points before plotting results. All experimental and observation results show a strong SST response to the TC. Comparing with the EX1, the SST cooling area in EX2 and EX3 is significantly increased, which is more consistent with the observations. A noteworthy phenomenon is that the maximum cooling area in observation as well as in EX3 is located to the R.H.S. of the TC path, which is well demonstrated in field observation [70,71] and multiple numerical experiments [25,72,73]. On the right side of TC track, the resonance between wind and current rotations and also between the wind and the TC translation result in larger waves which increase the turbulence and therefore enhance the SST cooling ultimately. However, in EX2, the maximum cooling area shows a left side shift tendency. After introducing the effect of WTFR (Figure 5c), the distribution of the maximum cooling area is considerably improved. demonstrated in field observation [70,71] and multiple numerical experiments [25,72,73]. On the right side of TC track, the resonance between wind and current rotations and also between the wind and the TC translation result in larger waves which increase the turbulence and therefore enhance the SST cooling ultimately. However, in EX2, the maximum cooling area shows a left side shift tendency. After introducing the effect of WTFR (Figure 5c), the distribution of the maximum cooling area is considerably improved. In general, the above comparisons indicate that the ability of the model to simulate SST during the TC gradually increased via the introduction of two different wave effects. The WITM remarkably improves the problem associated with an insufficient SST cooling area and cooling amplitude, although its impact on the distribution of the cooling area is relatively small. The phenomenon, i.e., left shift of the maximum cooling area, can be well modified after the introduction of WTFR.
To clearly show the variations of the upper ocean during TC and evaluate the role of wave mixing in these processes, observations from the two buoys, NTU1 and NTU2, are used [51]. Figure  6a and 7a show the comparison of time series wind and pressure data between observed and model used at NTU1 and NTU2, respectively. The constructed wind and pressure data show great consistency with the observation, although the wind data are flatter to the observed. The upper 25 m-averaged temperature (T25m) obtained from the observations and the simulations are compared (Figures 6b and 7b). The model output results in every 3 hours, and the time axis indicates the output time. The comparison shows that with the introduction of wave-induced mixing, there is a In general, the above comparisons indicate that the ability of the model to simulate SST during the TC gradually increased via the introduction of two different wave effects. The WITM remarkably improves the problem associated with an insufficient SST cooling area and cooling amplitude, although its impact on the distribution of the cooling area is relatively small. The phenomenon, i.e., left shift of the maximum cooling area, can be well modified after the introduction of WTFR.
To clearly show the variations of the upper ocean during TC and evaluate the role of wave mixing in these processes, observations from the two buoys, NTU1 and NTU2, are used [51]. Figures 6a and  7a show the comparison of time series wind and pressure data between observed and model used at NTU1 and NTU2, respectively. The constructed wind and pressure data show great consistency with the observation, although the wind data are flatter to the observed. The upper 25 m-averaged temperature (T25m) obtained from the observations and the simulations are compared (Figures 6b and  7b). The model output results in every 3 h, and the time axis indicates the output time. The comparison shows that with the introduction of wave-induced mixing, there is a marked drop in temperature. Compared with EX1, the results of EX2 have been improved to some extent, but they are still not sufficient; for example, the RSME is decreased from 1.92 • C to 1.43 • C at NUT1 and from 1.67 • C to 0.99 • C at NTU2 (Table 3). EX3's results are more consistent with the observations, and the RSME becomes 0.68 • C and 0.45 • C at NTU1 and NTU2, respectively. It is worth noting that the EX1 and EX2 give a negative E value at two buoy locations and the model simulation results are unsatisfactory. In EX3, the model performance is greatly improved, and the E value is 0.61 and 0.74, respectively. Figures 6c and 7c are the temperature profiles at the corresponding moment ticked in Figures 6b and 7b which show that the effect of wave-induced mixing is mainly concentrated on the upper 150 m. It can be seen that among the three experiments, the results of the EX3 are more in line with the observations. worth noting that the EX1 and EX2 give a negative E value at two buoy locations and the model simulation results are unsatisfactory. In EX3, the model performance is greatly improved ,and the E value is 0.61 and 0.74, respectively. Figures 6c and 7c are the temperature profiles at the corresponding moment ticked in Figures 6b and 7b which show that the effect of wave-induced mixing is mainly concentrated on the upper 150 meters. It can be seen that among the three experiments, the results of the EX3 are more in line with the observations.      Here, we compared the surface heat flux between the model calculated and observed at the buoy station NTU1 and NTU2 (Figure 8). Before and after the effect of TC, the difference of the calculated surface heat flux between the three experiments is relatively small, and the data are in good agreement with the observations. During the TC period, the calculated data are more accurate in the experiments which included the wave effect (EX2 and EX3).  Here, we compared the surface heat flux between the model calculated and observed at the buoy station NTU1 and NTU2 (Figure 8). Before and after the effect of TC, the difference of the calculated surface heat flux between the three experiments is relatively small, and the data are in good agreement with the observations. During the TC period, the calculated data are more accurate in the experiments which included the wave effect (EX2 and EX3). On the other hand, the response of the upper ocean to the TC is also reflected in the changes of the MLD and MLT. To clearly show this structure and its changes as well as the role of wave mixing in these processes, we use Argo float data (Figure 1b) to test the simulation results of the three experiments. The MLD is defined as the depth at which the temperature is 0.5 • C lower than the first layer. The MLT is then defined as the average temperature in this depth range.
The observational results show that at the Argo float WMO 2902666 (the cross in Figure 1b of the WITM (EX2), the maximum MLD reached 54 m, which is very close to the observed value. The result of EX3 accorded basically to the result by EX2 with only subtle differences in the subsequent oscillation stage. The observational and all three experimental results show that the MLD return to its pre-storm depth approximately 7 days after the passage of the TC (14 July). By comparing three experiments results with the observation, it is found that the WITM has an important effect on the formation of MLD, especially at the deepening stage.  Table 1 for the specification of the different experiments. Time is labeled when observation data is available.
The observed MLT decreased rapidly from 29.9 ℃ (before the TC impact) to 27.3 ℃ ( Figure  9b), with the maximum cooling exceeding 2 ℃. The pre-typhoon MLT of EX1 is 31.1 ℃, which is 1.2 ℃ higher than the observed. Affected by the TC, the minimum MLT declined to 30 ℃, with a maximum decrease of 1.1 ℃. After considering the effect of WITM (EX2), the simulated MLT cooling amplitude increased to 1.3 ℃, decreased from 30.3 ℃ to 29 ℃. With WTFR effect (EX3), the performance of the MLT simulation is further improved with a maximum decrease of 2.1 ℃. In EX3, the MLT warmed up to the pre-storm status until 24 July which shows a good agreement with the observation. However, the warming trend of the MLT in EX1 is not obvious and gradually returned to the pre-storm level until the last day of July (purple cross).
Below, three experimental results are compared with another 69 Argo observations during TC Nepartak to further study the effects of the two physical processes on the upper ocean ( Figure 10). Shallow observed MLD (10~30 m) and high observed MLT (warmer than 30 ℃ ) which we considered as storm-free status. All three experimental results are relatively consistent with the observation at that status. Affected by the TC, as observed MLD deepened and MLT decreased,  Table 1 for the specification of the different experiments. Time is labeled when observation data is available.
The observed MLT decreased rapidly from 29.9 • C (before the TC impact) to 27.3 • C (Figure 9b), with the maximum cooling exceeding 2 • C. The pre-typhoon MLT of EX1 is 31.1 • C, which is 1.2 • C higher than the observed. Affected by the TC, the minimum MLT declined to 30 • C, with a maximum decrease of 1.1 • C. After considering the effect of WITM (EX2), the simulated MLT cooling amplitude increased to 1.3 • C, decreased from 30.3 • C to 29 • C. With WTFR effect (EX3), the performance of the MLT simulation is further improved with a maximum decrease of 2.1 • C. In EX3, the MLT warmed up to the pre-storm status until 24 July which shows a good agreement with the observation. However, the warming trend of the MLT in EX1 is not obvious and gradually returned to the pre-storm level until the last day of July (purple cross).
Below, three experimental results are compared with another 69 Argo observations during TC Nepartak to further study the effects of the two physical processes on the upper ocean ( Figure 10). Shallow observed MLD (10~30 m) and high observed MLT (warmer than 30 • C) which we considered as storm-free status. All three experimental results are relatively consistent with the observation at that status. Affected by the TC, as observed MLD deepened and MLT decreased, three experiments show different performance. The response of the original ocean model (EX1) to the TC is insufficient and this phenomenon is gradually improved with the effects of the wave-induced mixing (EX2, EX3). Table 4 gives the comparisons of the MLD and MLT between Argo float observations and three experimental results during the TC period. In terms of MLD, EX1 presents ME of −8.97 m. In EX2, ME has fallen by more than 68% and EX3 shows a similar result. For the MLT, in the process of introducing the WITM (EX2) and WTFR (EX3), the model results are gradually improved. Especially, in EX3, the ME is decreased by 40% from 0.3 • C to 0.12 • C. WITM can improve the simulated MLD result observably, the E value is increased from 0.32 (EX1) to 0.90 (EX2), while WTFR plays a crucial role in the MLT changing, the E value is increased from 0.49 (EX1) to 0.75 (EX3). This is consistent with the previous results at Argo float WMO 2902666 (Figure 9). presents ME of -8.97 m. In EX2, ME has fallen by more than 68% and EX3 shows a similar result. For the MLT, in the process of introducing the WITM (EX2) and WTFR (EX3), the model results are gradually improved. Especially, in EX3, the ME is decreased by 40% from 0.3 ℃ to 0.12 ℃. WITM can improve the simulated MLD result observably, the E value is increased from 0.32 (EX1) to 0.90 (EX2), while WTFR plays a crucial role in the MLT changing, the E value is increased from 0.49 (EX1) to 0.75 (EX3). This is consistent with the previous results at Argo float WMO 2902666 (Figure 9).  Furthermore, we analyzed the effects of wave-induced mixing on the vertical distribution of temperature. Temperature profiles modeled by the EX2 (Figure 11a) and EX3 (Figure 11b) have been extracted perpendicular to the TC track (shown by a black dashed line in Figure 11b). The difference in temperature between with (EX2, EX3) and without (EX1) wave-induced mixing effects are presented in Figure 11c,d, respectively. The effects of wave-induced mixing show a hamburger-like construction with the upper and lower layer to be a cooling area and the middle layer to be a warming area (Figure 11c, d). The thickness of the upper layer extends to a depth of 50 m which is equivalent to the MLD approximately. These findings are well consistent with the previous study results [31,74]. The effect of wave-induced mixing is more intense on the right side of the track (the black dashed line showed the TC track), while it is weaker on the left side. Moreover, it is shown that the WTFR exerted a more dramatic cooling effect on the right side of track compared with that of the WITM effect.  Furthermore, we analyzed the effects of wave-induced mixing on the vertical distribution of temperature. Temperature profiles modeled by the EX2 (Figure 11a) and EX3 (Figure 11b) have been extracted perpendicular to the TC track (shown by a black dashed line in Figure 11b). The difference in temperature between with (EX2, EX3) and without (EX1) wave-induced mixing effects are presented in Figure 11c,d, respectively. The effects of wave-induced mixing show a hamburger-like construction with the upper and lower layer to be a cooling area and the middle layer to be a warming area (Figure 11c,d). The thickness of the upper layer extends to a depth of 50 m which is equivalent to the MLD approximately. These findings are well consistent with the previous study results [31,74]. The effect of wave-induced mixing is more intense on the right side of the track (the black dashed line showed the TC track), while it is weaker on the left side. Moreover, it is shown that the WTFR exerted a more dramatic cooling effect on the right side of track compared with that of the WITM effect.

Discussion
In this paper, using the MASNUM wave model and the POM, we studied the impact of two different wave-induced mixing effect processes, wave-induced turbulence mixing (WITM) and wave transport flux residue (WTFR), on the upper ocean structure especially during TCs in the Northwest Pacific Ocean. Unlike the traditional breaking wave studies, these two physical processes are both rely on the nonbreaking wave and are therefore independent of the wind stress.
By comparing the simulation results of the year 2016 with observation, it is shown that the original model tends to overestimate the SST, after including a consideration of the two wave effects, the simulation performance is remarkably improved. This improvement is most pronounced in summer and early autumn (Figure 3). We think this may benefit from the seasonal variation of the MLD. Observation data show that the global ocean has a deeper MLD in winter then in summer, while the spring and fall seasons act as transitional periods [75]. With the same intensity of wave-induced vertical mixing, a shallower MLD will make it easier for the colder water from the thermocline to be entrained into the well-mixed ocean surface and ultimately decrease the SST. The mean SST error shows a negative value in January, this may be caused by the inaccuracy of the surface heat fluxes used during the first part of the simulation (from the year 2013 to 2015). At the start of the year 2016, the SST remains a negative error at that stage, though the real time surface heat fluxes data were used. In addition, the three experiments overestimated SST to varying degrees for most of the year 2016 except January. The sea ice component is ignored during the simulation, and this may be one of the possible reasons causing the SST overestimation. The rainfall effect is unconsidered during surface net heat flux calculating which may also affect the SST simulation results.
During the TC period, the upper ocean always shows a significant response to the air forcing, such as SST cooling, MLD deepening and MLT decreasing [25,31,76,77]. Our study shows that with the WITM effect (EX2) the simulated sea surface cooling area and the cooling amplitude are more consistent with observations. The WTFR resulted in a better match with the observation of cooling area distribution. One of the most important features is that, by its action, the left deviation tendency of the maximum cooling area is suppressed ( Figure 5). Here, we give the surface Bsmt distribution ( Figure 12) in this area (indicated as a black dashed box in Figure 5b) on 9~11 July. The Figure shows that the left side of the typhoon track is dominated by positive value, indicating a warming effect which inhibits the left deviation tendency of the maximum cooling area. On the contrary, a negative value is observed on the right and means a cooling effect.

Discussion
In this paper, using the MASNUM wave model and the POM, we studied the impact of two different wave-induced mixing effect processes, wave-induced turbulence mixing (WITM) and wave transport flux residue (WTFR), on the upper ocean structure especially during TCs in the Northwest Pacific Ocean. Unlike the traditional breaking wave studies, these two physical processes are both rely on the nonbreaking wave and are therefore independent of the wind stress.
By comparing the simulation results of the year 2016 with observation, it is shown that the original model tends to overestimate the SST, after including a consideration of the two wave effects, the simulation performance is remarkably improved. This improvement is most pronounced in summer and early autumn (Figure 3). We think this may benefit from the seasonal variation of the MLD. Observation data show that the global ocean has a deeper MLD in winter then in summer, while the spring and fall seasons act as transitional periods [75]. With the same intensity of wave-induced vertical mixing, a shallower MLD will make it easier for the colder water from the thermocline to be entrained into the well-mixed ocean surface and ultimately decrease the SST. The mean SST error shows a negative value in January, this may be caused by the inaccuracy of the surface heat fluxes used during the first part of the simulation (from the year 2013 to 2015). At the start of the year 2016, the SST remains a negative error at that stage, though the real time surface heat fluxes data were used. In addition, the three experiments overestimated SST to varying degrees for most of the year 2016 except January. The sea ice component is ignored during the simulation, and this may be one of the possible reasons causing the SST overestimation. The rainfall effect is unconsidered during surface net heat flux calculating which may also affect the SST simulation results.
During the TC period, the upper ocean always shows a significant response to the air forcing, such as SST cooling, MLD deepening and MLT decreasing [25,31,76,77]. Our study shows that with the WITM effect (EX2) the simulated sea surface cooling area and the cooling amplitude are more consistent with observations. The WTFR resulted in a better match with the observation of cooling area distribution. One of the most important features is that, by its action, the left deviation tendency of the maximum cooling area is suppressed ( Figure 5). Here, we give the surface Bsmt distribution (Figure 12 For the mixed layer affected by TC, due to the intense and quick effect of the WITM, the MLD rapidly deepened to a similar level of observed which is severely underestimated in the original ocean current model (Figure 9a). WTFR has no great influence on the MLD deepening. However, by introducing the WTFR, the model shows a more realistic MLT cooling and warming period compared with that of the other two experiments (Figure 9b). The results show that the MLD return to its pre-storm depth approximately 7 days after the passage of the TC (14 July), and a longer period of almost twenty days is needed for the MLT recovery. This is similar to the findings of Wu and Chen [78], who also argued that the asymmetric effects of entrainment and detrainment on the ML are responsible for the prolonged time required to restore the MLT compared with that of the MLD. Comparing among three experiments, it is shown that the wave-induced mixing effects on the upper ocean temperature have a three-layer structure (Figure 11), which is cooling, warming, and cooling. This effect is more intense on the right side of the TC track. Especially, WTFR has a more dramatic cooling effect on the right side of track compared with that of WITM.
In the wake of TC, regarded as a "relaxation stage", the near-inertial waves are excited and the energy is dispersed in different vertical modes [79]. It is a non-negligible fraction of the wind power input and can be radiated down into the interior ocean [80,81]. The inertial response of the ocean to the TC is an important process for subsurface mixing and even for determining large-scale patterns of ocean circulation (like meridional overturning circulation) [82]. Thus, it is worthwhile to perform some analyses of that aspect here. The response of near-inertial currents is more pronounced to the right side of the TC track because of the resonance between the wind stress vector and inertial oscillations [72,83]. In Figure 13, we show the vertical section of the zonal velocity of upper 1000 m at Argo float 2902666 (Figure 1b, cross mask, on the right side of TC track) as a function of time, which is obtained from the present three numerical experiments. Modeled results show a near-inertial oscillation after the passage of TC, and the phase lead of the velocity oscillations clearly increased with depth, indicating downward radiation of energy. It can be seen that with the introduction of the WITM (EX2, Figure 13b) and WTFR (EX2, Figure 13c), the downward radiation velocity obviously increased. Moreover, EX3 shows the strongest oscillation amplitude among the three experiments. The near-inertial oscillation is gradually decayed after five inertial periods, which is consistent with previous findings [83]. Comparing the results of different experiments, For the mixed layer affected by TC, due to the intense and quick effect of the WITM, the MLD rapidly deepened to a similar level of observed which is severely underestimated in the original ocean current model (Figure 9a). WTFR has no great influence on the MLD deepening. However, by introducing the WTFR, the model shows a more realistic MLT cooling and warming period compared with that of the other two experiments (Figure 9b). The results show that the MLD return to its pre-storm depth approximately 7 days after the passage of the TC (14 July), and a longer period of almost twenty days is needed for the MLT recovery. This is similar to the findings of Wu and Chen [78], who also argued that the asymmetric effects of entrainment and detrainment on the ML are responsible for the prolonged time required to restore the MLT compared with that of the MLD. Comparing among three experiments, it is shown that the wave-induced mixing effects on the upper ocean temperature have a three-layer structure (Figure 11), which is cooling, warming, and cooling. This effect is more intense on the right side of the TC track. Especially, WTFR has a more dramatic cooling effect on the right side of track compared with that of WITM.
In the wake of TC, regarded as a "relaxation stage", the near-inertial waves are excited and the energy is dispersed in different vertical modes [79]. It is a non-negligible fraction of the wind power input and can be radiated down into the interior ocean [80,81]. The inertial response of the ocean to the TC is an important process for subsurface mixing and even for determining large-scale patterns of ocean circulation (like meridional overturning circulation) [82]. Thus, it is worthwhile to perform some analyses of that aspect here. The response of near-inertial currents is more pronounced to the right side of the TC track because of the resonance between the wind stress vector and inertial oscillations [72,83]. In Figure 13, we show the vertical section of the zonal velocity of upper 1000 m at Argo float 2902666 (Figure 1b, cross mask, on the right side of TC track) as a function of time, which is obtained from the present three numerical experiments. Modeled results show a near-inertial oscillation after the passage of TC, and the phase lead of the velocity oscillations clearly increased with depth, indicating downward radiation of energy. It can be seen that with the introduction of the WITM (EX2, Figure 13b) and WTFR (EX2, Figure 13c), the downward radiation velocity obviously increased. Moreover, EX3 shows the strongest oscillation amplitude among the three experiments. The near-inertial oscillation is gradually decayed after five inertial periods, which is consistent with previous findings [83]. Comparing the results of different experiments, EX1 has the shallowest radiation depth, limited in the upper 600 m. With the effect of the wave-induced mixing, the radiation depth is continually increased to almost more than 1000 m in EX3, which is more consistent with many observations [80,84,85]. For instance, Cuypers et al. [80] studied the impact of storm-induced near-inertial internal waves in the upper ocean and found that the near-inertial internal wave packets were identified down to 1000 m, the maximum depth of the measurements. The EX3's result is in better agreement with their findings, which indirectly reflects the effect of WTFR on the vertical momentum distribution and the mixing processes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 21 of 27 EX1 has the shallowest radiation depth, limited in the upper 600 m. With the effect of the wave-induced mixing, the radiation depth is continually increased to almost more than 1000 m in EX3, which is more consistent with many observations [80,84,85]. For instance, Cuypers et al. [80] studied the impact of storm-induced near-inertial internal waves in the upper ocean and found that the near-inertial internal wave packets were identified down to 1000 m, the maximum depth of the measurements. The EX3's result is in better agreement with their findings, which indirectly reflects the effect of WTFR on the vertical momentum distribution and the mixing processes. The inertial response of the ocean to the TC is an important process for subsurface mixing. By experimental comparison, it is found that the wave could also exert a considerable effect on this phenomenon. Generally speaking, the wave-induced mixing not only increased the oscillation amplitude and downward radiation depth but also speed the downward radiation velocity. The near-inertial response is significantly enhanced which is more consistent with the previous studies.
Near the surface, the momentum flux (total stress) transferred from the atmosphere to the ocean. Previous studies [86,87] show that the total stress can be divided into wave-and turbulence-generation components. Similarly, the Reynold flux of the temperature is also the combination of these two parts. Because of the experimental design, we will discuss the distribution of the surface temperature Reynold flux here. The total Reynold flux of temperature in this study can be expressed as follows: where the first term ( ) is the turbulence-induced component and the last two terms ( , )are the wave-induced components. The EX1, EX2, and EX3 contain the first term, the first two terms, and all three terms, respectively. The total temperature Reynold flux of the three experiments (Figure 14a,d, and h) show that the main cooling area is basically located in the right front of the typhoon. Figure 14a also represents the horizontal distribution of the turbulence-induced surface temperature Reynold flux, derived from the EX1. Within the typhoon eye, there exists a near-zero flux region due to the low wind speed [69]. The negative flux region that represents the cooling effect is fanned out toward the typhoon moving direction with the typhoon eye as the center, and the amplitude gradually decreases with the distance from the typhoon eye. In the EX2 and EX3 (Figure 14b,e), the turbulence-induced component is suppressed, and the surface temperature Reynold flux is dominated by the (Figure 14c,f); i.e., the WITM plays a central role in surface cooling that is consistent with our above comparison results (Figure The inertial response of the ocean to the TC is an important process for subsurface mixing. By experimental comparison, it is found that the wave could also exert a considerable effect on this phenomenon. Generally speaking, the wave-induced mixing not only increased the oscillation amplitude and downward radiation depth but also speed the downward radiation velocity. The near-inertial response is significantly enhanced which is more consistent with the previous studies. Near the surface, the momentum flux (total stress) transferred from the atmosphere to the ocean. Previous studies [86,87] show that the total stress can be divided into wave-and turbulence-generation components. Similarly, the Reynold flux of the temperature is also the combination of these two parts. Because of the experimental design, we will discuss the distribution of the surface temperature Reynold flux here. The total Reynold flux of temperature in this study can be expressed as follows: where the first term (F Khc ) is the turbulence-induced component and the last two terms (F Bv , F Bsmt ) are the wave-induced components. The EX1, EX2, and EX3 contain the first term, the first two terms, and all three terms, respectively. The total temperature Reynold flux of the three experiments (Figure 14a,d,h) show that the main cooling area is basically located in the right front of the typhoon. Figure 14a also represents the horizontal distribution of the turbulence-induced surface temperature Reynold flux, derived from the EX1. Within the typhoon eye, there exists a near-zero flux region due to the low wind speed [69]. The negative flux region that represents the cooling effect is fanned out toward the typhoon moving direction with the typhoon eye as the center, and the amplitude gradually decreases with the distance from the typhoon eye. In the EX2 and EX3 (Figure 14b,e), the turbulence-induced component is suppressed, and the surface temperature Reynold flux is dominated by the F Bv (Figure 14c,f); i.e., the WITM plays a central role in surface cooling that is consistent with our above comparison results ( Figure 5b). Comparing to the distribution of F Bv , the negative value area of the F Bsmt (Figure 14g) is relatively small, and the flux is also declined. However, the distribution of F Bsmt is more complex than a simple negative distribution. There are many positive flux areas among the negative and, in particular, an evident positive flux stripe behind the typhoon eye. These special flux distributions can result in a more realistic sea surface cooling area distribution in EX3 (Figure 5c), where the WTFR is considered.
Remote Sens. 2020, 12, x FOR PEER REVIEW 22 of 27 5b). Comparing to the distribution of , the negative value area of the smt (Figure 14g) is relatively small, and the flux is also declined. However, the distribution of smt is more complex than a simple negative distribution. There are many positive flux areas among the negative and, in particular, an evident positive flux stripe behind the typhoon eye. These special flux distributions can result in a more realistic sea surface cooling area distribution in EX3 (Figure 5c), where the WTFR is considered. The exchange of momentum and heat between air and sea is a very complicated process, where many different processes are involved especially during the TC period. The physical conditions are always beyond the simplified structure defined in our modeling. In this study, a wind field and wind pressure model is applied to construct the surface forcing data during the typhoon instead of coupling an atmospheric model such as those in the study [14,55]. Although the comparisons showed that the constructed wind and pressure data have a great consistency with the observation, coupled atmosphere-wave-ocean modeling studies might lead to a more comprehensive understanding of the air-sea interaction. In future studies, we will perform further in-depth research to test the effect of the two physical processes on the upper ocean with a fully coupled system.

Conclusion
We investigated the effect of two wave-induced mixing processes (WITM and WTFR) on the upper ocean structure under both general and TC conditions. Including these two physical processes generally produces a reduced mean error of SST under general condition, especially in summer (June, July, and August) and early autumn (September) time. After the introduction of the The exchange of momentum and heat between air and sea is a very complicated process, where many different processes are involved especially during the TC period. The physical conditions are always beyond the simplified structure defined in our modeling. In this study, a wind field and wind pressure model is applied to construct the surface forcing data during the typhoon instead of coupling an atmospheric model such as those in the study [14,55]. Although the comparisons showed that the constructed wind and pressure data have a great consistency with the observation, coupled atmosphere-wave-ocean modeling studies might lead to a more comprehensive understanding of the air-sea interaction. In future studies, we will perform further in-depth research to test the effect of the two physical processes on the upper ocean with a fully coupled system.

Conclusions
We investigated the effect of two wave-induced mixing processes (WITM and WTFR) on the upper ocean structure under both general and TC conditions. Including these two physical processes generally produces a reduced mean error of SST under general condition, especially in summer (June, July, and August) and early autumn (September) time. After the introduction of the wave effects, the performance of the model to simulate ocean response to TC could also be greatly improved.
Specifically, the WITM can effectively improve the simulation results of the SST cooling area and amplitude and has an important effect on the formation of MLD. WTFR plays a positive role in reproducing the distribution of the SST cooling area and is dominated in the MLT changing. The analyses show that the two physical processes could also affect the ocean near-inertial response to the TC. The wave-induced mixing processes should be properly considered during the upper ocean structure modeling.