1 Introduction

Since the middle of the last century, the Arctic sea ice has gradually been melting because of global climate change (Parkinson et al., 1999; Stopa et al., 2016). The Arctic Ocean is a key region in the global climate system, and variations in sea surface temperature (SST) are significant for climate change and ecosystems. The growth of sea fetch has a significant effect on ocean wave dynamics (Xing et al., 2023). As shown in the literature (Sutton and Allen, 1997; Proshutinsky et al., 2001), the SST and sea level have shown increasing trends. However, the characteristics of water temperature in the Arctic Ocean are worthy of further study because variation in water temperature is related to changes in oceanic dynamics (i.e., waves and currents) (Yao et al., 2023). Observations from moored buoys are valuable for oceanographic research. Collection of in situ data can be costly, making it difficult to measure extreme sea states such as the Arctic Ocean and tropical cyclones (TCs). Instead, remote-sensed products from altimeters (Liu et al., 2016) and synthetic aperture radar (SAR) (Shao et al., 2022a) offer valuable data regarding such conditions.

With the development of efficient computational methods and advancements in oceanography theories, simulations of oceanic dynamics by numerical models have been realized, increasing the possibility of long-term climate analysis. Currently, several numerical models, i.e., WAVEWATCH-Ⅲ (WW3) (The WW3DG, 2019; Sun et al., 2021) and Simulation Wave Nearshore (SWAN) (Booij et al., 1999; Wang and Jiang, 2012), are applied for wave simulations, especially in TCs (Xu et al., 2005; Hu et al., 2020). Both models follow the unique principle of the Wave Model (WAM) (The WAMDI Group, 1998). The difference is that WW3 uses a rectangular grid while SWAN uses an unstructured scheme. An unstructured grid is more flexible and has higher grid resolution, which can be realized in specific areas for more accurate simulations of small-scale phenomena. A structured grid can give greater computational efficiency and stability, which allows the usage of more efficient algorithms and data structures, demonstrating that it is suitable for large-scale global simulations. Thus, previously, WW3 was suitable for wave simulation over global oceans (Bi et al., 2015; Zheng et al., 2018; Shao et al., 2021b), and a triangular mesh in SWAN was usually used to illustrate the curved shorelines of coastal waters (Yang et al., 2020). Particularly, alternative parametrizations of ice-wave interactions (Williams et al., 2013; Li et al., 2019) are provided by WW3, allowing for application in polar regions (Hu et al., 2022; Shao et al., 2022b). There is less difference in the latest versions of the two models, i.e., SWAN (version 41.45) and WW3 (version 6.07), as both use rectangular and triangular mesh schemes. Currently, both WW3 and SWAN can use different types of grids. The latest public version of WW3 allows users to select a traditionally rectangular grid, a curvilinear grid, or an unstructured triangular grid. This is also true for the SWAN model. Moreover, the implementation of the unstructured grid in WW3 is now applicable to high-resolution coastal simulations and performs comparably with SWAN (Abdolali et al., 2020).

The baroclinity properties of seawater in the vertical direction result in the difficulty of comprehensively solving the Navier-Stokes equations for ocean circulation. In the 1960s, a low-order precision finite difference method was implemented to approximate these equations (Bryan, 1969). Particularly, structured horizontal grids and the z-level in the vertical direction were utilized in the numerical model. Following this principle, ocean circulation models have been developed by several groups, namely, the Harvard Ocean Prediction System (Robinson and Walstad, 1987) and the Parallel Ocean Program (Jones et al., 2010). However, when applying z-coordinate models to seas with extreme changes in topography, large simulation errors are unavoidable because vertical mass transport is sensitive to the terrain slope at the bottom of the ocean. Terrain-following σ-coordinates are employed in the updated models, such as the Princeton Ocean Model (POM) (Robertson et al., 2001), the Estuarine Coastal and Ocean Model, and the Regional Ocean Modeling System (Shchepetkin and McWilliams, 2005). Furthermore, isopycnic coordinates in the vertical direction have been applied in global circulation models, such as the Miami Isopycnic Ocean Model (Bleck and Smith, 1990) and the HYbrid Coordinate Ocean Model (HYCOM) (Chassignet et al., 2003). The advantage of isopycnic coordinates is the use of optional coordinates, i.e., the σ-coordinates in shallow water and the z-coordinates in the open or stratified upper ocean. In the horizontal direction, most of the above-mentioned models use orthogonal or nonorthogonal curvilinear grids, a scheme that does not fully fit coastal boundaries and thus cannot show the complex and irregular geometries of estuaries, islands, or tidal creeks. Thus, an unstructured grid similar to that in SWAN has been applied in several circulation models, such as the Finite-Volume primitive equation Community Ocean Model (FVCOM) (Chen et al., 2007) and the Semi-implicit Eulerian-Lagrangian Finite-Element (SELFE) (Zhang and Baptista, 2008). Furthermore, the FVCOM is widely used for the simulation of SST with good accuracy.

