Optimizing wave simulation to improve the ocean circulation model in the North Pacific Ocean

In this paper, the simulation performance of two source terms, ST4 and ST6, of the wave model WAVEWATCH III on significant wave heights is investigated in the North Pacific Ocean, and the root mean square error (RMSE), correlation coefficient (R) and mean absolute error (E MAE) were compared between them, and all indices showed that the waves was better modelled by using ST4 in the North Pacific region. Based on the ST4 source term scheme, the Stokes drift is calculated from the wave elements simulated by WAVEWATCH III, and is added to the momentum equation of the SBPOM circulation model with the Stokes drift calculated by spectral integration, and it is found that the sea surface temperature (SST) simulated by the former scheme is better than that by the latter scheme by comparing with that of the Argo buoys, in which the RMSE is reduced by 0.018°C, E MAE is reduced by 0.025°C, and R is improved by 0.007. Finally, by comparing the simulated Langmuir numbers of the two schemes, it is found that in most of the sea area (especially in the deep sea), the simulated Langmuir number of the former scheme is smaller than that of the latter scheme by about 0.1, which indicates the former scheme can be used to enhance seawater mixing, thus improving the accuracy of the coupled WAVEWATCH III-SBPOM model.


Introduction
Waves are an integral part of ocean dynamics, and their effects on the upper physical ocean elements have been widely studied [1].Waves are persistent and pervasive in space and time, and they significantly alter the vertical mixing in the upper ocean, which in turn affects the vertical distribution of temperature, salinity, potential vorticity, etc.Therefore, the study of wave-current interactions is of great significance to the study of ocean circulation and even climate change.Zhang [2] added the Stokes drift as an advection term to the temperature advection variation equation of the HYCOM model, and the results showed that the contribution of the Stokes drift to the temperature variation of the mixed layer is in a comparable order of magnitude with that of the mean current.In the real ocean, there are some deviations between the actual observations and the classical Ekman theory, the angle between the observed surface current velocity and the wind stress is smaller than the theoretical value of 45°, and the current velocity decays more rapidly with depth in the real ocean, the fundamental reason for this discrepancy is that the wave action is not taken into account in the Ekman theory [3], and Huang [4] also confirmed that the Coriolis-Stokes forces are responsible for changing the Ekman equilibrium and current velocity profiles within the mixed layer.By comparing the magnitude of the global Stokes transport with that of the Ekman transport, Bifan [5] found that the global Stokes transport can account for up to 50% of the magnitude of the Ekman transport, and McWilliams and Restrepo [6] obtained the same conclusion when estimating the proportion of the Ekman transport accounted for by the Stokes transport.It can be seen that wave-induced mass transport plays an important role in global mass transport, which in turn has an impact on the upper ocean.
It is shown that Stokes drift also plays a role in vertical mixing.Li et al [7] found that Stokes drift can enhance turbulent mixing and increase the energy dissipation rate in the mixing process by studying the effect of Stokes drift on the mixing layer in the upper ocean.Stokes drift plays a role in mixing mainly due to the existence of the Langmuir circulation, which is generated by the interaction of wave currents in the upper ocean and is a kind of longitudinal spiral motion axially parallel to the wind direction [8].Li and Cao et al [9][10] showed through their study that the presence of Langmuir circulation can enhance the vertical shear instability of the upper circulation, which makes the mixed layer mixing intensify and deepen.
It can be seen that with the increasing understanding of the ocean and climate system, we have a deeper understanding of the mechanisms of wave generation and evolution.Therefore, accurate modeling of waves is indispensable for scientific issues such as wave processes, parameterization of waves in the atmosphere-ocean model, and the effect of waves on upper-layer mixing.
So far, the following WAVEWATCH III (Abbreviation: WWIII) typical model source function term schemes are currently available to simulate large-scale and long-duration waves.The TC simulation is based on the one proposed by Chalikov [11], and the dissipation term is divided into long-wave and short-wave components.Hanson [12] used Ocean weather NRAQ wind field data to validate the TC simulation and pointed out that the hindcasting effect of the WWIII is strongly dependent on wave maturity, that is, the poor simulation of surge wave height reduces the overall simulation effect of the model; the BJA scheme is based on the source function term of WAM cycle4 and adopts the wave growth theory mainly based on Miles [13] and Janssen [14], with a further modification by Abdallaand Bidlot [15].The BJA is better in the simulation of both the wave height and the spectral peak period but is also considered that its adopted WAM4 dissipation scheme does not reflect the actual dissipation theory well, also considered that its adoption of the WAM4 dissipation scheme is not a good representation of the actual dissipation theory.The ACC350 scheme, first proposed by Ardhuin [16] and added to the WWIII, adopts the positive part of the WAM4 wind input term and modifies the expression of friction wind speed.A reduction in input energy due to surge dissipation is also incorporated.The dissipation terms include saturation breaking, accumulation breaking and wave turbulence interactions; the ST2 scheme, i.e., the Tolman and Chalikov 1996 scheme, was used to simulate waves in 2014 in the eastern Arabian near-shore waters by nesting SWAN in WWIII in 2016 by Amrutha [17] and the simulation results show that in the near-shore region the ST2 source term simulates the presence of mean wave period of the high sea state overestimation and low sea state underestimation.The ST4 scheme is a further modification based on the ACC350 scheme, where the expression of the surge dissipation term takes into account the effect of wind speed, and also reduces the energy input of the wind-directed waves in some cases due to the presence of surge dissipation.The ST4 scheme used in this paper refers to some of the parameter settings of the TEST451f version of Saha [18]; the ST6 scheme was first applied to the wave-time [19] and SWAN models [20], and later added to the WWIII [21].
It has been observed that both ST4 and ST6 schemes specifically parameterize the dissipation process of the swell and both are better than the others in simulating swell [22].In 2020, Umesh et al [23] assessed the sensitivity of different input dissipation parameters in the WWIII model in the Indian Ocean by comparing deep-water and shallow-water measurements and found that the ST4 source term simulated the best.However, the ST6 scheme attributes the dissipation to the interaction of waves with seawater turbulent processes [24], and Young [25] also pointed out that the two mechanisms are equivalent in the form of expressions, while other possible mechanisms exist.Even so, there are differences between the two schemes in the modeling of surges.

