Numerical simulation of the wave-induced Stokes drift effect on sea surface temperature in the North Pacific

The study area is located in the Pacific between 0° and 66°N latitude and between 100°E and 80°W longitude (Fig. 1). In the middle of the study area, the sea area is wide and the sea waves are large.

Figure 1

Topographic data used in the model were ETOPO5 data with a resolution of 5’ × 5’. Daily averaged wind data and heat flux data (shortwave radiation, longwave radiation, sensible heat flux and latent heat flux) with a resolution of 1.875° × 1.905° came from the National Centers for Environment Prediction (NCEP). The initial temperature field and the salinity field, as well as the lateral boundary were derived from HYCOM data with a resolution of 0.08° × 0.08°. The Optimum Interpolation Sea Surface Temperature (OISST) is a converged product with a spatial resolution of 25 km and a temporal resolution of 1 day, came from the Advanced Very High Resolution Radiometer (AVHRR). This product is derived from the AVHRR satellite-measured data and includes a large-scale adjustment of satellite biases with respect to the in situ data from ships and buoys and uses the optimal interpolation (OI) algorithm for the calculation (Sun et al. 2007; Reynolds et al. 2007; Yan et al. 2016). The European Center for Medium-Range Weather Forecast (ECMWF) reanalysis product ERA-40 with a resolution of 1.5° distinguishes between wind waves and swells, covers the globe and performs well in the North Pacific (Zheng et al. 2015). Buoy data were obtained from the National Data Buoy Center (NDBC) and included wind direction, wind velocity, significant wave height (SWH), wave period and sea temperature. TAO buoys are observational arrays located to the north and south of the equator in the Pacific Ocean, and their datasets are updated daily and are freely available to researchers (Mcphaden 1995; Kara 2006).

The spectral profile parameterization scheme is considered to be a relatively accurate Stokes drift parameterization scheme (Rascle et al. 2008; Webb & Foxkemper 2011; Rascle & Ardhuin 2013; Webb & Fox-Kemper 2015) and was therefore mainly used in this work, whereas the approximate parameterization schemes were introduced only to compare them with the former. The Stokes drift profile based on a one-dimensional wave spectrum is as follows (Breivik et al. 2016):

where g is the gravitational acceleration; ɷ is the angular frequency; k is the wavenumber; and z is the water depth, where the sea surface is 0 and upward is positive. S_ɷ is the one-dimensional frequency spectrum and is calculated as follows:

where E(ɷ,Ө) is the directional wave spectrum and Ө is the direction in which the wave component is traveling.

The sea surface Stokes drift velocity is as follows:

The Stokes transport is as follows:

For deep water, z → −∞, the profile of Stokes transport becomes

In Equations (1), (3) and (5), the spectral formula (EL spectrum) proposed by Elfouhaily et al. (1997) is used:

where k is the wavenumber and B_l and B_h represent the low-frequency spectrum (k < 10k_p, k_p is the wave wavenumber of the spectral peak) and the high-frequency spectrum (k > 10k_p), respectively. The formula for B_l is as follows:

where α_p is the generalized dimensionless Phillips-Kitaigorodskii equilibrium coefficient related to the longwave anti-wave age parameter Ω, i.e. Ω = U₁₀/c_p, where U₁₀ is the 10 m wind speed, c_p is the phase speed at the peak of the spectrum, c _l is the longwave phase speed, and F _p is the longwave effect function. The terms α_p and F_p are calculated as follows:

where k is the wavenumber of the spectral peak (kp=k0Ω2,p) $\left( {{k}_{p}}=\overset{p}{\mathop{{{k}_{0}}{{\Omega }^{2}},}}\, \right)$and k0=g/U102) ${{k}_{0}}\,=\,{{{g}/{U}\;}_{10}}^{2})$and L_PM and J_P are Pierson-Moskowitz-shaped spectra and peak enhancement functions, respectively. The terms L_PM and J_P are calculated as follows:

where J_P is related to U₁₀. The terms γ and Γ are calculated as follows:

The term B_h is calculated as follows:

where α_m is the generalized dimensionless Phillips-Kitaigorodskii equilibrium coefficient related to short wavelengths and c_m is the phase velocity at k_m ( cm=2g/km=0.23 ${{c}_{m}}=\sqrt{{2g}/{{{k}_{m}}}\;}=0.23$m s⁻¹, and k_m is assumed to be the wavenumber at the capillary gravity peak under the curvature spectrum); c_h is the shortwave phase speed; F_m is the shortwave effect function. The terms α_m and F_m are calculated as follows:

where u_* is the friction velocity of the sea surface.

The relationship between the frequency spectrum S_ɷ(ɷ) and the wavenumber spectrum S_k(k) is written as follows:

where V_g is the group velocity.

The spectral profile is inconvenient as it is numerically challenging and expensive to integrate the wave spectrum at every vertical level. Therefore, an approximate parameterization scheme is more suitable for the calculation of the Stokes drift. Three approximate parameterization schemes can be alternatives to the spectral profile.