From a global perspective, the SST distribution is influenced by the water transport associated with ocean circulation. Moreover, the improvement of sea fetch and sea level rise resulting from the sea ice reduction in the Arctic Ocean has led to the growth of waves. The literature (Sun et al., 2021; Yao et al., 2023) has demonstrated that wave-induced effects influenced the SST between 60˚S and 60˚N latitude. The aim of this study was to examine the contribution of waves to the SST in the Arctic Ocean. The waves in 1993 – 2022 were hindcasted by WW3 using European Centre for Medium-Range Weather Forecasts (ECMWF) wind data and HYCOM current field and sea level data. The wave-induced phenomena, i.e., wave breaking and mixing induced by nonbreaking waves resulting in changes in radiation stress and Stokes drift, were included in the wave simulation. The parallel version of the POM, the Stony Brook Parallel Ocean Model (sbPOM) (Jordi and Wang, 2012), was applied for water temperature simulation. Currently, the Arctic Ocean SST can be mainly obtained by satellite observations, in situ surveys, and model simulations. The satellite observations and in situ data have good precision; however, they are difficult to obtain in such extreme ocean states. Meteorological and oceanic models can predict and hindcast the SST in the whole Arctic Ocean at a relatively low cost. Moreover, previous studies have verified that the sbPOM model performs well in SST simulations and has good applicability (Sun et al., 2021; Yao et al., 2023). Therefore, in this current work, the sbPOM model is utilized to simulate SST in the Arctic Ocean. In addition, sensitivity analysis is performed with and without wave-induced terms in the modeling.

This article is structured as follows. Section 2 describes the numerical models, including their forcing fields and settings. Section 3 presents the validation of the simulated waves and the method for calculating the four wave-induced terms. Section 4 highlights the sensitivity analysis results to determine the contribution of waves to the SST in the Arctic Ocean. The conclusions are drawn in Section 5.

2 Descriptions of the Numerical Models

The models and settings of WW3 and sbPOM are briefly described in this section, and the data used in the two models are introduced.

2.1 Model Settings and Data Used in WW3

Based on the principles of the WAM model, the main step in modeling waves in three dimensions using the WW3 model is to solve the wave energy balance equation described below:

$\frac{\partial E}{\partial t}+\nabla \cdot\left[\left(c_{\mathrm{g}}+c\right) \times E\right]=F, $

(1)

where

$ F=F_{\text {in }}+F_{\text {bot }}+F_{\text {nl }}+F_{\text {ice }}+F_{\text {other }}, $

(2)

where E is the wave action density, ∇ is the Hamiltonian operator, c_g is the group speed; c is the velocity of the ocean current, and F is the input and dissipative source terms. The sources of F include the wind-wave growth term F_in, nonlinear wave-wave interaction term F_nl, wave bottom interaction F_bot, wave-ice interaction F_ice, and other terms F_other, i.e., wave-current interactions and wave breaking.

