Evaluation of Wave-Ice Parameterization Models in WAVEWATCH III® Along the Coastal Area of the Sea of Okhotsk During Winter

Ocean surface waves tend to be attenuated by interaction with sea ice. In this study, six sea ice models in the third-generation wave model WAVEWATCH III® (WW3) were used to estimate wave fields over the Sea of Okhotsk (SO). The significant wave height (Hs) and mean wave period (Tm) derived from the models were evaluated with open ocean and ice-covered conditions, using SO coastal area buoy observations. The models were validated for a period of 3 years, 2008–2010. Additionally, the impact of sea ice on wave fields was demonstrated by model experiments with and without sea ice. In the open ocean condition, the root-mean square error (RMSE) and correlation coefficient for hourly Hs are 0.3 m and 0.92, and for hourly Tm 0.97 s and 0.8. In contrast, for the ice-covered condition, the averaged RMSE and correlation coefficient from all models are 0.44 m (1.6 s) and 0.8 (0.6) for Hs (Tm), respectively. Therefore, except for the bias, the accuracy of model results for the ice-covered condition is lower than for the open water condition. However, there is a significant difference between the six sea ice models. For Hs, the empirical formula whereby attenuation depends on the frequency relatively agrees with the buoy observation. For Tm, the empirical formula that is a function of Hs is better than those of other simulations. In addition, the simulations with sea ice drastically improved the wave field bias in coastal areas compared to the simulations without sea ice. Moreover, sea ice changed the monthly Hs (Tm) by more than 1 m (3 s) in the northwestern part of the SO, which has a high ice concentration.


INTRODUCTION
The Sea of Okhotsk (SO) is a marginal ice zone (defined as the region of an ice cover that is affected by waves and swell penetrating into the ice from the open ocean) and is the southernmost sea with a seasonal ice cover in the Northern Hemisphere. Accurate forecasts of ocean surface waves in the SO are important for navigation planning of ocean transport because it can help identify hazardous areas and ensure safe shipping routes. In addition to its social importance, ocean waves in sea ice play a part in the interaction between sea ice, the ocean, and the atmosphere, and those attenuations are key to the success of the sea ice-wave coupled model (Roach et al., 2019). In winter, sea ice rapidly extends southeastward from November to March before receding (Figures 1A-G). Sea ice suppresses the wave-wind interaction by reducing fetch. It also modifies the wave dispersion relation, and the wave energy is attenuated through a conservative scattering and non-conservative dissipation phenomenon (Squire, 2020). Although the extent of sea ice in the SO has a large interannual variability, its maximum value has been reported by the Japan Meteorological Agency (JMA) to be decreasing at a rate of 3.9 %/decade 1 . Therefore, it is of great concern that a decrease in sea ice in the SO will result in an increase in the height of ocean surface waves in the future.
(Wavewatch III Development Group (WW3DG), 2019), one of the most widely used third-generation spectral wave models based on the radiative transfer equation for global and regional wave forecasts, implements several parameterizations for waveice interaction. In deep water, when currents are absent, the evaluation of wind-generated ocean waves is governed by: where, N = E/σ is the wave action density spectrum, which is a function of the wave number (k) or relative angular frequency (σ = 2πf ), direction (θ), space (x, y), and time (t), E is the wave energy spectral density, f is the frequency, and C g is the group velocity. For the ice-covered region, the source term on the righthand side of Eq.
(1) is defined as follows: where, S in is the input term by wind, S ds is the dissipation term induced by wave breaking, S nl is the nonlinear interaction term among spectral components, S ice is the wave-ice interaction term, and C i is the ice concentration. Both wind input and dissipation terms (S in and S ds ) are scaled by the open water fraction (1 − C i ), whereas S ice is scaled by the ice concentration. The effects of ice on ocean waves can be presented as a complex wavenumber k = k r + ik i , with the real part k r representing the physical wave number related to the wave length and propagation speeds, producing effects analogous to shoaling and refraction by bathymetry, and the imaginary part k i representing the exponential attenuation coefficient k i = k i (x, y, t, σ) which depends on the location, time, and radian frequency. k i is introduced in the WW3 model as: here, α is the exponential attenuation rate for wave energy, which is twice that of the amplitude (α = 2k i ). The above equation (Eq. 3) is used to calculate the dissipation by ice in WW3, denoted as IC1-5 (except for IC0). IC0 is based on Tolman (2003) and provides simple energy flux blocking depending on the local ice concentration. Thus, IC0 does not treat the effect as "dissipation" via the S ice source term. IC1 allows the user to provide an exponential attenuation rate of amplitude that is uniform in the frequency space (Rogers and Orzech, 2013). IC2 assumes dissipation by friction in the boundary layer below the ice cover (Liu and Mollo-Christensen, 1988). IC3 treats the ice cover as a linear viscoelastic layer based on the model by Wang and Shen (2010). IC4 was introduced by Collins and  and provides the wave energy dissipation by one of several simple, empirical, and parametric forms through direct fitting with field data. IC4 is different from the other models and needs seven empirical formulas denoted as IC4M1-M7. In addition, IC5 uses a viscoelastic model based on Mosig et al. (2015). k i is implemented for source functions in IC1-5. The estimation of k r requires IC2, IC3, and IC5 to provide a new dispersion relation. Descriptions of these models are provided in Supplementary Text 1.
WW3 wave-ice parameterization models were applied in several recent studies on field observations from the Arctic Sea, in regions such as the Barents Sea, Chukchi Sea, and Beaufort Sea (e.g., Cheng et al., 2020;Liu et al., 2020;Nose et al., 2020). Nose et al. (2020) evaluated the uncertainty of waveice parameterization models, using three theoretical models (IC2, IC3, and IC5) and the field observations obtained during November in the Chukchi Sea. In addition, Liu et al. (2020) validated the performance of three wave-ice parameterization models (IC2, IC3, and IC4M1-M4) using the field observation from April to May in the Barents Sea. The results suggest that IC3 and IC4M2 corroborate the observations the most. Although sea ice is expected to have a significant impact on the wave fields, no studies have evaluated the effect of sea ice on the wave field in the SO. This study evaluates the wave fields derived from six wave-ice parameterization models (IC0-5 in WW3) using the buoy observations on the north coast of Hokkaido (see Figure 1H). In this study, the wave fields derived from the models were also evaluated for both open ocean and ice-covered conditions. Moreover, the impact of sea ice on wave fields using model simulations with and without sea ice was also clarified. This study had two advantages over previous studies. The first is the evaluation of a considerable number of six wave-ice parameterization models. Six empirical models for IC4 (IC4M1-M7 except for IC4M5) were also evaluated in this study. The second advantage is the use of time-rich observation data. The SO exhibits periods of open ocean and ice-covered conditions, and the buoy observation fully covers both periods (3 years in this study). Therefore, considering the accuracy of the modeled wave field, it can be reliably used to compare the open ocean and ice-covered conditions in the same region.