(1) The monochromatic profile is given by Breivik et al. (2014)

where u₀ is the sea surface Stokes drift velocity calculated by Equation (3); K is the unit wave velocity vector; and z is the water depth, where the sea surface is 0 and the upward direction is positive. To ensure that the Stokes drift at the sea surface and the total Stokes transport agree with the full spectrum values, the wavenumber must be determined by km=u02Ust ${{k}_{m}}=\frac{{{u}_{0}}}{2{{U}_{st}}}$In this relation, U_st is the Stokes transport calculated by Equation (5).

(2) The exponential profile is given by Breivik et al. (2014):

where k_e = k_m/3. Breivik et al. (2014) consider that C ≈ 8 is the most appropriate. The terms u₀ and z are the same as those in the monochromatic profile.

(3) The Phillips spectrum profile is given by Breivik et al. (2016):

where β = 1 and u₀ and z are the same as those in the monochromatic profile. The term erfc(x) is the complementary error function, where erfc(x) = 1 −erf(x), while erf(x) is the following error function:

When x is large, erf(x) becomes the following:

For a vast number of calculations of the spectral profile, it has been customary to use the monochromatic profile (Skyllingstad & Denbo 1995; Polton et al. 2005; Tamura et al. 2012). The monochromatic profile requires the transport, and the surface Stokes drift velocity is easy to implement and computationally inexpensive, but it leads to an underestimation of the shear near the surface. Breivik et al. (2014) proposed an exponential profile based on the same parameters required to compute the monochromatic profile. This approach gives a more accurate shear with no added numerical cost. It offers a slight improvement over the monochromatic profile, however, it is not a good match for the shear. The Phillips spectrum profile was proposed based on the Phillips spectrum (Breivik et al. 2016). We can compute the profile from the same two parameters as the monochromatic profile, which gives a much better match, especially to represent the shear near the surface. Therefore, the Phillips spectrum profile is the best option.

The wave model WAVEWATCH III, developed mainly by Tolman (1996; 2009), is employed to provide the wave field data in this study. The model domain is set from 70°30'S to 66°N latitude and 0°15'E to 359°45'E longitude. The resolution is 0.5° in latitude and longitude. There are 274 × 720 grid points in the horizontal direction. The model be comes stable after a simulation for one month and the simulated time is from January 2014 to December 2014. The simulate d results include many characteristics of the waves, such as the significant wave height, wave direction, wavelength, etc. The simulated significant wave height was validated, and the Stokes drift velocity was calculated based on the simulated results.

The POM is a three-dimensional baroclinic ocean circulation model based on the original equations developed at Princeton University. The model adopts the σ coordinate in the vertical direction and the C grid as the calculation grid. The horizontal time difference is explicit and the vertical difference is implicit. The model is divided into two parts: the inner part and the outer part. The integration time step is shorter for the outer part and longer for the inner part, and it is solved by the time splitting algorithm (Mellor 1998). The Stony Brook Parallel Ocean Model (SBPOM) used in this study is a parallel version based on the POM model developed by Antoni and Wang (Jordi & Wang 2012). In the calculation, the whole area is divided into a plurality of two-dimensional local domains and the serial coding of each local domain is the same. The horizontal array allocated to each local domain is exchanged between the adjacent local domains through the Message Passing Interface (MPI) to achieve parallel computing (Li & Heap 2011). The SBPOM model has the same simulated area as the wave model, with a horizontal resolution of 0.5° and a vertical coordinate of 40 layers. The SBPOM model was spun up by reanalyzed daily averaged wind data from NCEP for 10 years. Then, the model is driven by the wind, heat flux (shortwave radiation, longwave radiation, sensible heat flux and latent heat flux) and the calculated Stokes drift with a time interval of 6 h. The simulated time is from January 2014 to December 2014.

SBPOM does not take sea waves into account. In this study, the Stokes drift is calculated based on the simulated results from WAVEWATCH III and introduced into the momentum equation of SBPOM in the form of velocity to analyze the effect of the Stokes drift.

For a viscous incompressible fluid, the momentum equation is as follows:

where u→ $\vec{u}$is the velocity vector and p is the pressure. The terms u→ $\vec{u}$and p can be decomposed into u→=u¯+u˜+u′ $\vec{u}\,=\bar{u}+\tilde{u}+{u}'$and p→=p¯+p˜+p′, $\vec{p}=\bar{p}+\tilde{p}+{p}',$respectively, where u is the mean current velocity, _T is the tide time scale, u ~<!-- ~ --> $\tilde{\text{u}}$is the wave-induced motion, T ~<!-- ~ --> $\tilde{T}$is the wave time scale, u' is the turbulent disturbance, and T' is the turbulent time scale, and T¯>T˜>T′. $\bar{T}>\tilde{T}>{T}'.$