In this work, the computing region of the wave model in the Arctic Ocean was located within 60˚ – 80˚N, 0˚ – 360˚E, and the simulated period was from January 1993 to December 2022. The results of the WW3 model simulation have a spatial resolution of 0.25˚ and a time interval of 1 h. The mode of input/dissipation terms set was ST2 + STAB 2, which is proven to have a better capability for simulating extreme sea states (Sheng et al., 2019). The four wave-wave interactions were calculated by the GMD2 package (Shao et al., 2018), and the ice-wave interactions were calculated by switch IC4_M1 (Shao et al., 2022b), which considers the marginal ice zone area. Furthermore, the wave frequency bin ranges from 0.04118 to 0.7186 at intervals of Δf/f = 0.01, and the direction resolution is utilized as a two-dimensional wave spectrum, which is decomposed into 24 regular azimuthal directions in 15˚ steps. The ECMWF reanalysis (ERA-5) wind field data, HYCOM current field and sea level data, and CMEMS sea ice data were employed as the forcing fields for the WW3 model. A previous study revealed that the input of HYCOM sea level and current data as forcing fields can enhance the accuracy of WW3 model simulation (Hu et al., 2020). The temporal and spatial resolutions of the ERA-5 wind were 1 h and 0.25˚ grid, respectively. The HYCOM sea surface current and sea level have a spatial resolution of 0.08˚ at 3 h intervals, whereas the CMEMS daily average sea ice has a resolution of 1/12˚. The settings and forcing fields of WW3 are summarized in Table 1.

For example, the water depth was taken from the General Bathymetry Chart of the Oceans (GEBCO) using a spatial resolution of 1 km. Fig.1 illustrates a GEBCO water depth map of the modeling region and the ERA-5 wind at 00:00 UTC on January 1, 2022, in the Arctic Ocean. The HYCOM sea surface current and sea level maps at 00:00 UTC and daily average CMEMS sea ice concentration and thickness on January 1, 2022, are presented in Fig.2 and Fig.3, respectively.

2.2 Model Settings and Data Used in sbPOM

The sbPOM is an updated version of the POM and has the advantage of parallel computation (Sun et al., 2021). In the vertical direction, σ-coordinates are used, and the following formula is used to convert the z-coordinates to σ-coordinates:

$ \sigma=\frac{z-\eta}{H+\eta}, $

(3)

where H is the seafloor topography, η is the sea surface height from the seafloor (z = −H) to the surface (z = η), and σ ranges from −1 to 0. In σ-coordinates, the basic equation can be expressed as

$ \frac{\partial D_u}{\partial x}+\frac{\partial D_v}{\partial y}+\frac{\partial \omega}{\partial \sigma}+\frac{\partial \eta}{\partial t}=0, $

(4)

where

$ D=H+\eta, $

(5)

where u, v, and ω are the velocities in the σ-coordinates.

The sbPOM employs different data sources to produce initial fields of SST, i.e., monthly average temperature and salinity data from the CMEMS, the ECMWF reanalysis (ERA-5) wind field, the total heat flux parameters from NCEP (latent heat, sensible heat, long-wave radiation, and short-wave radiation), and the four wave-induced parameters associated with wave breaking, nonbreaking waves, radiation stress, and Stokes drift are simulated by WW3. The simulations from WW3 and CMEMS data have a spatial resolution of 0.25˚. The temporal resolution of the NCEP data is 6 h, and the spatial resolution is 1.905˚ × 1.875˚ in the latitude × longitude direction. It is necessary to determine whether the sbPOM conducts the preprocessing of the initial field data to fit the grid and resolution of the model, i.e., interpolation, smoothing, or other data processing steps. Fig.4 presents the monthly average SST and sea surface salinity maps from CMEMS in January 2022. Similarly, the NCEP reanalysis total heat flux map at 00:00 UTC on January 1, 2022, is displayed in Fig.5, where red triangles refer to the geographic locations of the available Argos buoys. The sbPOM-simulated results have a spatial resolution of 0.25˚ and a time interval of 1 h. The settings and forcing fields for sbPOM are summarized in Table 2.

3 Methodology