Experimental Scheme
There are two ways of calculating Stokes drifts on sea surfaces, one using wave spectrum calculations [26]: where  is the gravitational acceleration;  is the angular frequency;  is the wave number;  is the water depth, where the sea surface is 0 and upward is positive.  () is the one-dimensional spectrum.
The other is the way of calculating Stokes drift using wave elements [27], which is more widely used and has the following expression: where  ⃗ ⃗ is the wave propagation direction,  is the wave number,  is the fluctuation frequency, and a is the wave amplitude.
To calculate the Stokes drift using the wave elements simulated by the WWIII wave model, it is necessary to transform equation ( 2) and equation (3).Generally speaking, the wave elements are connected through the dispersion relation, and for the small-amplitude gravity waves, the dispersion relation is: ) where h is the water depth, the wave velocity form considering the dispersion relation can be obtained by substituting equation ( 4) into the wave velocity expression c = σ/k and combining it with the equation k = 2π/λ for the wave number: It can also be expressed as the ratio of the wavelength (λ) to the period (T), so equation ( 5) can also be expressed as: Without considering single-frequency deep-water gravity waves, kh → −，tanh(kh) →1, the wave speed and wavelength can then be expressed as a parametric form of the wave period: Thus, through the connection between wave elements, a Stokes drift expressed using the significant wave height Hs and the mean period T can be derived: In this paper, we will take the WWIII wave model, based on the ST4 and ST6 schemes, compare the accuracies of the two source term schemes in the simulation area, select the scheme with the better results to simulate the wave elements, and obtain the Stokes drift based on the Stokes drift equation, and the Stokes drift using the wave spectral integral, and introduce them into the momentum of the SBPOM model, respectively.The effects of the Stokes drift on the SST and the global circulation are calculated in the SBPOM model momentum equations, respectively.Then, the simulated SST of the SBPOM ocean model is compared with the ECMWF data, and the Stokes drift computed in the two different ways is analyzed and added to the SBPOM ocean model, and finally the differences between the simulated sea surface temperature (SST) are compared, and the conclusions are given in the concluding section.

