Satellite retrieved sea ice concentration uncertainty and its effect on modelling wave evolution in marginal ice zones

Ocean surface waves are known to decay when they interact with sea ice. Wave-ice models implemented in a spectral wave model, e.g., WAVEWATCH III ® (WW3), derive the attenuation coefficient based on several different model ice types, i.e., how the model treats sea ice. In the marginal ice zone (MIZ) with sea ice concentration (SIC) < 1, the wave attenuation is moderated by SIC: we show that subgrid scale processes relating to the SIC and sea ice type heterogeneity in the wave-ice models are missing and the accuracy of SIC plays an important role in the predictability. Satellite retrieved SIC data (or a 5 sea ice model that assimilates them) are often used to force wave-ice models, but these data are known to have uncertainty. To study the effect of SIC uncertainty ∆SIC on modelling MIZ waves during the 2018 R/V Mirai observational campaign in the refreezing Chukchi Sea, a WW3 hindcast experiment was conducted using six satellite retrieved SIC products based on four algorithms applied to SSMIS and AMSR2 data. The results show that ∆SIC can cause considerable wave prediction discrepancies in ice cover. There is evidence that bivariate uncertainty data (model significant wave heights and SIC forcing) 10 are correlated, although off-ice wave growth is more complicated due to the cumulative effect of ∆SIC along an MIZ fetch. The analysis revealed that the effect of ∆SIC can overwhelm the uncertainty arising from the choice of model ice types, i.e., waveice interaction parameterisations. Despite the missing subgrid scale physics relating to the SIC and sea ice type heterogeneity in WW3 wave-ice models—which causes significant modelling uncertainty—this study found that the accuracy of satellite retrieved SIC used as model forcing is the dominant error source of modelling MIZ waves in the refreezing ocean. 15 Copyright statement. TEXT