After the Euler method is used for time averaging, u→ $\vec{u}$can be written as follows:

where 〈 u→〉 $\left\langle {\vec{u}} \right\rangle $and u ~<!-- ~ --> $\left\{ {\tilde{u}} \right\}$represent the average of the wave scale and the average of the disturbance scale, respectively, which implies that 〈u¯〉=u¯,〈u˜〉=0, $\langle \bar{u}\rangle =\bar{u},\langle \tilde{u}\rangle =0,$ 〈 u′〉=0,{ u¯}=u¯, $\left\langle {{u}'} \right\rangle =0,\left\{ {\bar{u}} \right\}=\bar{u},$and u ~<!-- ~ --> = u ~<!-- ~ --> , u ′ = 0. $\left\{ {\tilde{u}} \right\}=\tilde{u}\, , \,\left\{ {{u'}} \right\}=0.$ Substituting u→=u¯+u˜+u′ $\vec{u}=\bar{u}+\tilde{u}+{u}'$and p→=p¯+p˜+p′ $\vec{p}=\bar{p}+\tilde{p}+{p}'$into Equation (24) produces wave averaging, and the average momentum equation is as follows:

where ∇⋅(〈 −u˜u˜ 〉) $\nabla \cdot \left( \left\langle -\tilde{u}\tilde{u} \right\rangle \right)$−is the wave-induced stress term and ∇<!-- ∇ --> ⋅<!-- ⋅ --> −<!-- − --> u ′ u ′$\nabla \cdot \left( \left\langle -\left\{ {u}'{u}' \right\} \right\rangle \right)$is the Reynolds stress term. Several studies have focused on the wave-induced stress term (Longuet-Higgins 1953; Phillips 1977; Xu & Bowen 1994; McWilliams 1999; 2004; Wang et al. 2015). This term can be written as follows (Wang et al. 2015):

where u_st is the Stokes drift and fẑ × u_st is the Coriolis-Stokes force. The Coriolis-Stokes force is a large-scale wave stress induced by the interaction between the Stokes drift and the vorticity of the planet due to the displacement of the vertical wave after considering the rotation of the Earth (Hasselmann 1970; Weber 1983; Xu & Bowen 1994; Sun et al. 2004). u_st × ɷ is the Stokes-vortex force, which is a large-scale force generated by the interaction between the Stokes drift and a wind-driven horizontal current and is the main source of Langmuir circulation (Craik & Leibovich 1976; 1977). ∇ϕ is the term that adjusts the pressure gradient during the influence of the interaction between the Stokes drift and the horizontal current at the sea surface in terms of pressure (Mc Williams 2004; Smith 1998). ϕ can be written as follows:

Thus, the momentum equation with the Stokes drift has the following form:

Figure 2 shows the annual averaged significant wave height (SWH) distribution in the North Pacific Ocean simulated by WAVEWATCH III. The annual averaged SWH gradually decreases with decreasing latitude; SWH at high latitudes is higher than that at low latitudes (Zheng et al. 2015). The SWH on continental coasts is mostly less than 1.5 m. The SWH in high-value areas is greater than 2.5 m and the maximum value is 4.3 m, which is located in the middle of the North Pacific Ocean, corresponding to the prevailing westerly winds in the range of 30°N to 60°N. In the open ocean area corresponding to the northeast trade winds from 10°N to 30°N, the SWH is in the range of 2 to 2.5 m and the SWH is approximately 1.5 m near the Hawaiian Islands due to the island distribution. From the equator to 10°N, winds are weak in most sea areas and the SWH is small, approximately 1.5 m.

Figure 2

The simulation results for January, April, July and October 2014 were selected to represent winter, spring, summer and autumn, respectively. As shown in Figure 3, the reanalyzed data were used to validate the simulated results for these 4 months. The differences between the WAVEWATCH III-simulated results and the ECMWF SWH data demonstrate that the simulated errors are within the range of 0.5 m for most of the

Figure 3

sea area and thus, the simulation exhibits good performance. The simulated results are the best in July and the simulated error of most sea areas is within the range of 0.3 m. In the near-equatorial areas, the simulated error is within the range of 0.2 m. Compared with the validation data, the simulated results are low in most areas at low latitudes, but the western Pacific simulated result is high in January, with an error of 0.5 m. At high latitudes, the western Pacific simulated result is low, except for July, and the simulated error in April is the largest, reaching −0.6 m. The simulated result in the eastern part of the ocean is high and the simulation error in January is the largest, reaching 0.5 m.

The simulated SWH data are validated by data from 18 National Data Buoy Center (NDBC) buoy stations. The locations of the buoys are shown in Figure 4.

Figure 4

The accuracy of the simulated results is measured by calculating the mean error (ME), the root-mean-square difference (RMSD), the correlation coefficient (R), and the nondimensional skill score (SS) between the simulated and observational values. A negative SS indicates that the simulation is poor, but positive values are common. The larger the value, the better the simulation (Murphy 1995). The statistical measures can be written as follows:

where y_i and x_i represent the simulated values and observational values, respectively; y¯andx¯ $\bar{y}\,\,\text{and}\,\bar{x}$represent the average of the simulated values and observational values, respectively; σ_y and σ_x represent the standard deviations of the simulated values and observational values, respectively; and N represents the number of buoy values in 2014 (Table 1).

Table 2 presents the ME, RMSD, R and SS of the simulated values and observational values in 2014, and the results indicate that the simulated results of WAVEWATCH III are consistent with the observational values of the NDBC buoys. The simulated values are low: the ME of 15 buoys is negative, the maximum error is −13.64 cm and the minimum error is 1.03 cm. The RMSD is in the range of 7.52 to 32.19 cm. Except for the value of 46072, the RMSD values are all below 30 cm, with a minimum of 7.52 cm. The lowest R is

0.8401. Thus, all of the R values are greater than 0.8, indicating that the simulated SWH and the measured SWH are strongly correlated and that the simulation is good. The SS values are positive and greater than 0.7. Except for the values of buoys 46071, 46072, and 52404, the remaining SS values are greater than 0.8, indicating that the simulation is good. The reasons for simulation errors may be that swells are not taken into account or the parameterization scheme of wind stress does not fully reflect the transmission process of wind energy (Shi et al. 2013). Swail & Cox (2000) assessed the 10 m wind speeds from NCEP by using advanced third-generation (3-G) ocean wave prediction models and concluded that this approach performed best. Chelton & Freilich (2005) analyzed the accuracy of the 10 m wind speeds from NCEP and found that there was no evidence of bias in the 10 m wind speeds from NCEP. Therefore, the wave results simulated by WAVEWATCH III driven by wind data from NCEP are acceptable.

The sea surface Stokes drift was calculated with Equation (3) (Fig. 5). The distribution of the monthly averaged sea surface Stokes drift in the North Pacific shows a small spatial difference in July and a large spatial difference in other months. In January, April and October, the surface Stokes drift velocity increases with increasing latitude and the drift velocity of the 30°N to 60°N range is greater than that of the 0° to 30°N range. The highest sea surface Stokes drift velocity is at high latitudes, which is consistent with the conclusions of Tamura (2012), Zhang et al. (2014; 2015) and Wang et al. (2015). Areas with high surface Stokes drift velocity values are distributed north of 30°N and correspond to the position of the prevailing westerly winds. The maximum value is 0.2 m s⁻¹. The winds and waves near the equator (0°–5°N) are weak all year round and the Stokes drift is relatively weak. Therefore, the Stokes drift velocity is mostly less than 0.06 m s⁻¹ and varies slightly depending on the season. At high latitudes, the Stokes drift velocity varies greatly from season to season. In July, the Stokes drift velocity was the lowest and uniform, and the maximum velocity domain was located in the South China Sea, with a maximum value of 0.13 m s⁻¹. The maximum Stokes drift velocities in January, April and October were observed in the high-latitude area, and the area with velocities greater than 0.1 m s⁻¹ was large, especially in January and October, occupying almost all areas north of 40°N. In January and October, local high Stokes drift velocity areas were observed in the Philippine Sea and the South China Sea. In April, there were two areas with high drift velocity, and velocities at low latitudes were much lower than those at high latitudes. The Stokes drift velocity was small in July, mainly due to weak wind speed, small energy input and relatively weak

Figure 5

waves. The Stokes drift velocity is highly seasonal, with higher values in January and lower values in July associated with extratropical cyclones (Bell & Kirtman 2018).

To perform a more intuitive comparison of the sea surface Stokes drift velocity with the surface current velocity, the velocity ratio is defined as follows:

where u_v is the monthly averaged velocity of HYCOM. When b=−0.5, u₀/u_v=31.6%, meaning that the surface Stokes drift velocity is 31.6% of the surface current velocity; thus, b=−0.5 can be considered an indicator of the significant levels of the surface Stokes drift velocity in the ocean circulation.

Figure 6 shows the zonally averaged curve of the monthly mean Stokes drift velocity and the HYCOM velocity as well as the curve of the velocity ratio b. The surface Stokes drift velocity in July is low, and the change with latitudes is not significant. In January, April and October, the surface Stokes drift velocity first increases with increasing latitude and then decreases north of 50°N with increasing latitude. The value of b gradually increases and then decreases with increasing latitude. In the range of 30°N to 40°N, the HYCOM velocity is high and the value of b in this area is small. In the range of 0° to 10°N, the b values are below −0.5, indicating that the surface Stokes drift velocity has a relatively little impact on the surface ocean circulation. North of 10°N, b is greater than −0.5, except in July, when b is less than −0.5 near 35°N. These values indicate that the surface Stokes drift velocity has a relatively large impact on the surface ocean circulation and that the surface Stokes drift velocity is even greater than the surface current velocity. In the range of 40°N to 60°N, b is less than 0 only in July and the surface Stokes drift velocity is slightly lower than the surface current velocity. In April, b approaches 0