Materials
The simulation area is the North Pacific Ocean, with latitude ranging from 0° to 68°N and longitude ranging from 120°E to 60°W.The bathymetric topography is interpolated from the ETOPO1 topographic data, and in this paper, we adopt the offshore 10-meter-high wind field from NCEP reanalysis data four times a day, with a temporal resolution of 6 hours and a spatial resolution of 1.875° × 1.905°, and the land-water boundary adopts the solid-wall boundary conditions are used for the water-land boundary, and deterministic boundary conditions are used for the open boundary, given directly by the SODA data.The monthly mean temperature and salinity data from SODA are used as the initial field of the model, with a horizontal resolution of 0.5°×0.5°and 50 layers in the vertical direction.

Models and Settings
The WWIII wave model is a full-spectrum, third-generation wave model developed by the Ocean Modeling Group of the National Centres for Environmental Prediction (NCEP), also known as WWIII, which is more efficient and reasonable in describing the physical mechanisms of wind-wave, wavecurrent interactions, and wave-nonlinear interactions.The WWIII model was run from 1 January 2018 to 31 December 2019, and the results of 2019 were taken for analysis, the field output was once every 6 h, and the output wave elements included significant wave height, wave period, wavelength, wave direction, etc., and the output of the significant wave height point was once every 1 h.The Stokes drift velocity was calculated based on the simulation results.In addition, in the WWIII model, the spatial resolution of the wave spectrum is set to 24°×25°, with 24 wave directions and 25 frequencies, and the initial frequency of the model is 0.04118, the computational time step is 900 s, and the minimum time step of the source function is 300 s.Physical phenomena, such as the wave energy input, the nonlinear wave-wave interactions, the breaking of the white crown, and the bottom friction, are all taken into account in the model.The model takes into account the nonlinear wave-wave interaction, white crown breaking, and bottom friction.
The SBPOM (Stony Brook Parallel Ocean Model) model used in this paper is a parallel model developed based on the POM model [28].The shore boundary of its model adopts the solid-wall boundary, and the open boundary adopts the deterministic boundary conditions, which are directly given by the SODA data, with the wind field and heat flux as the upper boundary conditions of the model, and the SODA monthly mean temperature, salinity, and current field data as the initial field of the model.To ensure the stability of the results, the computation period is from 2010 to 2019, and the results of 2019 are taken for analysis, and the simulation results include temperature, salinity and flow field, and then the model is driven by winds, heat fluxes (short-wave radiation, long-wave radiation, sensible heat flux and latent heat flux) and computed Stokes drift with a time interval of 6h.

Wave simulation results
From Fig. 1 and Fig. 2, it can be seen that the significant wave heights simulated by the ST4 and ST6 schemes are more consistent with the ECMWF reanalysis data, and in the high-latitude waters of the North Pacific Ocean, the two schemes have a greater deviation from the ECMWF reanalysis data, with a deviation of 0.5m.The simulated significant wave heights are more consistent with the ECMWF reanalysis data in the western part of the ocean.In the western part of the ocean, there is a larger deviation within 0.6m from the ECMWF reanalysis than in the eastern part of the ocean.Despite the deviation of the modeled RAW heights from the ECMWF reanalysis data, the overall agreement between the two modeled RAW heights and the ECMWF reanalysis data is good.From Fig. 3, it can be seen that the simulation effect of the ST4 scheme is better than that of the ST6 scheme, especially in the areas with large simulation deviations in Fig. 1 and Fig. 2, the simulation accuracy of the ST4 scheme is improved.To further compare the difference between the simulated and measured significant wave heights of the WWIII model under the two schemes, point outputs of the buoy locations were selected as in Fig. 4.