Model Design
Two model domains were created using a nesting process for a horizontal resolution of 0.25 • (domain 1) and 0.08 • (domain 2) ( Figure 1G). The outer domain (domain 1) covers the entire SO (42 • -63 • N, 135 • -165 • E). The inner domain (domain 2) was used to validate the wave fields in the coastal area (43 • -48 • N, 141.5 • -146 • E). The directional resolution was 10 • , and the frequency  range was 0.035-1.1 Hz, which was logarithmically discretized into 30 increments. GEBCO, 2020 2 was used to provide the bottom topography and coastlines.
The simulation of domain 1 incorporated 6-hourly surface wind data from the 55-year JMA Reanalysis (JRA55) (Kobayashi et al., 2015). This product is approximately 55 km in latitude and longitude. In addition, the wind data for domain 2 were obtained from the JRA55 dynamic regional downscaling product (DSJRA55) (Kayaba et al., 2016) developed by JMA,  which has a spatial resolution of 5 km and a temporal resolution of 1 h. Daily ice concentration was obtained from NOAA Optimum Interpolation (OI) sea surface temperature (SST) version 2 high-resolution dataset with a 0.25 • × 0.25 • spatial grid (Reynolds et al., 2007). Ice thickness was incorporated from the Climate Forecast System Reanalysis (CFSR) produced by the National Centers for Environmental Prediction (NCEP) (Saha et al., 2010). The CFSR product has a spatial resolution of 0.25 • at the equator, extending to a global 0.5 • beyond the tropics, with a temporal resolution of 6 h. The wind, ice concentration, and thickness data were linearly interpolated to the same spatial grid in the wave simulation of both domains. On the north coast of Hokkaido, a southeastward Soya warm current exists along the coast throughout the year. However, this study does not include the influence of ocean currents in model simulation because the strength of the ocean current fluctuates seasonally, and is weak during winter (Ohshima et al., 2017) which is the focal season of this study. In addition, the buoy observation is located in deep water (see section "Buoy Observation" for observation depth), and the effect of tides is also not considered for our model simulation.
In this study, six models for S ice , IC0, IC1, IC2, IC3, IC4, and IC5 were used (see Supplementary Text 1). In addition, in order to investigate the impact of sea ice on the wave field, simulations that did not incorporate ice concentration were conducted. Hereinafter, the model results without ice concentration are denoted as "Non-ICE." As shown in Supplementary Table 1, the kinematic viscosity (ν) is required for IC2, and the ν and the effective shear modulus (G) are required for IC3 and IC5 as input parameters. Some theoretical ice parameters (ν and G) for the theoretical models (IC2, IC3, and IC5) have been proposed by previous studies (see Supplementary Table 1). In this study, three, four, and two simulation cases were conducted for IC2, IC3, and IC5, respectively, using the theoretical ice parameters summarized in Supplementary Table 1. In addition, various parameters have also been proposed for the binomial fitting of IC4M2 and the step function of IC4M6 based on different observational data (Supplementary Tables 2, 3). Therefore, eight and six simulation cases were performed for IC4M2 and IC4M6, respectively. Although IC4M5 and IC4M6 both provide a step function in the frequency space, IC4M6 has more steps than IC4M5. Therefore, IC4M5 was excluded from validation in this study. A simple diffusive scattering model (denoted as IS1 in WW3) was used for these simulations. Another scattering model (denoted as IS2) was implemented in WW3. However, the difference between the scattering models (i.e., IS1 and IS2) was small compared to the difference between the dissipation models (IC0-IC5) (not shown). For terms S in and S ds , we used both ST4 (Ardhuin et al., 2010;Rascle and Ardhuin, 2013) and ST6 (Rogers et al., 2012;Zieger et al., 2015;Liu et al., 2019). In this study, the model results with ST6 are presented in the main text, while the model results with ST4 are presented in Supplementary Material.
All simulations were performed over a 3-year period from 2008 to 2010. The significant wave height (H s ) and mean wave period (T m01 ) were both recorded every hour during the computation period. To simplify the notation, T m01 is denoted as T m . To obtain the modeled value at the buoy position, we bilinearly interpolated the fields to the buoy position using the surrounding four grid values from domain 2.