Figure 6

and the surface Stokes drift velocity is approximately equal to the surface current velocity, and the impact is slightly significant. In October, b is greater than 0 and the surface Stokes drift velocity exceeds the surface current velocity, and the impact is significant. North of 60°N, as the latitude increases, b begins to decrease, approaching −0.5. Therefore, in most of the North Pacific, the surface Stokes drift has a significant impact on the surface ocean circulation and even dominates at high latitudes.

According to Equation (4), the Stokes drift is vertically integrated to obtain the Stokes transport in the North Pacific Ocean (Fig. 7). The distributions of the Stokes transport and Stokes drift are almost

Figure 7

the same. The high-value area of the Stokes transport coincides with the high-value area of the Stokes drift. The Stokes transport distribution shows a small spatial difference in July and a large spatial difference in other months. In January, April and October, the Stokes transport velocity increases with increasing latitude and the transport velocity of 30°N–55°N is greater than that of 0°N–30°N, except for July. The high Stokes transport velocity is observed at high latitudes, which is consistent with the conclusions of Deng et al. (2007), Zhang et al. (2014; 2015) and Wang et al. (2015). Near the equator (0°–5°N), the Stokes transport velocity is low due to the low Stokes drift and does not vary significantly seasonally. In terms of its size distribution, the transport velocity is higher at high latitudes in January, April and October: the Stokes transport velocity in the area of 30°N–60°N is mostly above 0.35 m² s⁻¹, reaching its maximum value of more than 0.95 m² s⁻¹ around 45°N in January. Peaks occur in the area around 45°N. In July, in most areas of the North Pacific, the distribution of Stokes transport velocity was uniform, with values below 0.3 m² s⁻¹. The seasonal variation in the Stokes transport velocity in the North Pacific is apparent. The Stokes transport velocity in July is low, while the Stokes transport velocity in January, April and October is high, which is consistent with the seasonal variation characteristics of the surface Stokes drift.

Previous research (Deng et al 2007; Zhang et al. 2009; Zhang 2013; Bi 2013) has shown that Stokes transport accounts for a large proportion of ocean transport. McWillianms et al. (1999) and Bi (2013) compared the magnitudes of Stokes transport and Ekman transport and concluded that Stokes transport accounts for up to 50% of the Ekman transport. Polton et al. (2005) defined the Ekman-Stokes number according to the work of Mcwilliams & Restrepo (1999), and this parameter represents the proportion of Stokes transport in the net transport. The Ekman-Stokes number is calculated as follows:

where T_E is the total Ekman transport composed of T_Ex and T_Ey, τ_x and τ_y are the x direction and y direction of the wind stress, ρ_w is seawater density, and f is the Coriolis parameter. E_s can be simplified to an expression involving U₁₀ (Polton et al. 2005):

where U₁₀ is the 10 m wind speed, κ is the Kaman constant, κ = 0.4, and C_D is the drag coefficient. The expression of C_D (Wu 1982) has the following form:

Figure 8 shows the variation in the latitudinal mean Ekman-Stokes number with latitude obtained according to Equation (39). The Ekman-Stokes number is zero at the equator due to a Coriolis value equal to 0, whereas the Ekman-Stokes number in high-latitude areas is greater than 0.1 and higher than that in low-latitude areas. The Stokes transport in the sea near 50°N in January accounts for more than 20% of the net transport, and the Stokes transport in other months is less than 20%. In July, from 20°N to 55°N, the Stokes transport accounts for approximately 10% of the net transport, and the Stokes transport in other months features a peak near 50°N. North of 60°N, the Ekman-Stokes number in January and October decreases gradually, while the Ekman-Stokes number in April and July increases gradually. The influence of Stokes transport on the mixing layer in the ocean is concentrated in the high-latitude sea area, which occupies a large proportion of the net transport. The divergence of seawater caused by transport changes the structure of the mixed layer in the ocean, leading to the enhancement of the mixing effect in the upper ocean, which accounts for the SST changes. Compared to other items in the standard Ekman balance, the Ekman-Stokes number (Polton & Belcher 2007) also reflects the distribution of the relative magnitudes of the Coriolis-Stokes force in the North Pacific. At high latitudes, the Coriolis-Stokes force is relatively large.

Figure 8

E_s clearly increases with increasing wind speed and latitude. At 45°N, when the wind speed is 10 m s⁻¹, E_s = 0.34, indicating that the Stokes transport at mid-high latitudes is not negligible in terms of ocean transport.