Figure 4. The geographic locations of NDBC
Due to the space limitation of this paper, we selected one buoy each in the nearshore (46089) and deep sea (51002) for the validation of the significant wave height throughout the year, and the results of the comparison between the simulated significant wave heights and the measured values are shown in Fig. 5.Although the simulated values and the measured values have a large difference at some peaks, with the maximum difference of up to 2 m, and the difference is smaller when the significant wave heights are smaller, it can be seen that the trends of the simulated significant wave heights and the measured values have maintained a better consistency whether it is in the nearshore or the deep sea., where the RMSE reflects the deviation between the predicted value and the true value, the mean deviation reflects the overall deviation between the two, and the correlation coefficient embodies the correlation coefficient reflects the degree of correlation between the two.The smaller the RMSE, the smaller the BIAS, and the closer the absolute value of R is to 1, the smaller the error is and the better the model works.The calculation formula is shown below: Where x i and y i denote simulated and observed values, respectively, and y and x denote the mean of the simulated and observed values.
From the comparison, it can be seen that the significant wave heights simulated by the ST4 and ST6 schemes are correlated with the measured NDBC data.Table 1 shows that the correlation coefficients of the significant wave heights simulated by the two schemes are above 0.8, among which the correlation coefficient of the 46066 buoy is the best, with the correlation coefficient exceeding 0.9, which indicates that the simulation results and the measured values of the buoys belong to the strong correlation type.RMSE is within 0.5 m for most of the buoys, except for some buoys where the error is more than 0.5 m.Among them, the 46089 buoy has the largest RMSE, which is more than 0.6 m, and the 46066 buoy has the smallest RMSE, which is 0.267 m.In terms of simulated deviation, the simulated significant wave heights are lower than the measured values, and the simulated deviation of the 46089 buoy reaches 0.5 m.The simulated deviation of the 46069 buoy reaches 0.6 m.In terms of simulation deviation, it can be seen that the overall simulated significant wave height is lower than the measured value, with a deviation of 0.31 m at the 46089 buoy, which may be related to the fact that the forced wind field is smaller than that of the actual wind field, or the parameterization scheme of the wind stress in the model does not fully reflect the wind energy transfer process, thus leading to the deviation of the simulation results.In terms of the mean absolute error, the errors at each buoy are within 0.4 m, with a minimum of 0.12 m.The above comparative analyses show that the simulated significant wave heights of the WWIII wave model and the NDBC buoys match well with the overall trend, the correlation coefficients are high and the errors are small, which suggests that it is plausible to use the NCEP wind field as the driving force of the WWIII model for the significant wave heights at the spatial fixed points.This indicates that the simulated significant wave heights at the spatial fixed points are credible with the NCEP wind field driving the WWIII model.
In addition, it can be calculated that the average RMSE of the significant wave heights at each buoy simulated by the ST4 source term scheme is 0.4089 m, and that of the ST6 source term scheme is 0.4398 m; the  MAE of the ST4 source term scheme is 0.132 m, and that of the ST6 source term scheme is 0.259 m; the average R of the ST4 source term scheme is 0.8507, and that of the ST6 source term scheme is 0.8368.Taken together, the accuracy of the ST4 source term scheme is better than that of the ST6 source term scheme, and therefore the following section will be based on the ST4 source term scheme for the Stokes drift calculation.

SST simulation results
To investigate the effects of the two calculated Stokes drifts on the simulated values of SST in the North Pacific Ocean, the simulation results of the WWIII model were used to add the Stokes drift calculated from the wave elements and the Stokes drift obtained from the integration of wave spectra to the SBPOM model, respectively (for convenience, the study related to the addition of Stokes drift calculated from the wave elements to the SBPOM model is called Scheme 1 and the study related to the addition of Stokes drift from the wave spectral integration to the SBPOM model is called Scheme 2 in the following).From Fig. 6 and Fig. 7, it can be seen that adding the simulated Stokes drift to the SBPOM circulation model can have an obvious cooling effect.Especially in the high-latitude waters of the North Pacific Ocean, since the Stokes drift and Stokes transport are strong throughout most of the year, it is easy to generate radiative dispersion behind the Stokes transport to cause upwelling, which compensates for the surface water and at the same time enhances the convection in the upper mixing layer, and the enhancement of the connection between the upper and lower seawater will further reduce the SST, and in the low-latitude waters, the magnitude of the cooling induced is lower than that of the low-latitude waters.In the low-latitude waters, the cooling caused is lower than that in the high-latitude waters, and the distribution characteristics of this cooling are in good agreement with the distribution of Stokes drift and Stokes transport.It can be seen that the influencing factors of SST are not only limited to the traditional direct heat transfer and the interaction between sea surface air temperature (SAT) and SST [29], but also the Stokes drift through its induced volumetric transport will change the temperature, which in turn will affect on SST.By comparing Fig. 6, Fig. 7, and Fig. 8, it can be seen that the experimental results of Scheme 2 are much improved compared to Experiment 1, with most of the sea area having significant cooling, especially in the middle and high latitudes, where the maximum cooling reaches 0.3 °C.