In this section, the validation of WW3-simulated waves in 2020 – 2022 against measurements from the HY-2B altimeter is shown, as well as the description of the method used for calculating the four wave-induced terms.

3.1 Validation of the WW3 Simulation

The Haiyang-2B (HY-2B) Chinese marine monitoring satellite was successfully launched in October 2018. Sea surface wind and wave measurement products are officially released by the National Satellite Ocean Application Service. The HY-2 altimeter can measure sea surface height, significant wave height (SWH), and wind speed worldwide with high accuracy (Jiang et al., 2012). Shao et al. (2021a) reported a root mean squared error (RMSE) of 0.29 m for SWH with a correlation (COR) of 0.98 for geophysical data record data. To validate the WW3 simulations, the SWH measurement data for 2020 – 2022 from the HY-2B altimeter were taken. The WW3-simulated SWH map overlaid by the measurements from the HY-2B altimeter at 12:00 UTC on June 15, 2022, is presented in Fig.6a, where the color footprints indicate the SWH from the HY-2B altimeter. There are more than 100000 matchups for validation, as shown in Fig.6b. Statistical analysis obtained an RMSE of 0.45 with a COR of 0.96 and scatter index (SI) of 0.18.

$ SI = \frac{1}{{\bar Y}}\sqrt {\frac{{\sum\nolimits_{i = 1}^n {{{\left[ {({X_i} - \bar X) - ({Y_i} - \bar Y)} \right]}^2}} }}{n}}, $

(6)

where n-elements X_i and X are the values of simulated SWH and its average value, respectively, and Y_i and Y are the values of SWH from HY-2B and its average value, respectively.

In this work, WW3-simulated waves were suitable for calculating the four wave-inducing parameters used in the sbPOM.

3.2 Four Wave-Induced Terms

Inside a unit area of water, the rate of energy dissipation due to wave breaking is controlled by the following equation (Guan et al., 2007):

$ R_{d s}=2.97 \gamma p_w g \beta_*^{-2} \omega_p E, $

(7)

where

$\beta=\frac{g}{u_* \omega_p}, $

(8)

$ E=\int F(k) \mathrm{d} k, $

(9)

$\omega_p=\frac{2 \mathtt{π}}{T_p}, $

(10)

where R_ds refers to the downward input of turbulent kinetic energy flux due to wave breaking on the sea surface, ρ_w refers to the seawater density, g refers to the acceleration due to gravity that is equal to 9.8 m s⁻², γ refers to the energy dissipation per unit white crown as a percentage of the total wave energy dissipation that is equal to 0.1, β refers to the age of the wave, E refers to the density of the spectrum, k refers to the wave velocity vector, F(k) refers to the two-dimensional wave spectrum, ω_p refers to the spectral peak angular frequency, u_∗ refers to the frictional speed (=$\sqrt{C_d U_{10}} $), C_d refers to the drag coefficient, and T_p refers to the dominant wave period.

The parametric equations for the vertical eddy viscosity coefficient K_m and the vertical diffusion coefficient K_h associated with the mixing induced by nonbreaking waves are (Hu and Wang, 2010)

$K_m=\frac{2 a k^2 \lambda}{\mathtt{π} T} \mathrm{e}^{\frac{2 \mathtt{π} z}{\lambda}}, $

(11)

$K_h=\frac{2 P k^2}{g} \delta \beta^3 W^3 \mathrm{e}^{\frac{g z}{\beta^2 W^2}}, $

(12)

where k is the Karman constant, a is the amplitude, which is half the SWH, T is the fluctuation cycle, λ is the mean wavelength, z is the distance from the surface to a certain depth on the sea floor, P is a dimensionless variable related to the Richardson number coefficient, and δ is the wave steepness that is equal to 0.1.

The components of radiation stress (Wu et al., 2011), namely, S_xx, S_yy, and S_xy, are calculated using

$ S_{x x}=k E\left(\frac{k_x k_x}{k_x^2+k_y^2} F_{c s} F_{c c}-F_{s c} F_{s s}\right)+E_D, $

(13)