Buoy Observation
To validate the model results of the wave field, we used buoy observation data from the Nationwide Ocean Wave Information Network for Ports and Harbors (NOWPHAS) 3 , provided by the Ports and Harbors Bureau, Ministry of Land, Infrastructure, Transport, and Tourism (MLIT). Significant wave height and mean wave period data obtained every 20 min were used for the observation depth of 52.6 m in the Monbetsu (south) Station (red dot in Figure 1H). The distance from the coast of the buoy is 8,200 m. Buoy data for 3 years, from 2008 to 2010, were used, same as the modeling period. In the present study, the observation data were averaged from 20 min to 1 h and compared with the simulation results. Figure 2 shows the frequency distribution of H s and T m observed by the buoy during 2008-2010. The averages for H s and T m are 0.83 m and 4.80 s, respectively (gray bars in Figure 2). The total number of hourly observation data points was 24909. In the present study, to evaluate the model simulations for the ice-covered condition, we utilized values only when the ice concentration in the coastal area (44 • -46 • N, 142.5 • -145.5 • E) around the buoy was 10% or more. The number of observation data points was 3277 for the ice-covered condition (blue bars in Figure 2). The number of data points for the ice-covered conditions is significantly reduced compared to that for all data but it covers the entire observation area (Figure 2).
The sea ice in the SO is a one-year ice type and is thinner than that of the Arctic Sea (e.g., Nihashi et al., 2018). Although different from the buoy observation points of this study, the floe size of sea ice up to 5 m accounts for 90% of the total in the northeastern coast of Hokkaido (Kioka et al., 2020).