Langmuir circulation is generated by the wave-current interaction in the upper ocean and is a longitudinal spiral motion that is axially parallel to the wind direction (Langmuir 1938). The Langmuir circulation strengthens the vertical shear instability of the upper ocean circulation, so the effect of mixing in the upper mixed layer is strengthened and the depth of mixing is increased, and the Langmuir circulation also produces small-scale effects in the upper ocean, e.g. by affecting the turbulent kinetic energy. McWilliams et al. (2004) argued that the increase in turbulent energy caused by the Langmuir circulation is the main mechanism leading to the enhancement of mixing. The Langmuir circulation can enhance the upper turbulence effect and the induced vertical convection can increase the vertical water exchange and play a “rolling” role, which causes the lower layer of cold water to cool the upper layer of water. McWilliams et al. (1997) used the Langmuir number as a measure of the magnitude of the Langmuir turbulence effect and defined it as follows:

where u _* is the sea surface friction velocity, u*=| τ |/ρw ${{u}_{*}}=\sqrt{{\left| \tau \right|}/{{{\rho }_{w}}}\;}$, u₀ is the sea surface Stokes drift velocity, τ is the wind stress, and ρ_w is the seawater density. The Langmuir number indicates the relative magnitude of the sea surface friction velocity and the sea surface Stokes drift and is a measure of the relative magnitude of the Langmuir turbulence effect and the Reynolds effect (McWilliams et al. 1997). Li et al. (2005) concluded that if only the dynamic effect is taken into account when the Langmuir number is less than 0.3, the Langmuir turbulence effect has a mixing effect, and the smaller the Langmuir number, the more dominant is the Langmuir turbulence effect in the ocean.

The monthly averaged Langmuir number in the North Pacific is calculated by using Equation (41) (Fig. 10). The Langmuir number is small at high latitudes

Figure 9

Figure 10

and large at low latitudes. Therefore, the Langmuir turbulence effect is more evident at high latitudes than at low latitudes. In the range of 0° to 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 et al. (2015). In the range of 30°N to 60°N, the Langmuir number is less than 0.3. In July, the Langmuir number in the eastern North Pacific was greater than that in the western North Pacific, and in January, April and October, the Langmuir number was less than 0.2 around 30°N. In January and October, the distribution of the Langmuir number is very similar, the differences between areas with Langmuir numbers less than 0.3 are small, and areas with high Langmuir numbers are distributed along the southern coast of North America.

The Stokes drift decreases exponentially with water depth and the maximum influence depth in the vertical direction is referred to as the Stokes drift influence depth. The profile is as follows (Zhang et al. 2014):

where k is the wavenumber. According to the distribution of the Stokes drift influence depth in the North Pacific (Fig. 11), it appears to be deep on the eastern side and shallow on the western side. The Stokes drift influence depth gradually increases from west to east with longitude, and the results are consistent with those reported by Zhang et al. (2014). This is mainly due to the existence of the “swelling pool”. Swells are stronger in the eastern areas than in the western areas, the average wave period is longer in the eastern areas than in the western areas, and the wavenumber is smaller in the eastern areas than in the western areas. According to Equation (42), the Stokes drift influence depth is larger in the eastern areas than in the western areas (Chen et al. 2002; Zhang et al. 2011; Zhang et al. 2014 ). The Stokes drift influence depth in most areas of the North Pacific is in the range of 3 to 20 m. The zonal difference is small and the latitudinal difference is large. The maximum influence depth can reach 20 m and it occurs on the coast of Mexico in April. The largest difference between the eastern and western Pacific is observed in April and reaches 18 m. The latitudinal change in July is minimal. In the Okhotsk Sea, the Japan Sea, the East China Sea, the South China Sea and the sea near the Philippine Islands, the Stokes drift influence depth is less than 5 m, especially in the coastal areas, where the shallowest depths are only 2 to 3 m. In these areas, the influence depth does not change significantly with the season.

Figure 11

The OISST product was used to validate the simulated SST without the Stokes drift (Fig. 12). When the Stokes drift is not taken into account, the simulated

Figure 12

SST values in most sea areas are high and the error is in the range of ±2°C. In the Japan Sea, the simulated SST is low and the seasonal variation is small. At high latitudes, the simulated SST is high, and the simulation errors in April and October are larger than those in January and July. In the range of 30°N to 45°N, the simulated SST in the middle Pacific is low in January and that in the western Pacific is also low in July, and the simulated SST of the surrounding seas is high. At low latitudes, the simulation error is smaller than that at high latitudes, and the simulation error of most sea areas is within the range of ±0.5°C. In the sea near the equator, the simulated SST in January is high, with a maximum error of 1.5°C, while the simulated SST in October is low, with a maximum error of −1°C.

To study the effect of the Stokes drift on the SST simulation, the Stokes drift calculated using the WAVEWATCH III-simulated results was introduced into the SBPOM model. To analyze the influence of the Stokes drift on the SST simulation in the North Pacific and the difference between the simulated SST with the Stokes drift and that without the Stokes drift, the optimal value α is defined:

Figure 13 shows that the Stokes drift mainly plays the role of cooling in the simulation of SST, especially at high latitudes, and the maximum absolute value of α does not exceed 1.2°C. In the sea area north of 30°N, the Stokes drift has less influence in July than in other months, and the absolute value of α is less than 0.3°C. In January, April and October, the Stokes drift has a greater impact, and the absolute value of α is up to 1.2°C. In January, the area affected by the Stokes drift is mainly in the eastern and western parts of the North Pacific. In April, the area affected by the Stokes drift is mainly in the central and eastern parts of the North Pacific. In October, the area affected by the Stokes drift occupies almost the entire area north of 30°N. South of 30°N, the absolute value of α for the four months is less than 0.6°C and the Stokes drift effect is relatively small. In July, the effect of the Stokes drift on the SST simulation is almost negligible. In January, the simulated SST with the Stokes drift in the eastern North Pacific near the equator is lower by approximately 0.3°C. In April, the Stokes drift causes partial cooling in the central North Pacific, with a temperature reduction of 0.5°C. In October, the simulated SST with the Stokes drift in the eastern North Pacific near the equator rises, with a maximum temperature variation of less than 0.5°C.

Figure 13

The error between the simulated SST with the Stokes drift and OISST data is calculated (Fig. 14). According to Figures 12, 13 and 14, the simulation effect of SST is improved after taking the Stokes drift into account. The simulation error of most sea areas is in the range of ±0.5°C and the simulated

Figure 14

SST is closer to the OISST data. The simulated results in January, April and October are good. In July, the simulated results are good south of 30°N, whereas north of 30°N, the influence of Stokes drift is small and the simulation error is in the range of ±0.8°C. The Stokes transport is prone to producing a divergence, which causes upwelling (Smith 1998). The upwelling enhances the convection in the upper mixed layer, which affects the surface water, and the enhancement of convection between the upper and lower layers of seawater cools the SST. According to Figures 5, 6 and 7, the Stokes drift and Stokes transport in January, April and October are stronger at high latitudes than at low latitudes of the North Pacific, causing a significant drop in temperature. At low latitudes, the Stokes drift and Stokes transport are weak, causing a weaker temperature drop than at high latitudes. The distribution of this cooling is consistent with the distributions of the Stokes drift and Stokes transport. In July, the Stokes drift and Stokes transport are weak throughout the North Pacific, resulting in a weaker cooling of the SST.

The Coriolis-Stokes force, one of the wave-induced stress terms, changes the energy input, the velocity profile and the sea surface current distribution of the mixed layer (Hasselmann et al. 1970), and changes in the current distribution directly affect the energy conservation equation, further affecting the SST. The Stokes drift changes the temperature and density distribution of the upper ocean and the current velocity shear through the Coriolis-Stokes force. The Enhancement of instability causes the sea surface current to change. Changes in the sea surface current lead to changes in transport. The upwelling generated behind the transport lowers the SST. Figures 8 and 9 show that at high latitudes, the Coriolis-Stokes force is strong and the Stokes drift has a strong cooling effect, while at low latitudes, the Coriolis-Stokes force is weak and the Stokes drift has a weak cooling effect. Figure 10 shows that the smaller the Langmuir number, the more obvious the SST cooling effect. As the Stokes drift becomes stronger, the Langmuir number decreases, the Langmuir turbulence effect becomes dominant and the turbulent energy increases. This increased turbulence enhances the convection of the seawater, bringing water from depth to the surface and lowering the SST. In the eastern North Pacific near the equator in October, after taking the Stokes drift into account, the simulated SST rises. This simulated SST rise may be explained by the effect of wave mixing. Changes in the ocean temperature and salinity resulting from the wave mixing cause the pressure gradient force to change and weaken the upwelling. The lower seawater cannot rise to the sea surface. As a result, the simulated SST increases (Song 2006; Qiao 2010; Song 2012; Liu 2017).

Eighteen NDBC buoys and seven TAO buoys presented in Figure 4 are used as validation points. The RMSD, R and SS between the daily averaged buoy time series of SST and the daily averaged simulated SST with and without the Stokes drift in 2014 were calculated (Table 3). Table 4 shows the RMSD, R and SS between

the daily averaged buoy time series of sea subsurface temperature and the daily averaged simulated sea subsurface temperature.

As shown in Table 3, there is a large error between the simulated SST without the Stokes drift and the buoy data. The maximum RMSD is 2.06°C and the R values for four buoys (51000, 21419, 52401, and 52402) are less than 0.8. The R value for 46005 is the largest, 0.97, and the values for all buoys are greater than 0.75. After taking the Stokes drift into account, the RMSD is reduced; the maximum RMSD is 1.87°C, the minimum RMSD is only 0.52°C, and the average RMSD is reduced from 1.08°C to 0.95°C. Except for buoy 51000, the buoys have R values greater than 0.8, and the maximum R value is 0.97. Compared with the simulated results without the Stokes drift, the average R of the result with the Stokes drift increased by 0.01, and the average SS increased from 0.82 to 0.84. Therefore, the introduction of the Stokes drift into the ocean circulation model can improve the accuracy of the SST simulation.