Discussion
The presence of Langmuir circulation enhances the vertical shear instability of the upper circulation, leading to enhanced mixing and deepening of the upper mixed layer, and Langmuir circulation also affects small-scale effects in the upper ocean, such as turbulent kinetic energy, etc. McWilliams suggests that the increased turbulent kinetic energy term induced by Langmuir circulation is the main mechanism leading to enhanced mixing [30].Due to the additional turbulent kinetic energy generation term in the equations by considering the Stokes drift, the Langmuir circulation can enhance the upper layer turbulence effect, and its induced vertical convection can increase the vertical water exchange, which leads to the upturning of the cold water in the lower layer, resulting in the cooling effect.The Langmuir number is used as a measure of the magnitude of the Langmuir turbulence effect and its expression is defined [30]: ) 1 2 (16) where  represents the Langmuir number,  * = √||/  , is the sea surface friction velocity,  is the sea surface wind stress, and   is the seawater density.
The average Langmuir number of the North Pacific Ocean was simulated by this formula, and as shown in Fig. 11 and Fig. 12, the Langmuir turbulence effect is more obvious in the high-latitude waters.Between 30°N and 60°N, the Langmuir number is basically lower than 0.3, and between 0° and 30°N, the Langmuir number in the eastern North Pacific is large, and that in the western North Pacific is small, which is consistent with the conclusion of Wang Z [31].As can be seen from Fig. 11 and Fig. 12, the Langmuir number is inversely related to the magnitude of SST reduction, the smaller the Langmuir number, the more pronounced the SST cooling.As can be seen from Fig. 13, in most of the North Pacific Ocean, the Langmuir number calculated by Scheme 2 is about 0.1 smaller than that of Scheme 1, and the Stokes drift calculated by Scheme 2 is better able to increase the vertical water exchange, which leads to a more significant SST cooling.

Conclusion
In the North Pacific region, the significant wave heights simulated by the ST4 source term scheme are better than those by the ST6 scheme.The Stokes drift is calculated from the wave elements simulated by the ST4 scheme, and the Stokes drift is also calculated by spectral integration, they are added to the SBPOM in the form of momentum equations, respectively, and the comparisons by Argo buoys show that the cooling effect under the Stokes drift calculated by spectral integration is better than that by wave elements, and the calculated results are closer to the measured values.After comparison with the Argo buoy, it is found that the cooling effect under the Stokes drift calculated by spectral integration is better than that calculated by wave elements, and the calculated results are closer to the measured values.

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Difference between the significant wave height simulated by the ST4 scheme and the ECMWF data (m)

Figure 5 .
Figure 5.Comparison of simulated significant wave heights with measured values This evaluation of the simulation is considered in terms of the root mean square error (RMSE in the following figure), the mean deviation (BIAS in the following figure), the mean absolute deviation ( MAE in the following figure), and the correlation coefficient (R in the following figure), where the RMSE reflects the deviation between the predicted value and the true value, the mean deviation reflects the overall deviation between the two, and the correlation coefficient embodies the correlation coefficient reflects the degree of correlation between the two.The smaller the RMSE, the smaller the BIAS, and

Figure 6 .
Figure 6.Difference between simulated SST with and without considering Stokes drift under Scheme 1 (℃)

Figure 9 .Figure 10 .
Figure 9.The geographic locations of ArgoTo quantitatively analyze the simulation of the sea surface temperature (SST), the simulation results and Argo measurements are plotted as a scatter plot, as shown in Fig.10, and the fitted line (blue line) is compared with y=x (black line), as can be seen in Fig.10(A), the blue line is found to be below the black line, which indicates that the SST simulated in Scheme 1, has a larger error.From Fig.10(B), it

Table 1 .
Comparison of simulated and measured values