assist ships in polar waters to circumnavigate hazards such as high winds and waves, collision with perennial sea ice, and sea-spray icing; however, Jung et al. (2016) describe that the existing polar prediction systems need to be urgently enhanced to effectively manage the risks and opportunities associated with growing human activities. In this regard, the Polar Prediction 25 Project (PPP) has contributed to advancing the predictive capabilities. While wave forecasting in polar oceans is still in its early years, the need for advancing wave forecast capacity will only grow in the emerging Arctic Ocean. This paper focuses on the effect of sea ice concentration (SIC) uncertainty on third-generation spectral wave model simulations in and near a marginal ice zone (MIZ). WMO (2014) defines the MIZ as "the region of an ice cover which is affected by waves and swell penetrating into the ice from the open ocean". This study is primarily focused on the MIZ region at the interface between the open ocean 30 and sea ice field.
Documented academic work on wave-ice interactions has a long history dating back as far as Greenhill (1886) (V.A. Squire, 2007;Mosig et al., 2015). When wind waves propagate through/under sea ice cover, the dispersion relation is modified and wave energy is attenuated due to non-conservative dissipation and a conservative scattering phenomenon. Standalone contemporary spectral wave models simulate wave-ice interactions using sea ice as forcing; in this space, the intensive field measure-35 ments of the Arctic Sea State and Boundary Layer Physics Program  have made a solid contribution to the recent advance of The WAVEWATCH III ® Development Group (WW3DG) (2019) wave-ice interaction parameterisations. Rogers et al. (2016); Cheng et al. (2017); Ardhuin et al. (2018);  describe the development and optimisation of the latest WW3 parameterisations for wave evolution in sea ice cover. Despite the progress, Squire (2018); Thomson et al. (2018) qualify accurately quantifying the wave decay and connecting the associated mechanisms over a large domain still 40 remain a challenge because sea ice fields are notoriously heterogeneous; therefore, the wave-ice interaction parameterisation is a source of uncertainty when simulating wave evolution in MIZs. Recent developments of coupled wave-ice-ocean models on a pan-Arctic scale (Boutin et al., 2019;Roach et al., 2019;Zhang et al., 2020) reflect the growing interest in the surface waves' role in the atmosphere, ocean, and sea ice dynamics: perhaps this indicates that advancing the wave-ice interaction physics is becoming a more pertinent issue to broader scientific communities. 45 Besides the model interior, e.g., wave-ice interaction parameterisations, the 0 th order uncertainty pertains to sea ice forcing accuracy such as SIC and sea ice thickness (SIT). In particular, SIC retrieved from satellite radiometers (or sea ice models that assimilate satellite observations) forms the most fundamental input into wave-ice models and should have a profound effect on sea state predictions. Spatial distributions of SIC in the Arctic Ocean can be mapped daily based on satellite microwave radiometry, and they have been the primary source of sea ice trend and climatological studies; however, discrepancies among 50 retrieval algorithms have been a long-known issue, and numerous intercomparison studies have investigated the effects of retrieval algorithms, and to a lesser extent instruments, on SIC estimates (Comiso et al., 1997;Meier, 2005;Andersen et al., 2007;Notz, 2014;Ivanova et al., 2015;Comiso et al., 2017;Chevallier et al., 2017;Roach et al., 2018;Lavergne et al., 2019).
To date, there is no robust validation of any algorithm, so users are urged to understand strengths and weaknesses of the algorithms when using and interpreting the data (Ivanova et al., 2015;Comiso et al., 2017). The long-known SIC discrepancies 55 imply there is uncertainty in the knowledge of true sea ice coverage (Notz, 2014). The uncertainty is potentially greater for MIZs in the refreezing ocean as satellite derived SIC estimates are known to underestimate thin ice less than 35 cm (Heygster et al., 2014;Ivanova et al., 2015). Because the satellite retrieved SIC has uncertainty, the choice of a SIC product a modeller and model developers select is an error source. Since the latest WW3 wave-ice parameterisation developments have all used different sea ice forcing products, understanding the effect of SIC uncertainty on wave predictions is a relevant contribution 60 and is the primary objective of this paper.
The expedition that inspired this study is introduced to close the preliminary section: the R/V Mirai Arctic Ocean observational campaign in the refreezing Chukchi Sea during November 2018 (JAMSTEC, 2018). A 12 day MIZ transect observation was conducted to capture daily changes in the sea ice field and the associated environmental conditions at the same geographical location. The observation showed firsthand how surface waves propagate through a heterogeneous MIZ sea ice field. We 65 began to inquire how the sea ice field heterogeneity may affect a wave-ice model and how the satellite retrieved SIC represents the observed sea ice field, which led to the subject of satellite retrieved SIC. The ensuing Section 2 introduces the methods employed to analyse the wave model uncertainties associated with SIC forcing including the R/V Mirai observation details.
Section 3 discusses the SIC from a wave modelling perspective using the snapshot images of a sea ice field obtained during the MIZ transect observation. A wave hindcast experiment is conducted using various SIC products as forcing for which the 70 model results are compared with limited available in situ wave observations and two independent predictions as described in Section 4. The analysis is extended to the refreezing Chukchi Sea to examine the bivariate uncertainty data (model significant wave heights and SIC forcing) from a physical view point of modelling wave decay and growth, which is discussed in Section 5. In this section, we also investigate the relative significance of the SIC uncertainty compared with the wave-ice interaction parameterisation uncertainty. Section 6 concludes and discusses the study findings.   Figure 1. Shipboard measurements used in this study include surface wind, sea surface temperature (SST), air temperature, and surface wind waves (WM-2 and Piper-C#15). The details of the R/V Mirai measurement systems are provided in Appendix A.

Drifting buoy wave measurement
Two drifting type wave buoys were also deployed during the R/V Mirai observational campaign. One failed within hours, but the other buoy, Piper#13, survived for 19 days after being deployed on 6 November 2018 22:18; it was remotely switched to a sleep mode to preserve battery on 26 November, and the remote connection ceased on 5 December for some unknown reason. Hardware and on-board data processing were mostly the same as Nose et al. (2018) except Piper#13 produced bulk 95 parameters at 15 minute intervals, which were transmitted near real time via Iridium satellite communication. Piper#13 was deployed at 73.32°N, 201.09°E, and its track between 6 and 27 November is presented in Figure 1 overlaid on the NRCS mosaic. The wave height is calculated as H m0 = 4 √ m 0 within the analysed range of a spectrum between the low and high cut-off frequencies, f_low and f_high, respectively. m 0 = f _high f _low S(f )df where S is the variance density spectrum.

100
SIC estimates from Earth-orbiting satellites are an indirect measurement calculated from microwave brightness temperatures.
Although passive microwave radiation has low energy, brightness temperatures between sea ice and open ocean are distinguishable due to the difference in surface emissivity and physical temperatures. Microwave brightness temperatures measured from different frequency channels can account for the spatial and temporal variations of the ocean surface, so retrieval algorithms can be applied to produce SIC field estimates .

105
Since the 1970's, a number of multichannel passive microwave radiometers have been in operation, and the sensors currently in operation (that are most used) for sea ice analysis are SSMIS and AMSR2. The key difference between the two sensors to derive the SIC spatial distribution is footprint resolution as the latter instrument has around 3-4 times higher resolutions for frequencies near 19, 37, and 89 GHz. For these two sensors, a large number of SIC retrieval algorithms have been developed primarily because different algorithms produce considerably different SIC estimates. This is evidenced by a long list of inter-110 comparison studies (Comiso et al., 1997;Meier, 2005;Andersen et al., 2007;Notz, 2014;Ivanova et al., 2015;Comiso et al., 2017;Chevallier et al., 2017;Roach et al., 2018;Lavergne et al., 2019).
A total of eight SIC products were selected for this uncertainty study based on four algorithms applied to SSMIS and AMSR2 data. Hereinafter, uncertainty of satellite derived SIC c i for a set of data products is defined as follows: -NASA-Team-the algorithm has been used for sea ice trend and climatological studies since the beginning of the satellite radiometry era. The SIC data used in this study are the original NASA-Team algorithm applied to SSMIS data and the enhanced NASA-Team2 algorithm applied to AMSR2 data.
-Bootstrap-like NASA-Team, the Bootstrap algorithm has been used for sea ice trend and climatological studies for many years. The SIC data used in this study are the Bootstrap algorithm applied to SSMIS and AMSR2 data.

125
-OSISAF-this algorithm was selected because the highly reputable European Centre for Medium-Range Weather Forecasts (ECMWF) wave model (ECWAM) uses sea ice forcing based on the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system (Donlon et al., 2012) for which SIC is retrieved using the OSISAF algorithm applied to SSMIS data. We also analyse the AMSR2 data in this study.
-ARTIST-Sea-Ice-this algorithm uses the 89 GHz frequency signal to produce high resolution SIC estimates. This 130 algorithm was selected as accurate higher resolution forcing is generally desirable for numerical models. For this product, we only use the AMSR2 data but analyse two different grids: the pan-Arctic data with 6.250 km resolution and the regional Chukchi-Beaufort data with 3.125 km grid resolution.
The principle of all sea ice algorithms as described in Comiso (2007) is that measured radiative flux can be expressed as  (Comiso, 2007). The accuracy of SIC is then dependent on the closeness of tuning brightness temperature tie points to the ice-free and fully ice-covered ocean surface. The selection of frequency channels to derive polarisation ratios or differences (V and H) and gradient ratios (V polarisation) to retrieve SIC also dictates strengths and uncertainties of each algorithm. Technical details of the respective algorithms are described in the Table 1

WAVEWATCH III ® spectral wave model
The effect of SIC uncertainty on MIZ wave predictions was investigated by a hindcast experiment using The Arctic Ocean wave model developed at the University of Tokyo (TodaiWW3-ArCS) based on WW3, which was introduced in Nose et al.
(2018). Third-generation spectral wave models simulate the numerical evolution of waves as energy budgets based on the action density balance equation, The left hand side concerns wave kinematics where N is the wave action density spectrum, which is a function of frequency σ, direction θ, x and y space, and time t, and c describes the propagation velocities in spatial and spectral coordinates. In deep water when neglecting currents, c is the group velocity c g . Source terms are on the right hand side and the ones relevant to this study include the following: the wind input term s wind , the wave dissipation term s dissipation , the non-linear interaction term 150 s non-linear interactions , and the wave-ice interaction term s ice . The sum of these source terms s is expressed based on the following default scaling in ice-covered waters: Specifically to this study, c i relates to the satellite retrieved SIC and s ice to the ice type, i.e., how the model treats sea ice.
The effect of sea ice on waves are represented via the modified dispersion relation σ = σ(k) where |k| = k = k r + ik i . The real 155 part k r is the physical wavenumber and alters the propagation speed of waves in a sea ice field (analogous to effects of shoaling and refraction by bathymetry), and the imaginary part k i is the exponential decay coefficient. k i is introduced in the model as s ice = −2c g k i N for fully ice-covered sea, i.e., c i =1, and the solution to dN dt = s ice is N 0 e −2cgkit . There are five options for treating sea ice in WW3 denoted as IC1-5; c i provides the scaling in the linkage between s ice and ICX as where p n is the sea ice properties, e.g., effective shear modulus and effective viscosity. Therefore, the rate of attenuation depends on the wave period and sea ice properties and is moderated by c i , i.e., N 0 e −2cicgkit .
The wave-ice models implemented in WW3 that calculate k r to model k i are IC2, IC3, and IC5. IC2 calculates dissipation due to basal friction in the boundary layer below an ice sheet, which is modelled as a continuous thin elastic plate based on the work of Liu and Mollo-Christensen (1988). IC3 treats sea ice as a visco-elastic layer based on Wang and Shen (2010) During the R/V Mirai cruise, sea ice in the MIZ was mainly grease, nilas, and pancake ice, so the hindcast experiment was conducted using the IC3 package (with the default parameters) as it has been designed for these ice types (Rogers et al., 2016;Cheng et al., 2017). Scattering is not expected to be the dominant process in this type of ice fields , so it was not considered in the experiment. Regarding SIT forcing, a homogeneous input option with a value of 10 cm was applied; TodaiWW3-ArCS used in this study has a horizontal resolution of 4 km, and its domain covers most of the Pacific side of the Arctic Ocean including the East Siberian, Chukchi, and Beaufort Seas. The model boundaries connected to the seas of the Arctic Ocean was enclosed by ice cover during the November 2018 modelling period (corresponding to the R/V Mirai observation), so nesting was unnecessary. Similar to Rogers et al. (2016), we neglected swell penetration through the Bering Strait. The technical details of TodaiWW3-ArCS's geographical and spectral grids are provided in Appendix B.
ASI-3km and OSISAF-AMSR2 data were excluded for the wave hindcast experiment. The former has a regional coverage 190 that is too small for the TodaiWW3-ArCS domain, and the OSISAF-AMSR2 data have noise in the open ocean, which yield erroneous wave simulation results when they are used as model forcing (as described in Appendix C). The remaining six satellite retrieved SIC products in Table 1 were used as model forcing to examine wave modelling uncertainty, so the ∆c i wave hindcast experiment dataset has {H m0 ci1 , ..., H m0 ci6 } where c ix denotes the satellite retrieved SIC forcing.
Step-like changes of SIC from daily intervals are not ideal as forcing, so the SIC data were interpolated to match the model output frequency of It should be noted that when satellite derived SIC data are used as forcing, the heat and momentum fluxes are distorted in the marine atmospheric boundary layer because the lower atmosphere and the ocean surface are no longer coupled. Inoue et al.
(2011) evaluated surface heat transfer from three reanalysis products by focusing on how the models treat sea ice; they found 200 the treatment of SIC is a key factor for the estimation of surface turbulent heat fluxes. Guest et al. (2018) have elucidated the ice-edge jet generation mechanism based on the in situ data obtained in the refreezing Beaufort Sea. Undoubtedly, altering the sea ice field would feedback to the wind, but this is not captured in this wave hindcast experiment.
3 Sea ice concentration: definition, characteristics, and the use in wave-ice models WMO (2014) defines SIC as "the ratio expressed in tenths describing the amount of the sea surface covered by ice as a fraction 205 of the whole area being considered". The so-called "area considered" presumably varies for different objectives. The length scale of O(10) km may be adequate for sea ice extent climatology, but for wave-ice interactions, the wave provides a scale in a phase-resolved sense. Satellite retrieved SIC represents the fraction of ice-covered water over a large area, sufficiently large enough that the SIC represents a property of a continuum. In reality, the sea ice in the MIZ is granular, and ice floes jam due to horizontal convergence by Langmuir circulation, internal waves, and wind variability, resulting in a formation of features such 210 as ice bands and wind streaks-with which waves likely interact distinctively.
On 14 November 2018 during the MIZ transect observation, R/V Mirai encountered moderate on-ice waves with an H m0 up to around 2.00 m propagating towards the ice edge (this H m0 estimate is consistent from both the shipboard wave data described in Appendix A and hindcast models as shown later). Figure 2 presents a series of snapshot images of the sea ice field during the encounter. R/V Mirai traversed over 10 km in the MIZ from the ice edge, and each image area extends at scaled according to dN dt = c i s ice (Equation 4), the subgrid scale physics is completely missing. It is plausible the subgrid scale distribution of SIC and ice types can be treated in a stochastic manner to provide meaningful mean values to the grid-scale model. On the other hand, SIC c i also affects the WW3 wave-ice model by means of scaling (Equation 4). Figure 2 shows SIC 220 data from eight satellite retrieved products described in Section 2.2 during this event. The SIC estimates interpolated at the R/V Mirai positions largely deviate among the products, characterising the uncertainty of the satellite retrieved SIC. Moreover, the entire time series of the 12 day MIZ transect observation depicts ∆c i is persistent ( Figures A4 to A6 of Appendix D).
Hereafter, we show how large the effect of ∆c i on modelling MIZ waves can be, so much so that it overwhelms the choice of s ice , e.g., ICX. 225 4 ∆c i effects on wave modelling at the observation sites Our goal is to understand ∆c i effects on wave-ice models, but adequate model accuracy, at least qualitatively, is needed for a meaningful uncertainty analysis. Because we lack a sufficient duration of robust in situ wave measurements, two independent numerical wave models that produce predictions in the Arctic Ocean, namely ERA5 ECWAM and the Arctic Monitoring and  Furthermore, when the open ocean sea state is energetic, daily peak ∆H m0 occurs as R/V Mirai sailed into and out of the 250 ice-covered water, indicating the representation of an ice edge in the model forcing plays an important role.
Regarding the comparison of three different base models, ERA5 ECWAM consistently has a positive bias compared with other models except for the on-ice wave event on 14 November when they all agree reasonably well. For this event, both WM-2 and Piper-C#15 have comparable estimates of the measured peak H m0 as R/V Mirai was sailing out of the ice cover, which was around 2.00 m. The ERA5 ECWAM positive bias compared with other models is exacerbated when the   Figure A7 of Appendix D. The buoy data are discussed first. It measured H m0 > 1.00 m for two days after being deployed at the finals hours on 6 November 2018. H m0 tapered off to around 0.20 m, but briefly rose to 0.60 m when the wind speed increased at around 12 November 2018 00:00. Then, it drifted in ice cover with less wave penetration; 270 no measurable wave signals propagated to Piper#13 as the measured spectra indicate instrument noise except during the on-ice wave event between the late hours on 14 November and the early hours of 15 November. Peak wave periods, T p , which is a frequency inverse corresponding to the maximum variance density, were consistently around 9 s; this is likely a true wave signal even though H m0 was only 0.05 m.
Regarding the wave hindcast, TodaiWW3-ArCS ∆H m0 was the largest in open water and decreased with mean(H m0 ). In 275 general, Piper#13 H m0 during 7 and 8 November are underestimated by all numerical models. When the wave energy tapered off between 8-11 November, both ERA5 ECWAM and ARCMFC model H m0 are overestimated. However, the episodic increase of wave energy on 11 and 12 November is reproduced in these models, albeit with a positive bias, whereas the TodaiWW3-ArCS simulations did not show any increase in the H m0 at the Piper#13 location. There are no data shown for the ERA5 ECWAM and ARCMFC models after 12 November at Piper#13 presumably due to the ice masks. There are three 280 occasions when all TodaiWW3-ArCS simulations slightly overestimates H m0 compared with the buoy data, which indicate the model attenuation rates may be too weak or SIC forcing is inaccurate. 5 ∆c i and wave modelling in the refreezing Chukchi Sea The mosaic depicts sea ice edges have a wavy, but highly nonlinear, jagged form. For most of the ice edges, Figure 5 shows OSISAF-SSMIS derived 0.15 contours are smoother whereas the BST-AMSR2 and ASI-3km contours appear to follow the sea ice edge concavity and convexity with a varying degree of closeness; ASI-3km appears to be qualitatively more consistent with the ice edges detected in the NRCS data. Figure  The above suggests the satellite retrieved SIC uncertainty can be considerable on a regional scale. The preceding Section 4 demonstrated that simply changing SIC forcing alone produces considerable ∆H m0 in the wave hindcast experiment at the observation sites. In this section, we extend the wave hindcast analysis to the refreezing Chukchi Sea. For the on-ice wave case, relatively strong small-scale south westerly winds as depicted in Figure 6a generated waves with an H m0 of about 2.00 m propagating towards the ice edge. When on-ice waves encounter ice cover, rapid attenuation is expected within O(10) km Squire, 2018), so the ice edge locations and the SIC variability near it affect the model simulation of wave decay. For the selected off-ice wave case, a low pressure system over Alaska and a high pressure system north west of the Chukchi Sea generated sustained north easterly winds over much of the Chukchi Sea as depicted in Figure 6b, which generated open water H m0 > 3.00 m. In this case, the ice edge and SIC field determine the fetch on which waves are 310 generated, and as such, ∆c i introduces ∆H m0 .

On-and off-ice wave evolution in the refreezing Chukchi Sea MIZs
Owing to the combination of the non-homogeneous nature of wind that generates waves and the SIC field heterogeneity, there was no statistical association for the bivariate uncertainty data (∆H m0 and ∆c i ) even when the H m0 was normalised to wind forcing. In an attempt to elucidate the model uncertainties in the context of physical processes, a scatter plot is produced for both cases with the following visualisation technique: marker sizes are scaled to mean(H m0 ) as a bubble plot and each 315 marker is colour coded according to mean(c i ) like a colour-coded scatter plot. The former aims to emphasise the model data near the ice edge where the effects of wave attenuation and growth associated with ∆c i are most prominent, and the colour coded markers indicate mean(c i ) among all forcing. For simplicity, we refer to these figures herein as an enhanced scatter plot. Figure 7a depicts the spatial distribution of the model ∆H m0 for the on-ice wave case with 0.01, 0.50, and 0.85 mean(c i ) contours. Not all ice edges are aligned to the on-ice wind orientation because of the ice edge geometry. On-ice wave analysis 320 was, therefore, conducted for a strip of the model grid points roughly 100 km in width along the south westerly on-ice wind orientation as depicted in a grey dashed quadrilateral. An enhanced scatter plot shown in Figure 7b depicts bivariate uncertainty data corresponding to the model grid points along the on-ice wind orientation. There is a strong indication of a correlation between the ∆c i and model ∆H m0 for this on-ice wave case. The correlated uncertainties imply that as on-ice waves approach and decay due to wave-ice interactions, larger discrepancies in the representation of SIC as forcing causes greater model  Analysis of the off-ice case is carried out in a similar manner. The data bound by Quadrilateral 2 as shown in Figure 8a reflect 335 the model grid points with the north easterly off-ice wind conditions. Although the correlation is not as high as the on-ice wave case (with higher scatter for ∆c i > 0.10), the enhanced scatter plot in Figure 8c shows that the bivariate uncertainty data are correlated for the off-ice wave case as well. Analogous to the on-ice wave case, high ∆H m0 can occur near the ice edge when only one or two simulations have the SIC forcing representing open water conditions, while the wave growth is suppressed for the remaining simulations due to higher SIC. This effect is depicted in Figure 9c, which shows ∆H m0 and ∆c i along a transect 340 of the Quadrilateral 2 long axis. Along this transect, ASI-6km and BST-AMSR2 have the most north east ice edge, and the waves rapidly grow under the strong north easterly wind forcing whereas the higher SIC of the other simulations suppress the wave growth. A distinct difference between the off-and on-ice wave cases regarding the ∆c i effect on wave-ice models is that  Off-ice wave evolution is a complex process because the fetch is not only controlled by the location of the ice edge, but also wave-ice interactions as implemented in WW3. The current numerical approach to simulate wind pumping energy into waves in ice cover is dictated by c i because waves grow when (1 − c i )(s wind + s dissipation ) > c i s ice . Whether wave evolution in ice cover follows the Equation 3 scaling has been discussed in Rogers et al. (2016); Thomson et al. (2018); the latter cites Li et al. (2017) who confirmed wind input to high frequency wave energy in the Antarctic Ocean. The off-ice ∆H m0 is apparently also 350 influenced by the cumulative effect of the ∆c i along the fetch distance affected by the wave-ice interaction parameterisations as implemented in WW3.
Lastly, for both on-and off-ice wave cases, significant ∆H m0 extends to the waters where the wind forcing is orientated along the ice edge; so the model data are briefly examined in the region of MIZs north east of Wrangel Island, which is shown as Quadrilateral 1 in Figure 8a along the sea ice edge and north easterly wind forcing orientation. This region has considerable 355 ∆c i (not shown here) in a similar manner to Figure 5, and the model ∆H m0 is just as sizeable under the influence of high wind forcing. There is evidence of correlated bivariate uncertainty data in Figure 8b, and a combination of on-and off-ice wave features for the respective enhanced plots discussed in the previous paragraphs are apparent. Deciphering the physical processes is complicated; however, the bivariate uncertainty data along a transect illustrates how ∆c i and ∆H m0 are related; Figure 9b shows these results for a cross section oriented along the ice edge (the long axis of Quadrilateral 2) on 21 November 2018 18:00. The ∆c i hindcast experiment conducted in this study intentionally adopted the default IC3 source term parameters. As mentioned in Section 1, considerable progress in the WW3 wave-ice interaction parameterisations has been made owing to the Thomson et al. (2018) measurements. However, Ardhuin et al. (2018) also explain that observation of wave attenuation could 365 also be reproduced with many model forcing and parameter combinations, which as stated in the article is not unexpected because different wave-ice interaction processes are taking place along the wave propagation path. As discussed in Section 2.3, WW3 offers three physics-based s ice parameterisations: IC2, IC3, and IC5. They are based on different attenuation mechanisms and dispersion relations (see Appendix B), but their default ice rheological parameters (as given in the manual or the WW3 regression test cases) also vary as apparently the sea ice parameters are context-based according to the manual. The 370 three main parameters used to tune the IC2, IC3, and IC5 attenuation rate are as follows: eddy viscosity (m 2 s −1 ) ν, ice density (kg 2 m −3 ) ρ, and effective elastic shear modulus (Pa) G, although G is not used in IC2. The default ρ values are consistent for these s ice parameterisations. The default ν values are 153.6e−3, 1.0e+0, and 5.0e+7 for IC2, IC3, and IC5, respectively, and the default G values are 1.0e+3 and 4.9e+12 for IC3 and IC5, respectively. Rogers et al. (2016) adopted G = 0 assuming it is negligible in nilas and pancake ice fields. Based on the enormous range of tunable parameters and other factors described 375 in this paragraph, the s ice parameterisation is a notable uncertainty source. To evaluate the relative significance of ∆c i , we compared the results with the wave-ice interaction source term uncertainty (s ice uncertainty) for the modelling period between 5 and 25 November 2018. The s ice uncertainty wave hindcast experiment was conducted using IC2, IC3, and IC5 based on the same model setup described in Section 2.3 except all three simulations used BST-AMSR2 as SIC forcing.
A pair of ∆H m0 datasets were derived for the ∆c i and s ice uncertainty wave hindcast experiments. They were collated for Arctic Ocean, we defer to Ardhuin et al. (2018) to determine the test case and selected 50 cm. From a wave-ice modelling perspective, SIT effectively serves as a tuning parameter when forced as a homogeneous field. For example, the attenuation rate k i of IC2 as shown in Appendix B has SIT in the form of (1 + k r M ) in the denominator: M = ρhi ρw (Liu and Mollo-Christensen, 1988) where h i is the SIT, and ρ and ρ w are the ice and sea water density. If we take a deep water wavelength corresponding to 7 s wave period, changing the SIT from 10 cm to 50 cm increase the k i by at most 3 %. k i sensitivity on SIT 390 examined in Wang and Shen (2010); Mosig et al. (2015) (IC3 and IC5) appears more sensitive; as such, sensitivity analysis was conducted for our model. Repeating the ∆c i and s ice uncertainty experiments with 50 cm SIT, max(∆H m0 ) values increased respectively to 2.34 m and 1.95 m, but the ∆c i remains as the dominant error source. Further, IC3 was most affected by the SIT change for the equivalent transects of Figure 9 (not shown here). Even though there was no event during the study period when scattering were expected to be the dominant process (the implication of this is given in Section 6), sensitivity of the 395 finding to scattering was also examined by combining IS2 scattering with IC2 and IC3 with the default parameters. The results remained robust for the study period. Lastly, sensitivity to the choice of SIC forcing used in the s ice uncertainty experiment was examined by using ASI-6km instead of BST-AMSR2: the experiment also resulted in the same outcome. Table 2 is a list of SIC forcing used in the recent WW3 s ice developments as well as the forcing of the wave models analysed in Section 4. Some studies employed numerical sea ice models or considered various sources and assessed their 400 accuracy/suitability. It is interesting to learn that each wave-ice modelling study used different satellite retrieved SIC data.

Conclusions and discussions
The WW3 wave-ice models represent the exponential decay of waves in the presence of sea ice as dN dt = −2c i c g k i N (Equation 4). We investigated the effect of the satellite retrieved SIC uncertainty ∆c i on modelling waves in the refreezing Chukchi Sea MIZ using six SIC data sets based on the four commonly used retrieval algorithms: NASA-Team, Bootstrap, OSISAF, 405 and ARTIST-sea-ice. The wave hindcast experiment reveals ∆c i causes model wave height uncertainty ∆H m0 , and there is evidence that bivariate uncertainty data (∆H m0 and ∆c i ) are correlated, although off-ice wave growth is more complicated due to the cumulative effect of ∆c i along an MIZ fetch We compared the ∆H m0 distribution of the ∆c i experiment with that of the s ice uncertainty experiment. Both uncertainties are found to be considerable during the simulation period with maximum ∆H m0 values of 1.95 m and 1.44 m, respectively.

410
This result is found to be robust based on the sensitivity analyses that tested the SIT forcing and the inclusion of scattering.
Despite the s ice parameterisations being derived from different concepts and the WW3 wave-ice models completely missing the subgrid scale physics relating to sea ice field heterogeneity, the accuracy of satellite retrieved SIC used as model forcing is the primary error source of modelling MIZ waves in the refreezing ocean. The study outcome suggests wave-ice model tuning may not be as effective at this time when the knowledge of the true SIC field is too uncertain. It is worthy to note that swell 415 waves that propagate O(100) km into the ice-covered water where the scattering would likely be the dominant process were not observed during the study period. As such, the effect of ∆c i for such waves remains to be resolved.
Future improvements on the wave-ice models should come from two ends; continual developments of parameterised physics on the regional and pan-Arctic scale and working on a subgrid scale physical model on the other end. Solid and robust observational evidence through remote sensing and shipboard measurements is likely the key to connecting these two ends.

420
Appendix A: R/V Mirai measurement system R/V Mirai is equipped with two anemometers that were located on the foremast at 25 m elevation, and indicative wind conditions at the ship positions were derived from 10 minute vector moving averages of 6 s interval instantaneous true wind speed and direction. SST was measured −1 m below the sea surface with further 5 m inlet to the gauge while air temperatures were measured on the foremast at 23 m elevation.

425
Shipboard waves were obtained based on two methods: a microwave radar system (WM-2) (TSK Tsurumi Seiki Co., 2019) at the bow and stern of the ship, and nine-axis Inertial Moment Unit (IMU) (Piper-C#15), which is a device similar to the one used by Kohout et al. (2015). WM-2 has a sampling frequency of 2 Hz and collects raw sea surface elevation for 1,152 seconds at 35 minutes past each hour, and its integrated analogue system removes hull agitation and carries out Doppler correction. Bulk parameters like the significant wave height and period are produced based on the zero-crossing method. Wave 430 observations during the campaign from the WM-2 integrated analog system (TSK Tsurumi Seiki Co., 2019) were significantly affected by Doppler correction errors. Collins III et al. (2015) have shown shipboard measurements are less affected by this effect when ship speed is < 3 ms -1 . Applying a 2 ms -1 ship speed threshold greatly reduced conspicuously spurious data, and these data were used as indicative wave heights in this study (e.g., Figure 3). Piper-C#15 on board the vessel relies on an IMU. The processing method is consistent with Kohout et al. (2015) except 15 minute intervals were used instead of 1 hour.

435
Waves in the Chukchi Sea during the study period was dominated by wind seas, which have shorter wavelength relative to the ship dimensions. These waves are impeded by R/V Mirai's hull, so the shipboard Piper-C#15 has limitations on measuring wind seas. Response Amplitude Operator of R/V Mirai and the WM-2 data can be combined in theory to transfer IMU's high frequency signals to true surface elevations, but post-processing remains ongoing work. Although most of the Piper-C#15 data did not reflect the true wave field, the peak Piper-C#15 H m0 of 2.00 m during the on-ice wave event on 14 November 2018 440 agreed with the peak WM-2 H m0 ; this value is also comparable with the ERA5 H m0 as well. This provides confidence that the waves observed during this event was at least around 2.00 m.

Appendix B: Supplementary information of TodaiWW3-ArCS and the dispersion relation of WAVEWATCH III ® wave-ice interaction models
The TodaiWW3-ArCS geographical grid at high latitudes was based on the curvilinear grid implemented by Rogers and 445 Campbell (2009) with a polar stereographic projection of 75°N latitude (produced by Mathworks Matlab's polarstereo_inv function). The model domain coverage on the polar stereographic grids is as follows: the easting extent between −1,800 km and 1,512 km, and the northing extent between 520 km and 2,904 km. The geographical grid was defined using the International Bathymetry Chart of the Arctic Ocean bathymetry (Jakobsson et al., 2012) and the Global Self-consistent, Hierarchical, Highresolution Geography shoreline data (Wessel and Smith, 1996), and there are approximately 301,535 sea point cells. The 450 spectral grid was configured with 36 directional and 35 frequency bins with the latter ranging from 0.041 Hz to 1.052 Hz.
The WW3 wave-ice model dispersion relations for IC2, IC3, and IC5 are provided here. IC2 is based on the work of Liu and Mollo-Christensen (1988), and the dispersion relation is defined as follows according to the WW3 manual: h i is the SIT, ρ and ρ w are the ice and sea water density, ν is the viscosity, E is the Young's modulus of elasticity, φ is the Poisson's ratio, and P is the compressive stress in the ice pack.
G is the shear modulus.

465
IC5 is based on Mosig et al. (2015), and the dispersion relation is defined as follows according to the WW3 manual: ice observation is grouped in four phases as distinct ocean surface features were captured from four meteorological conditions the ship encountered. Figure A3 presents the shipboard measured wind and SST data as well as bilinearly interpolated (in space) ERA5 10 m wind speeds. During the first few days between 9 and 13 November 2018 (Phase 1), gradual sea ice 485 growth was observed both in extent and ice cake/floe sizes under generally calm surface conditions. On 14 November, the most significant on-ice wind event during the transect period occurred. The peak wind speed measured by R/V Mirai was 18 ms −1 , and > 10 ms −1 winds persisted for roughly 18 hours at the ship location. WM-2 and Piper-C#15 H m0 both peaked > 2.00 m in ice cover indicating energetic sea state of this event. The MIZ was mostly broken up ice fields on the following day, which was followed by the most apparent sea ice advance on 16 November as a seemingly dense ice field was encountered. We've  FIGURES.