Validation With Buoy
The wave fields for the open water condition were evaluated before comparing them with the wave fields for the ice-covered condition. Figure 3 shows the scatter diagrams for the open water condition for the model simulation and the buoy observation. The bias, RMSE, and correlation coefficient for H s were 0.02, 0.3, and 0.92 m, and for T m were 0.14, 0.97, and 0.8 s, respectively. The model results are in close agreement with the buoy observations, and are consistent with the results of a previous study (e.g., Shimura and Mori, 2019). This comparison shows the model results of IC1 because there is no significant difference between the other six models  in the open water condition (not shown). The ST4 model simulation was also calculated with close accuracy to the ST6 model simulation (see Supplementary Figure 1).
As shown in Supplementary Tables 1-3, previous studies have proposed various parameters for the three theoretical models (IC2, IC3, and IC5), and two empirical models (IC4M2 and IC4M6). Therefore, these five models are evaluated prior to the comparison between the six wave-ice parameterization models. There are no remarkable differences in IC2 depending on the theoretical parameter, but there are significant differences in IC3 and IC5 (see Supplementary Table 4 for IC2, Supplementary  Table 5 for IC3, and Supplementary Table 6 for IC5). Hereafter, IC2 is the simulation result of using ν by Liu et al. (1991), and the results are almost the same (Supplementary Table 4). In addition, IC3 and IC5 are the model results based on ν by Liu et al. (2020) and G by Mosig et al. (2015), and these results demonstrated relatively better accuracy (see Supplementary  Tables 5, 6). The empirical parameter of Meylan et al. (2014) was employed for IC4M2 because the model results mostly agree with buoy observations (Supplementary Table 7). The results of IC4M6 are not presented here as it can be almost reproduced by the binomial fitting of IC4M2 (see Supplementary Tables 7, 8).
To visualize the standard deviation (STD), root mean square error (RMSE), and correlation coefficient between the model simulations and NOWPHAS buoy data, Figure 4 displays Taylor diagrams between the two fields (Taylor, 2001). In addition, Table 1 lists the statistical analysis results between the model simulation and buoy observations for H s and T m . Additionally, Supplementary Figures 2, 3 display the scatter diagrams for the model simulations and the buoy observations for H s (T m ). The model accuracies for the ice-covered condition are lower than that for the open water condition except for bias, regardless of the six wave-ice parameterization models (Figures 3, 4 and Table 1). The averaged RMSE and correlation coefficient from all models for the ice-covered condition are 0.44 m and 0.8, for H s , respectively, and 1.6 s and 0.6, respectively, for T m , respectively. The RMSE and correlation coefficient of H s (T m ) for the icecovered condition are 0.14 m and 0.12 (0.63 s and 0.2) and are less accurate compared with those for the open water condition.
For H s , all model simulations indicate a correlation coefficient greater than 0.75, a normalized STD between 0.99 and 1.4, and a normalized RMSE (NRMSE) and RMSE of less than 0.83 and 0.52 m, respectively ( Figure 4A, Supplementary Figure 2, and Table 1). The bias for all H s simulation was within ± 0.3 m except for IC0 (Table 1). In particular, the bias of IC2 and IC4M2 were within 0.03 m, which were simulated with high accuracy ( Table 1). In addition, the RMSE and the correlation coefficients of IC1 and IC4M2 were 0.4 m and >0.81, respectively, and were better than those of other simulations (Table 1), although the normalized STD of both simulations is slightly overestimated as presented in Figure 4A. In contrast, the bias, NRMSE, and RMSE of IC0 were 0.47, 0.82, and 0.51 m, respectively, and were poorly estimated as compared with other simulations (Figure 4A,  Supplementary Figure 2, and Table 1).
Overall, for T m , all simulations provided corresponding correlation coefficients of less than 0.72, which was relatively worse than those of H s , same as the open water condition ( Figure 4B, Supplementary Figure 3, and Table 1). In addition, the differences in the statistical values between simulations were large compared to those of H s (Figure 4 and Table 1). IC2, IC4M1, IC4M3, and IC4M7 provided poor simulation results, as their RMSE and correlation coefficients for T m were greater than 1.62 s and less than 0.6, respectively ( Figure 4B, Supplementary  Figure 3, and Table 1). In contrast, IC1 and IC4M4 yielded simulations with least amount of error and indicated a NRMSE (RMSE) less than 0.75 (1.34 s) and a correlation coefficient of greater than 0.66, although the normalized STDs for both simulations were less than 0.8 and tend to be underestimated ( Figure 4B, Supplementary Figure 3, and Table 1). Moreover, the bias for IC4M2 and IC5 was ± 0.1 s, smaller than those of other simulations ( Table 1).
The statistical results, except as shown above for the bias, depend on the averaging interval. Supplementary Table 9 lists the statistical values between the model simulation and buoy observations for the daily mean. The averaged RMSE for H s and T m from the all simulations is 0.38 m and 1.24 s, respectively. The correlation coefficient for H s and T m is 0.82 and 0.7, respectively.  Compared with those for hourly data, the RMSE and correlation coefficient of H s (T m ) are 0.06 m and 0.02 (0.36 s and 0.1) more accurate, respectively. Moreover, we validated the model simulations as a function of ice concentration (Figure 5). All H s simulations were overestimated at low ice concentrations (C i < 20 %) (Figure 5A). At high ice concentrations (C i > 20%), the trend was dependent on the simulation (Figure 5A). IC0, IC3, and IC4M4 overestimated, while IC4M1 and IC4M3 tended to underestimate (Figure 5A). IC1 and IC4M2 were relatively close to the buoy observations and were simulated with high accuracy, consistent with the comparison results shown in Figure 4A and Table 1. For T m , differences between the simulations became remarkable as the ice concentration increased (Figure 5B). IC4M2, IC4M4, and IC5 were in good agreement with the observations (Figure 5B). IC1, IC4M1, IC4M3, and IC4M7 underestimated, especially for IC4M1 and IC4M7 at C i > 40% (Figure 5B). On the other hand, the T m of IC0 and IC3 were overestimated, regardless of ice concentration ( Figure 5B).
Overall, there were no significant differences in the wave fields between ST4 and ST6 in the buoy location (Supplementary  Figures 4, 5 and Supplementary Table 10). However, the normalized STD for H s with ST4 was remarkably reduced (approximately 0.15) compared with that of ST6 (Supplementary Figure 4A). When examined as a function of ice concentration, H s and T m with ST4 were slightly smaller than those with ST6, but the trend in both simulations remained the same (Supplementary Figure 5).
As shown in Figure 5, the differences between the simulations became significant at high ice concentrations, especially for T m . Figure 6 and Table 2 show the comparison results between the model simulations and the buoy observations for the high icecovered condition (ice concentration > 50%). For H s , IC1 and IC4M2 were relatively agreed with the buoy observations, similar to their comparison results for the ice-covered condition (ice concentration > 10%) shown in Table 1. In addition, the bias, RMSE, and correlation coefficient of IC5 for H s were -0.06, 0.35, and 0.91 m, respectively, and were also close to those of IC1 and IC4M2 simulations ( Figure 6A and Table 2). For T m , overall, the qualitative results in Table 2 were similar to the results shown in Table 1, except for quantitative values. IC4M4 results were better than those of other simulations ( Table 2). Moreover, IC5 results mostly agree with the observations, compared with the other two theoretical models (IC2 and IC3) ( Figure 6B).