In Table 4, the simulated sea subsurface temperature with the Stokes drift is slightly different from that without the Stokes drift. The RMSD of 11 buoys with the Stokes drift is greater than that without the Stokes drift. The maximum RMSD is 2.11°C. The average RMSD with the Stokes drift is smaller than that without the Stokes drift by only 0.1. The R values of these results are all above 0.75, and the average R of the simulated result with the Stokes drift decreases by 0.1. The average SS of the simulated result with the Stokes drift increases by 0.1. For the simulation of subsurface temperature, similar performance can be obtained with the Stokes drift and without the Stokes drift.

In this study, four buoys located in different areas were selected to analyze the simulated results without the Stokes drift, the simulated results with the Stokes drift and the buoy data (Figure 15). The temperature profile simulated with the Stokes drift is similar to the observational temperature profile of the buoys. Based on a comparison of the vertical simulation results with and without the Stokes drift, even though the Stokes drift depth is only 20 m, it indirectly affects the entire mixed layer through the interaction of the Langmuir circulation and the downwelling (Skyllingstad & Denbo 1995; Li et al. 2005). The reason is that the Coriolis-Stokes force and the Stokes-vortex force result in a reduced mean shear and the Langmuir circulation intensifies the downwelling (Polton & Belcher 2007; Tamura et al. 2012; Li et al. 2013). The temperature profiles in Figure 15 show that after taking the Stokes drift into account, the upper boundary of the thermocline can be determined more accurately, and the analysis of the thermocline characteristics can be improved.

Figure 15

The Stokes drift velocity is calculated by three Stokes drift approximation parameterization schemes, introduced into the SBPOM model in the same way. The four experiments: E-1, E-2, E-3 and E-4 correspond to the spectral profile, the monochromatic profile, the exponential profile and the Phillips spectrum profile, respectively. Using the daily averaged buoy time series from 25 buoys presented in Figure 4 to calculate the RMSD of the simulated results of E-2, E-3, and E-4, the effects of different Stokes drift approximate parameterization schemes are compared and the results are shown in Table 5.

Table 5 shows that there is a difference between the Stokes drift parameterization schemes in the SST simulation. The comparison of the results of the four experiments reveals that the RMSD of E-1 is small, the RMSD of E-2 is greater than that of E-1, and the RMSD for 46005 in E-3 is the only value in E-3 that is smaller than that in E-1. The RMSD for 46005 in E-4 is also smaller than that in E-1, but the RMSD values for the remaining buoys are greater than those in E-1. In terms of the average RMSD, the average RMSD of E-1 is the smallest. The average RMSD of E-4 is the closest to that of E-1, followed by E-3, and finally E-2. The RMSD values for two buoys (46075 and TAO2) in E-4 are greater than those in E-3. The RMSD values for other buoys are lower than those in E-3. The RMSD values of all buoys in E-4 are lower than those in E-2. The RMSD values of E-3 slightly differ from those in E-2, but the RMSD values for 22 buoys in E-3 are less than those in E-2. Of the four experiments, the best simulated results were obtained in E-1, with RMSD of 0.95°C, and the simulated results of E-4 are the closest to those of E-1, with RMSD of 1.02°C. The RMSD of E-3 is 1.05°C, thus following E-4. The RMSD of E-2 is 1.06°C, indicating that the simulation is not as good as that of E-1, E-3, or E-4. In general, the simulated results with the Stokes drift are better than those without the Stokes drift.

The temperature profiles of E-1, E-2, E-3, and E-4 were analyzed by selecting four different buoys from those presented in Figure 15. When analyzing the vertical temperature profiles in the mixed layer, the smallest error between the simulated results and the buoy data is that of E-1, followed by E-4. As illustrated by B, C, and D in Figure 16, compared to E-2 and E-3, E-4 performed best in the vertical temperature simulation in the mixed layer. Below the mixed layer, the simulations of E-2, E-3, E-4 and E-1 are not satisfactory. The improvement of SST simulation with the Stokes drift parameterization schemes is mainly reflected in the mixed layer.

Figure 16

Applying the wave spectrum results of ocean wave models to ocean circulation models is a great challenge. The computational complexity of the spectral parameterization scheme is very high, especially the integration of the wave spectrum on each vertical layer. As evidenced by Table 5 and Figure 16, the Phillips spectrum profile approximate parameterization scheme is the best among the three approximate parameterization schemes. With a lower numerical cost than the spectral profile parameterization scheme, these three approximate parameterization schemes require only the transport and the surface Stokes drift velocity. Therefore, in the case of numerical costs, we can use the Phillips spectrum profile approximate parameterization scheme in the ocean circulation model.