$ S_{y y}=k E\left(\frac{k_y k_y}{k_x^2+k_y^2} F_{c s} F_{c c}-F_{s c} F_{s s}\right)+E_D, $

(14)

$ S_{x y}=S_{y x}=\sqrt{k_x^2+k_y^2} E \frac{k_x k_x}{k^2} F_{c s} F_{c c}, $

(15)

$ F_{s c}=\frac{\sin h k(z+h)}{\cos h k D}, \quad F_{c c}=\frac{\cos h k(z+h)}{\cos h k D}, $

(16)

$ F_{s s}=\frac{\sin h k(z+h)}{\sin h k D}, \quad F_{c s}=\frac{\cos h k(z+h)}{\sin h k D}, $

(17)

$ E=\frac{1}{16} \rho_w g H_s^2, $

(18)

where H_s is the SWH, and k_x and k_y are the wave numbers in the x- and y-directions, respectively. D = h + η and $ \int_{-h}^{\eta+} E_D \mathrm{d} z=E / 2$ when z ≠ η and E_D = 0.

The formula for mass transport U_s induced by Stokes drift on the sea surface can be expressed as (Mesquita et al., 1992)

$ U_s=U_{s s} e^{\frac{8 \mathtt{π}^2 z}{g T^2}} k,$

(19)

where

$ U_{s s}=\frac{2 \mathtt{π}^3 H_s^2}{g T^3}, $

(20)

where H_s is the SWH, k is the wave number vector of fluctuation, T is the dominant wave period, and z is the water depth (z = 0 at the surface; z > 0 above the surface).

4 Results and Discussion

The SST in the Arctic Ocean is simulated for 1993 – 2022 using the sbPOM with and without the influence of waves simulated by WW3. Measurements from Argo buoys are utilized for validation. Moreover, the contribution of waves to the SST is discussed.

4.1 Validation of the SST Simulation

To investigate the influence of waves on the SST, two types of simulations from sbPOM were analyzed. In this study, there were more than 20000 matchups between SSTs calculated using the sbPOM and Argo floats. Comparison without wave-induced effects is shown in Fig.7a, while Fig.7b displays the results, including the four wave terms. The RMSE was 0.81℃ with wave-induced effects, which is less than the 1.11℃ RMSE reached without waveinduced effects. We used the sbPOM settings in our previous work (Hu et al., 2022) for the water temperature simulation of the Antarctic surface layer. It should be noted that the underestimation of the ERA-5 wind may result in the behavior of systematic deviation in Fig.7a, and this is also found in the SST simulation in the Antarctic Ocean. However, the accuracy of the SST by sbPOM was enhanced by including the four wave terms, i.e., a COR of 0.97 and SI of 0.13. The SST simulations with wave-induced effects are more realizable, probably because waves affect the temperature of surface water by mixing the ocean. This mixing may cause deep cold water to rise to the surface or sink warmer surface water to the depth.

4.2 Discussion

The spatial distribution of the seasonal average difference in SST (values with wave effects minus values without wave effects) was further analyzed in 2022 in spring, summer, autumn, and winter, as shown in Fig.8. A negative difference in SST demonstrates that waves cause a decrease in SST. This is reasonable because the wave will accelerate heat exchange at the sea-air boundary layer. The maximum difference in SST is up to −1℃ in the Baffin Bay and the Greenland and Norwegian Seas in winter, probably caused by the large wave in winter in the North Atlantic, as illustrated in Fig.9. Similarly, the spatial distribution of the seasonal average WW3-simulated SWH in 2022 is presented in Fig.9. Interestingly, the pattern of the difference in SST is consistent with that for waves. Moreover, the difference in SST generally increases with increasing SWH in autumn and winter. The spatial distribution of seasonal average ice concentration from CMEMS in 2022 is shown in Fig.10. The presence of sea ice limits the exchange of heat between seawater and the atmosphere. This results in the relative stability of the SST. During the winter, larger waves in the North Atlantic may expose the surface of the sea to the atmosphere more easily, contributing to faster heat exchange and increasing the difference in SST.