Sea Ice Impact for Wave Field
To evaluate the impact of sea ice on ocean waves in the SO, we compared the non-ICE simulations and IC4M2, which are relatively accurate, especially in H s . Figure 7 shows time series and scatter plot of the monthly averaged H s and T m for the simulation and observation at the buoy position. In general, large H s and T m were observed in the coastal areas of Hokkaido during winter (Figures 7A,B). Interestingly, IC4M2 simulations remarkably mitigated the overestimation of H s and T m for non-ICE simulations from January to April, when sea ice existed (i.e., ice concentration is not 0%) (Figures 7A,B). In addition, the improvement of H s and T m for IC4M2 is confirmed from the statistical values (Figures 7C,D). Figure 8 shows the spatial distribution of averaged H s and T m in February from IC4M2, and the differences between IC4M2 and non-ICE. As expected, the wave fields were strongly dependent on the sea ice field, and H s (T m ) became smaller (larger) as the ice concentration increased (Figures 8A,B). The difference in H s (T m ) between IC4M2 and non-ICE is greater (less) than 1 m (3 s) at high ice concentrations (C i > 70%) (Figures 8C,D).

DISCUSSION
As shown in section "Materials and Methods, " various coefficients have been proposed for the binomial fitting of IC4M2 based on different observational data (Supplementary Table 2). In the eight simulation results of IC4M2, the biases of IC4M6H1, WA3  Table 7). In fact, the attenuation rates of IC4M6H1, WA3 UK, and WA3 NIWA were lower (Supplementary Figure 6). In addition, we validated that the accuracy of IC4M1, IC4M3, and IC4M7 was poor (especially in T m ), and the attenuation rate was significantly different from that of IC4M2 (Supplementary  Figure 7). This is probably because the sea ice conditions based on these parameterizations are different from those of the SO. The sea ice thickness is less than 10 cm in the coastal area of Hokkaido (Figure 1). As mentioned in the previous section, 90% of the floe size of sea ice is within 5 m in the /n, where X s is the simulation result, X o is the observation data, and n is the number of data points. The ice-covered condition result shows the average from all model simulations (i.e., ten modeled results in Figure 4 or Table 1). The color shading indicates the standard deviation of ten modeled simulations. The IC1 simulation result is used for the open water condition. 24 h is highlighted by the gray broken line.
coastal area of Hokkaido. In contrast, IC4M1 uses field data with sea ice floe sizes ranging from 20 to 30 m. In addition, IC4M3 is based on an ice thickness between 0.5 and 3 m, which is much thicker than the ice conditions in this study. Moreover, IC4M7 is based on observations of only the pancake ice region, although both pancake and frazil ice may exist in the SO. On the other hand, IC4M2 was in good agreement with observations for H s , and IC4M4 was relatively close to the observations for T m . The IC4M2 (i.e., Meylan et al., 2014) and IC4M4 were based on parameterization data from the same field observation in the Antarctic Sea with ice thickness ranging from 0.5 to 1 m; its ice conditions are relatively close to those of SO. The empirical formula of IC4M4  assumes the attenuation is a function of only H s as shown in Supplementary Text 1. Recently, Kohout et al. (2020) suggested attenuation considers wave period and ice concentration in addition to H s , based on the observation with < 0.5m ice thickness in the Antarctic Sea. In the future, implementation of this model in WW3 is expected to further improve the accuracy of simulation in thin ice thickness areas such as SO.
Although the simulations of IC5 were relatively worse compared with those of the accurate empirical model (IC4M2 for H s , and IC4M4 for T m ), IC5 were better than those of two theoretical models (IC2 and IC3), especially for the high icecovered condition (ice concentration > 50%). The results of these models remarkably depend theoretical parameters (ν for IC2, and ν and G for IC3 and IC5), as shown in Supplementary Tables 4-6 (especially in IC3 and IC5). In this study, constant theoretical parameters were used for IC2, IC3, and IC5. However, these parameters are affected by the ice conditions (e.g., Cheng et al., 2017); thus, these are likely changing both spatially and temporally in the real ocean. Theoretical models have the potential to further improve accuracy in the future, although it is not possible to know the specific types or floe sizes of ice in the entire SO. In other words, the simulation of wave fields under ice-covered conditions depends the accuracy of sea ice used as a model forcing in addition to the accuracy of the model itself, as shown below.
Recently, Nose et al. (2020) revealed that the uncertainty between ice concentration products is greater than the uncertainty between theoretical models (IC2, IC3, and IC5). Thus, it should be noted that our results depend not only on the parameterization for source terms such as S in , S ds , and S ice , but also on the ice concentration used as forcing. In fact, the results of this study are significantly dependent on the temporal resolution of the ice concentration. For example, Figure 9 shows the statistical analysis results for H s and T m as functions of the interval of time averaging. As shown in Supplementary Table 9, the NRMS (correlation coefficient) decreases (increases) as the interval of time averaging increases (Figure 9). However, up to approximately 24 h, the rate of change of NRMS and the correlation for both wave fields with the ice condition are larger than those for the open water condition. This is probably due to the daily (24 h) ice concentration used in this study. In addition, the theoretical models IC2, IC3, and IC5 also depend on the ice thickness. Moreover, differences in wind data may also be one cause of uncertainty in wave fields.

CONCLUSION
In this study, we evaluated six WW3 wave-ice parameterization models (IC0-IC5) using buoy observations located in the southern part of the Sea of Okhotsk for 3 years from 2008 to 2010. In this comparison, H s and T m from the model were evaluated with open ocean and ice-covered conditions. Overall, the accuracy of model results for the ice-covered condition is lower compared to the open water condition, except for the bias. However, in the ice-covered condition, IC4M2 appears to relatively agree with buoy observations for H s , and IC4M4 is closest to the observations for T m . We also clarified the impact of sea ice on wave fields in the SO. In the coastal areas, the simulation with sea ice drastically improved the bias of the wave fields (H s and T m ) compared to that without the simulation without sea ice. In addition, the difference between the simulations with and without sea ice is more than 1 m (3 s) for the monthly mean H s (T m ). The results of the present study can help researchers and engineers perform calculations of wave fields in the Sea of Okhotsk.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.