To quantitatively estimate the contribution of waves to the SST in the Arctic Ocean, linear regression analysis was performed. As shown in Fig.11b, there was a significant negative correlation between the difference in SST and SWH in summer with less sea ice, where the coefficient of determination (R² = 0.42) is greater than that in other months. Particularly, R² was quite low in winter, although the SWH was up to 5 m. We believe that this low R² is due to the rapid expansion of sea ice. In other words, the SST in the Arctic Ocean decreases considering the four effects of waves and sea ice freezing; however, the relationship between the SWH and the difference in SST is relatively weak in winter and spring. The interannual variation of SWH and difference in SST during 1993 – 2022 of the study area is shown in Fig.12, where the red line is the monthly-averaged difference in SST (values with wave effects minus values without wave effects) and the black line is the monthly-averaged SWH, indicating that the difference in SST oscillation has a negative response to the SWH.

5 Conclusions

SST is expected to rise under global climate change (El-Geziry, 2021). Moreover, wind stress has a similar increasing tendency. As stated by Yao et al. (2023), the growth of waves can slow the amplitude of sea surface warming. In our work, the objective was to examine the contribution of waves to the SST in the Arctic Ocean due to the enlargement of the sea fetch caused by sea ice melting (Li et al., 2021). Moreover, sensitivity analysis with and without wave terms for 1993 – 2022 was performed to conduct quantitative predictions.

The ERA-5 wind, HYCOM sea surface current, HYCOM sea level, CMEMS sea ice concentration, and CMEMS thickness were utilized to force the WW3, where the switch IC4_M1 calculated the ice-wave interaction (Shao et al., 2022b). The WW3-simulated SWH were compared with those of the HY-2B altimeter in 2020 – 2022, obtaining an RMSE of 0.45 for SWH with a COR of 0.96 and SI of 0.18. Four wave-induced effects, i.e., wave breaking and mixing induced by nonbreaking waves resulting in changes in radiation stress and Stokes drift, were calculated from the WW3-simulated waves. The fields used to initialize the sbPOM included the CMEMS monthly average sea temperature and sea salinity. The ERA-5 wind and NCEP total net heat flux were the forcing fields in the sbPOM. Moreover, sensitivity analysis for sbPOM with and without wave-induced effects was performed. The sbPOM-simulated SST was validated against the Argos measurements, demonstrating an RMSE of 0.81℃ with wave-induced effects less than the RMSE of 1.11℃ achieved without wave-induced effects. This suggests that waves contribute to SST in the Arctic Ocean. The seasonal distribution of the difference in SST (values with wave effects minus values without waves) in 2022 indicated that the pattern of difference in SST was consistent with that for waves. Moreover, a negative correlation (R² = 0.42) between the difference in SST and SWH in summer was observed; however, this behavior was not obvious in winter or spring due to the presence of sea ice. The interannual variation during 1993 – 2022 also indicated a negative relationship between SWH and the difference in SST.

Mesoscale eddies are common in global seas (Kubryakov et al., 2021). In the future, FVCOM-SWAVE, a coupled marine model, will be employed to investigate the interaction between mesoscale eddies and waves in the Arctic Ocean.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Nos. 42076238 and 42376174), and the Natural Science Foundation of Shanghai (No. 23ZR1426900). We thank the National Oceanic and Atmospheric Administration (NOAA) for the provision of the WAVEWATCH-Ⅲ (WW3) model. The original code of the Stony Brook Parallel Ocean Model (sbPOM) was accessed at http://www.ccpo.odu.edu. The wind field from the European Centre for Medium-Range Weather Forecasts (ECMWF) was downloaded via http://www.ecmwf.int. The General Bathymetry Chart of the Oceans (GEBCO) water depth was collected via ftp.edcftp.cr.usgs.gov. Sea surface current and sea level data from HYCOM were collected via https://www.hycom.org. The monthly average temperature, salinity, sea ice concentration, and thickness were from the Copernicus Marine Environment Monitoring Service (CMEMS) and were accessed via https://marine.copernicus.eu.