Spectral attenuation of gravity wave and model calibration in pack ice

. We investigate an instance of wave propagation in the fall of 2015 in thin pack ice (<0.3 m) and use the resulting attenuation data to calibrate two viscoelastic wave-in-ice models that describe wave evolution. The study domain is 400 km by 300 km adjacent to a marginal ice zone (MIZ) in the Beaufort Sea. From Sentinel-1A synthetic aperture radar (SAR) imagery, the ice cover is divided into two regions delineated by the first appearance of leads. According to the quality of SAR retrievals, we focus on a range 20 of wavenumbers corresponding to 9~15 s waves from the open water dispersion relation. By pairing directional wave spectra from different locations, we obtain wavenumber-dependent attenuation rates, which slightly increase with increasing wavenumber before the first appearance of leads and become lower and more uniform against wavenumber in thicker ice after that. The results are used to calibrate two viscoelastic wave-in-ice models through optimization. For the Wang and Shen (2010) model, the calibrated equivalent 25 shear modulus and viscosity of the pack ice are roughly one order of magnitude greater than that in grease/pancake ice reported in Cheng et al. (2017). These parameters obtained for the extended Fox and Squire model are much larger than laboratory values, as found in Mosig et al. (2015) using data from the Antarctic MIZ. This study shows a promising way of using remote sensing data with large areal coverage to conduct model calibration for various types of ice cover. thinner ice and lower concentration than that above (after) the FAL. 90 Based on the distribution of SAR retrievals, we select the study domain by azimuth from 450 to 750 km and range from 0 to 400 km. Figures 1(b)(c) show the distributions of ice concentration from AMSR2 and ice thickness from SMOS (Soil Moisture and Ocean Salinity, https://icdc.cen.uni-hamburg.de/1/daten/cryo-sphere/l3c-smos-sit.html), respectively, in the region of interest. This work is to retrieve wave attenuation and use the result to calibrate viscoelastic models in thin pack ice (<0.3 m) dominant region. the gravitational acceleration, 𝑘 .E indicates the wavenumber in open water. Note the actual range of frequencies in the present case might be different due to the possible change of dispersion relation in pack ice. The dispersion relation cannot be measured from the instantaneous SAR data. Using accelerometers deployed on ice floes (Fox and Haskell, 2001) and from marine radar or buoy data (Collins et al., 2018), it was found that the wavenumber in ice cover region near 140 the open sea ( ≲ 10 km) was close to that in open water. more observations with higher quality data will improve modeling of the wave-ice interaction and the robustness of the models.


Introduction
The rapid reduction of Arctic ice in recent decades has become a focal point in climate change discussions.
The fact that this reduction exceeds model predictions (Stroeve et al., 2007;Comiso et al., 2008) emphasizes 35 the need to better understand the complex interaction between the sea ice, the ocean, and the atmosphere.
One of these interaction processes is between ocean waves and sea ice. Ocean waves help to shape the formation of new ice covers (Lange et al., 1989;Shen et al., 2001), break existing ice covers (Kohout et al., 2016), modify the upper ocean mixing (Smith et al., 2018), or potentially compress sea ice through wave radiation stress (Stopa et al., 2018a). In turn, ice covers alter the wave dispersion and attenuation, redistribute 40 wave energy, suppress wave-wind interaction and wave breaking (Squire, 2007(Squire, , 2018. To account for the ice effects in wave propagation, the spectral wave model WAVEWATCH III ® (WW3) has included several dispersion/dissipation (IC0, IC1, IC2, IC3, IC4, and IC5) and scattering (IS1 and IS2) parameterizations, called "switches" in WW3, to estimate the ice effect on waves (WAVEWATCH III ® Development Group, 2019). Both IC3 and IC5 theorize that sea ice can store and dissipate mechanical energy, 45 hence they model the ice cover as a viscoelastic material. Therefore, the energy dissipates and the dispersion changes when waves propagate through an ice field. The storage property is reflected in the potential and elastic energy, and the dissipative property is in the equivalent viscous damping. The difference between IC3 and IC5 is that IC3 is an extension of the viscous ice layer model with a finite thickness (Keller, 1998) by including elasticity into a complex viscosity (Wang and Shen, 2010), while IC5 is an extension of the thin 50 elastic plate model (Fox and Squire, 1994) introduced by Mosig et al. (2015) by adding viscosity into a complex shear modulus (equivalent to the complex viscosity via the Voigt model). Below we will refer to these two viscoelastic models as WS and FS respectively.
To use these models for wave forecasts in ice-covered seas, one needs to determine these equivalent elasticity and equivalent viscosity for all types of ice covers. To derive these parameters using first principles is 55 challenging, as demonstrated by de Carolis et al. (2005), who obtained the viscosity of grease ice using principles in fluid mechanics of a suspension. Alternatively, an inverse method has been adopted to parameterize IC3. Using in-situ data from the R/V Sikuliaq field experiment , the WS model was calibrated to match the observed wavenumber and attenuation (Cheng et al., 2017). This al. (2018b)). A large wave event during this time with wave heights exceeding 4 m in the captured region 70 provided quality wave data. Ardhuin et al. (2017) developed a method to invert from wave orbital motion to directional wave spectra from SAR images of waves in sea ice based on the velocity bunching mechanism. Stopa et al. (2018b) refined this methodology for the partially and fully ice-covered regions by substantially reducing data contaminated by ice features with similar length scale as the wavelength. Using SAR retrieved wave spectra, the authors found the significant wave height attenuated steeply prior to the first appearance of 75 ice leads (denoted as FAL hereafter), with milder attenuation rate after the FAL.
An overview of ice and wave condition with available observations in the studied region is presented in Figure 1. More details of the retrievals are given in Stopa et al. (2018b). Significant spatial variability is observed in both the wave and ice conditions. Figure 1(a) presents a subregion captured by the SAR images covering from a portion of Alaska to the largest azimuth position where waves are last observed visually in 80 the images. Colors indicate the dominant wavenumber distribution of the retrievals, and arrows indicate the dominant wave direction defined in section 2. Cell size is coarsened to 12.5 km × 12.5 km to enhance visualization. The ice edge is indicated by the contours of ice concentration (<0.4) from AMSR2 (Advanced Microwave Scanning Radiometer 2, http://doi.org/10.5067/AMSR2/A2_SI12_NRT). An in-situ buoy: AWAC-I (a subsurface Nortek Acoustic Wave and Current buoy, moored at 150 o W, 75 o N) is marked by a 85 black asterisk. Except for the in-situ observations from the Sikuliaq ship (green diamond) and several drifting buoys (blue dots) near the ice edge, ice morphology information including ice types and their partial concentrations are absent. Nevertheless, the FAL (black dots) presumably marks the separation between discrete floes and a "semi-continuous" ice cover with dispersed leads. The ice condition of a subregion below (before) the FAL was more complex with thinner ice and lower concentration than that above (after) the FAL. 90 Based on the distribution of SAR retrievals, we select the study domain by azimuth from 450 to 750 km and range from 0 to 400 km. Figures 1(b  This paper is organized as follows: site description and wave characteristics are depicted in section 2. We use the directional wave spectra retrieved from the Sentinel-1A SAR imagery to derive the wavenumberdependent attenuation in section 3. The methodology of model calibration and results are presented in section 105 4. The methodology of calibrating IC3 used in Cheng et al. (2017) based on wave energy spectrum ( , ) in terms of frequency needs to be modified for the present wavenumber-based spectra. Section 5 discusses the characteristics of wavelength and attenuation in pack ice, calibrated viscoelastic parameters, and wave attenuation modeling. The final conclusions are given in section 6.

Wave characteristics in ice 110
The retrieved wavenumber-direction spectrum is divided into multiple bins of wavenumber ) from 0.011 m -1 to 0.045 m -1 with an increment of 0.002 m -1 , and 360 bins in directions with bin width as 1 o . The main wave direction + , for each wavenumber ) is defined as the mean of a Gaussian function fitted to ( ) , ), with an example given in Figure S1 of the supplemental material. Furthermore, we define the dominant wavenumber of the wave spectra, ),-./01213 , as the one corresponding to the maximum of the directional 115 sum ∫ ( ) , ) . In Figure 1, ),-./01213 generally declines crossing the FAL towards the north. Before the FAL, ),-./01213 varies with ice concentration but is insensitive to ice thickness variation and wave directions. After the FAL, ),-./01213 decreases in the wave propagating direction associated with the increase of ice thickness, where the ice field is presumably a semi-continuous cover populated with leads. 120 Figures 2(a)(b) show two-dimensional histograms of the main wave direction for each wavenumber + ) before and after the FAL, respectively, where wave direction is defined in the meteorological convention.
Grayscale indicates the occurrence frequency of + ) in each ) bin. We observe a significant change of + , crossing the FAL: + ) before the FAL spreads from 160 o to 190 o , while + , is more tightly clustered from 180 o to 200 o after the FAL. The difference before and after the FAL is most significant for ) < 0.035 m -1 125 For the subsequent spectral analysis, we further restrict ) to 0.019 m -1 ≤ ) ≤ 2 / ; , where ; is the azimuth cutoff indicating the minimum resolvable wavelength from the SAR imagery (Stopa et al., 2015;Ardhuin et al., 2017). Below this wavelength, the patterns of ice-covered ocean surface roughness from SAR imagery are more related to ice features rather than waves (Stopa et al., 2018b). For ) < 0.019 m -1 , energy density < ) , + , = is small with high spatial variation ( Figure A1 in the Appendix), hence treated as noise 130 band and removed from further study. The corresponding frequency (wave period ) range is estimated as 0.067 to 0.11 Hz (9 to 15 s) using the 135 open water dispersion relation, i.e., (2 ) B = .E , where is the gravitational acceleration, .E indicates the wavenumber in open water. Note the actual range of frequencies in the present case might be different due to the possible change of dispersion relation in pack ice. The dispersion relation cannot be measured from the instantaneous SAR data. Using accelerometers deployed on ice floes (Fox and Haskell, 2001) and from marine radar or buoy data , it was found that the wavenumber in ice cover region near 140 the open sea (≲10 km) was close to that in open water. We define an apparent wave attenuation for each ) by assuming an exponential decay of the wave spectral densities from location to location : 145 is the average of the main wave directions at and B; and \] are the distance and direction from A to B in the longitude-latitude coordinates, respectively. A selected pair of A and B is named as a pair hereafter. To reduce the uncertainties of naturally present ice and wave variability, a set of quality control criteria are applied to a pair before further analysis: 150 1) Ignore substantially oblique waves. Difference between ̅ and the vector from location to location is restricted to M ̅ − \] M ≤15°.

2) Avoid strong spatial variations of ice condition between A and B. Distance between A and B is
restricted to ≤ 60 km. 3) As the wave state change rapidly across the FAL mentioned in section 2, no pair across the FAL is 155 selected, in order to compare the influence of the ice morphology on wave attenuation. 4) The Pearson correlation coefficient between energy spectra \ ( ) , ) and ] ( ) , ) is greater than 0.9. The similarity of wave spectra patterns indicates that the wave fields from A to B are dominated by the same wave system. 5) Exclude outliers of the spectral attenuation where wave energy of B is higher or close to that of A. 160 ( ) ) > 10 -6 m -1 is required. 6) A selected pair has at least 10 data points of ( ) ) to do calibration in section 4.
We obtain 2634 pairs (2194 pairs before the FAL, and 440 pairs after the FAL) through the above criteria to calculate the apparent wave attenuation ( ) ) by Eq. (1). The empirical thresholds of criteria 2) and 4) are determined based on the sensitivity of the number of selected pairs depending on these values. Figure  show two-dimensional histograms of against ) before and after the FAL, respectively, where the domain is divided into equal 30 bins in log scale from 10 -6 to 10-4 m -1 . Grayscales represent the occurrence of at each combined bin of and ) . Red curves indicate the most frequent occurrence of ( ) ) against ) . It is observed that more data are collected before the FAL, and a slightly increasing trend of ( ) ) versus ) before the FAL, while ( ) ) obtained after the FAL are mostly lower and independent of ) . 170 The range of this apparent spectral attenuation is in agreement with Stopa et al. (2018b), in which, the authors selected multiple tracks throughout the ice region and focused on the overall decay of the significant wave height over hundreds of kilometers. The reduction of attenuation crossing the FAL shown in Figure   Because of the large study domain and the apparent difference of ),-./01213 in the east-west direction before FAL as shown in Figure 1(a), it is worthwhile examining the regional variability. Figure  , and its color indicates the mean ice thickness between the two. In the middle column, the values corresponding to ) = 0.019, 0.025, 0.031 and 0.039 m -1 obtained from these pairs are separated into two 190 groups: before and after the FAL. In each subdomain delineated by the longitude and the FAL, the distribution of for each ) is presented by a violin plot, whose width indicates the probability density distribution of and a circle marker inside indicates the median. These violin plots show that before the FAL is larger than that after the FAL for all ) in all three subdomains. We note that the sample size after the FAL in the easternmost subdomain (146 o W, 147 o W) is the lowest, corresponding to the largest variability of the violin plots. In 195 the right column, the power spectral densities (PSDs) retrieved at locations of the ends of the selected pairs in the left column are given. Curves indicate the averaged PSD per 0.1 degree in latitudes, with latitude marked by colors. As defined earlier, ) associated with the peak of a PSD curve is ),-./01213 . Before the FAL the magnitude of PSD drops rapidly while ),-./01213 varies slightly as the increase of latitude. In contrast, the PSDs vary slowly after the FAL while ),-./01213 decreases quickly. Note at high latitudes, 200 PSD at ),-./01213 (red) is higher than that of the same wavenumber at low latitudes. We will revisit this phenomenon in the discussion section.

Wave attenuation due to ice effect 210
Following Eq. (1), the apparent attenuation obtained above is determined by the energy difference between two locations. This apparent attenuation is the result of multiple source terms, including the wind input 01 , the dissipation through wave breaking -L , the energy transfer due to nonlinear interactions among spectral components 1d , and the dissipation/scattering of wave energy due to ice cover 0;e . The radiative transfer equation for surface waves concerning all the above effects is 215 The other source terms 01 and -L are estimated using formulations from Snyder et al. (1981) and Komen et al. (1984). For 1d , we select the Discrete Interaction Approximation (DIA, Hasselmann et al., 1985a, b) to estimate its value. These formulations and associated coefficients are described in Cheng et al. (2017).
Wherever needed in these formulations, is approximated by the open water dispersion relation with the measured ) . Ice concentration is from AMSR2, ice thickness is from SMOS, and wind 235 data is from the Climate Forecast System Reanalysis (CFSR, Saha et al., 2010).
Ice-induced wave attenuation is known from dissipation of wave energy and scattering of waves (e.g., Wadhams et al., 1988;Squire et al., 1995;Montiel et al., 2018). Here we attribute the attenuation entirely to the dissipative process with the following arguments. Ardhuin et al. (2018) reported that in the studied region, ice floe scattering is a weaker effect on wave attenuation compared with other processes, including the 240 boundary layer beneath ice, inelastic flexing of ice cover, and wave-induced ice fracture. The inelastic flexing and ice fracture may be considered as part of the dissipative mechanism within the ice cover already included in the viscous coefficient, but scattering is a re-distribution of energy, which must be isolated from the apparent attenuation before using the data to calibrate the viscoelastic models. We estimate the scattering effect based on the study of Bennetts and Squire (2012) as follows. In that study, wave attenuation by floes, 245 cracks and pressure ridges were examined. In the absence of in-situ observations, we assume that few and small ridges are present in the studied ice cover (< 25 cm), hence the effect of ridges is negligible. In our case of long waves propagating through such thin ice cover, Bennetts and Squire (2012) where 0 / is the attenuation rate due to the ice cover. Figure 5 shows two-dimensional histograms of against ) before and after the FAL, as in Figure 3. Not surprisingly, since wind effect is low due to the low open water fraction in the study region, and the relatively short distances between the pairs for the nonlinear 260 transfer of energy to accumulate, we observe that 0 / is very close to the apparent attenuation , and will be further discussed in the discussion section.

Wave-in-ice model calibration 265
In this section, we present the calibration of the WS model and the FS model. In both models, ) , 0 are solved from the respective dispersion relations for given ice thickness and viscoelastic properties. An optimization procedure is introduced to select the viscoelastic parameters that minimize the difference between the model 0 and the 0 / obtained from the SAR over the entire range of ) .

Dispersion relation 270
The dispersion relations of the two models are written as For the WS model, In this study, we use = 1000 m for deep water, E23e) = 1025 kg/m 3 , 0;e = 922.5 kg/m 3 and = 0.3 for ice.

Calibration methodology
For each pair selected in section 3, we optimize , by fitting modeled 0 to the measured 0 / over the ) / domain. Note that 0 , ) in the models are solved simultaneously for an array of . We need to convert the 285 theoretical 0 -relation to 0 -) / . Specifically, for given and , we solve arrays of ) and 0 through Eq. No weighting factor could produce a fitting over the entire spectral attenuation curve. For the present study, the variation of spectral attenuation is small as shown in Figure 5 with no particular region of emphasis, we choose =¨∫ ( , ) ) which has a broad band around the peak energy of the wave field.
We choose a search domain (10 -7 Pa ≤ ≤ 10 10 Pa and 10 -4 m 2 /s ≤ ≤10 4 m 2 /s) for the WS model. This search domain covers all other reported viscoelastic values for ice covers (e.g., Newyear and Martin, 1999;300 Doble et al., 2015;Zhao and Shen, 2015;Rabault et al., 2017). For the FS model, it is known that very large , are needed to obtain the level of 0 observed (Mosig et al., 2015). We thus choose a large search domain 10 Pa ≤ ≤ 10 20 Pa and 10 ≤ ≤10 15 m 2 /s. This parameter range is far beyond the measured data from solid ice (Weeks and Assur, 1967). The global optimization procedure is performed by the genetic algorithm using function ga in MATLAB and Global Optimization Algorithm Toolbox (R2016a). 305

Results
For the WS model, the optimized , are gathered into a small cluster around ≅ 10 P" Pa and a large cluster around ≅ 10 ª Pa with large variation in . The existence of two separate clusters of calibrated results was also found in grease/pancake mixtures near the ice edge (Cheng et al., 2017). As multiple solutions could be obtained in solving optimal problem of the nonlinear system Eq. (5), constraints are applied to select 310 solutions that are physical. In the optimization, we reject the solutions around = 10 P" Pa, and solutions with near 10 4 m 2 /s that the residual from Eq. (6)  apply the bivariate Gaussian distribution to the rest of points to obtain the 90% probability range of ( , ) 315 before and after the FAL, respectively. Means and covariances of the related probability density functions are given in Table 1. Slight difference of distributions of , before and after the FAL is noticed. The mean of elasticity (viscosity) is slight lower (higher) before the FAL than that after the FAL. The results imply that the ice layer behaves more elastic in the inner ice field and more viscous towards the ice edge. It is worth noting that the ranges of both and are about one order of magnitude greater than those of the 320 grease/pancake ice (Cheng et al., 2017). As expected, the effective elasticity in the WS model for more solid ice is closer to the elastic modulus of sea ice. For the FS model, the calibrated ( , ) are clustered, where scatter plots of ( , ) are given in Figure S3 in the supplemental material. Hence all data points are used in the bivariate Gaussian fitting. The resulting mean 325 and covariance values are also given in Table 1. Scatter plots of ( , ) are given in Figure S3 in the supplemental material. Notice that the mean values of , from the FS model are extremely larger than the intrinsic values of ice. The distribution of calibrated values is further discussed in the discussion section.
Hereafter we only elaborate on results of the WS model. Figures 6(a) shows the overall comparison of 0 -) between the measurements (Gray) and the WS model (black). Both the median (solid curve) and 90% 330 boundaries (dash curves) are in good agreement. In Figure 6(a), we superimpose an empirical model from Meylan et al. (2014), which is also included in WW3 to account for the ice effect as switch IC4. By fitting the wave buoy data from the Antarctic MIZ obtained in 2012, Meylan et al. (2014)   respectively. Contours represent ice thickness distribution from SMOS. The FAL is marked as red dots. The 345 spatial domain is divided into 12.5 km × 12.5 km cells to enhance visualization. Cell color indicates the averaged value of log J½ and log J½ of all pairs with midpoints inside the cell. Ignoring the outliers, it shows that is generally larger after the FAL (thicker ice with leads) than before the FAL, while has an opposite trend. Note that the outliers could be from multiple sources, such as noise in the retrieved wave data, spatial variability of ice condition, as well as assumptions made in selecting pairs and calculating attenuation. The 350 counterpart of the FS model is given in Figure S5 in the supplemental material.

Discussion
In this section we discuss some key wave characteristics mentioned in the analysis, and the behaviors of calibrated model parameters. Some thoughts about wave-in-ice modelling are provided at the end.

Evolution of wave characteristics
The behaviors of ),-./01213 and PSD in Figures 1 and 4 may come from multiple mechanisms. The decrease 360 of ),-./01213 towards the interior ice is in agreement with a similar study in the MIZ (Shen et al., 2018), as well as the lengthening of dominant wave periods commonly reported in many field observations (e.g., Robin, 1963;Wadhams et al., 1988;Marko, 2003;Kohout et al. 2014;Collins et al., 2015). The phenomenon is commonly explained by the low pass filter mechanism in literature. That is, higher frequency waves are preferentially attenuated, thus the power spectrum peak shifts to longer waves. However, though this 365 mechanism is supported by the increasing trend of 0 with increasing wavenumber before the FAL, it is insufficient to explain the observations after the FAL, where the 0 -) curve is practically flat ( Figure 5) and the peak of PSD shifts toward lower ) as latitude increases (Figure 4). Two other mechanisms also may cause a downshift of the wave energy towards lower wave numbers. One is the change of dispersion relation, so that the wavelength grows deeper into the ice cover. The other is the nonlinear transfer between wave 370 components which move the energy in the main wave frequency into the low frequency components. How to quantify these mechanisms is still an open question with little data especially in pack ice, where the skills in observations are still developing. A critical information is the wave period in wave spectra, especially in the region after the FAL.
The attenuation rate obtained in this study (~10 -5 m 2 /s) against wavenumber shows a slight increasing trend 375 before the FAL, while nearly flat trend and lower after the FAL. Similar magnitude of attenuation is found in the MIZ in both Arctic and Antarctic as reported in previous studies in the same period range but larger for shorter periods (e.g., Meylan et al., 2014;Doble et al., 2015;Rogers et al., 2016;Cheng et al., 2017). The analysis of attenuation against wavenumber in section 3 is obtained for pairs of observations over relatively short distances (<60 km). In the Appendix, we present another method to determine the overall attenuation, 380 similar to the analysis of the attenuation of the significant wave height shown in Stopa et al. (2018b). By analyzing wavenumber by wavenumber over the entire space where the SAR data covers, the apparent attenuation coefficient is obtained by fitting hundreds of data over long distance, ignoring the higher variation of ice thickness (0.01~0.3m) and the shift of dominant wavenumber from the PSD curves. The results are shown in Figure A2. In some cases, we can even have negative values of this attenuation, meaning energy 385 increases as wave propagates. What we would like to emphasize by showing this result is that over a very long-distance ice condition can change significantly, in addition to nonlinear energy transfer and wind input/dissipation, thus long-distance spectral analysis may become meaningless.

model calibration
The covariance values of the calibrated viscoelastic parameters are greater than those found in grease/pancake 390 ice using buoy data (Cheng et al. 2017). In processing the SAR data, there are many challenges to be dealt with. It is extremely difficult to separate 1) sea ice variability, 2) wave height variability, and 3) instrument variability (speckle and in-coherent SAR noise). All of these influence the variations we see in the wave spectral data. The mask used in determining qualified data used in this study was described in Stopa et al. https://doi.org/10.5194/tc-2019-290 Preprint. (2018b). Still, the scatter of the spectral attenuation shown in Figure 3 is significant even after the filtering 395 described in section 3.1. This data scatter results in a large range of calibrated model parameters. We believe that the scatter of calibrated , in both WS and FS models could be narrowed down by reducing the uncertainty of measured attenuation data. We also believe that for the FS model the calibrated values and spread could be reduced by modifying the objective function with additional constraint on the wavenumber, which presently shows an extremely large range as shown in Figure S4. Regardless, the FS model needs 400 much larger , to produce attenuation rates comparable to the observed data. This is the case not only for the pack ice but also for the MIZ. When analyzing data from the Antarctic MIZ reported in Kohout and Williams (2013), and by further constraining ) = .E /1.7 , Mosig et al. (2015) obtained = 4.9 × 10 JB Pa, = 5 × 10 Ê m 2 /s. Though far below the present case, these values are still orders of magnitude above the intrinsic values measured from solid ice (Weeks and Assure, 1967). Meylan et al. (2018)  . The leading term of which in small ) yields ≈ ( ) Pª ).

Thoughts of modelling wave-ice interaction
Damping models play a crucial role in spectral attenuation. At present, to use any specific model to describe wave attenuation is tentative. Identified damping mechanisms are many. For instance, boundary layer under the ice cover (Liu and Mollo-Christensen, 1988;Smith and Thomson, 2019), spilling of water over the ice 415 cover and interactions between floes (Bennetts and Williams, 2015), jet formation between colliding floes (Rabault, 2019), have all been reported in laboratory or field studies. A full waves-in-ice model that clarifies all important mechanisms is not yet available. While theoretical improvements are needed to better model wave propagation through various types of ice covers, practical applications cannot wait. Calibrated models that are capable of reproducing key observations must be developed in parallel to model improvements. The 420 present study provides a viable way to calibrate two such models available in WW3. These models lump all dissipative mechanisms in the ice cover into a viscous term. This type of model calibration studies has two obvious utilities. One, with proper calibration, models can capture the attenuation of the most energetic part of the wave spectrum. Two, the discrepancies may be used to motivate model development that includes missing mechanisms, thereby help future model development. Because different physical processes may play 425 different roles under different ice morphology, collecting wave data to calibrate these models under various ice types is necessary. Finally, more observations with higher quality data will improve modeling of the wave-ice interaction and the robustness of the models.

Conclusions
The wave spectra were retrieved from SAR imagery in the young pack ice predominant area in Beaufort Sea 430 in 12 October 2015, further north of the grease/pancake ice region along the ship and buoy tracks. According to the analysis of data retrieved over several hundred kilometers, the observed decrease of wave energy and lengthening of dominant waves towards the interior ice are consistent with earlier in-situ observations. We investigate wave attenuation of dominant spectral densities per wavenumber between two arbitrary locations in the region. Similar attenuation rates are observed for all wavenumbers from 0.019 to 0.045 m -1 (estimated 435 wave period from 9 to 15 s). We then use the wave attenuation rates to calibrate two viscoelastic-type wavein-ice models. Both models can generally match the observed attenuation corresponding to the energetic portion of the wave spectra, with large difference on dispersion. For the WS model, the calibrated shear modulus (viscous parameter) in the region beyond the first appearance of leads with thicker ice is slightly larger (smaller) than that of the region closer to ice edge before the first appearance of leads. The resulting 440 wavenumber is within 5% of that produced from the open water dispersion.
Wave-ice interaction is complicated due to many co-existing physical processes. At present no models have fully integrated all identified processes. However, marine operations in the ice-covered seas need wave forecasts. The present study demonstrates a method based on measured data to calibrate existing models, so that they can be applied to meet the operational needs. As the models improve, further calibration exercises 445 may be performed accordingly. For example, the eddy viscosity model (Liu and Mollo-Christensen, 1988) attributed the wave attenuation entirely to the water body under the ice cover. This model is also included in WW3, with the eddy viscosity as a tuning parameter. The present analysis can also be used to calibrate that model. From a physical point of view, dissipative mechanisms may be present simultaneously inside the ice cover and the water body underneath. However, calibration of complex models with multiple co-existing 450 processes is a difficult task which requires much more dedicated study.
We conclude by noting that high resolution spatial data from remote sensing provide new opportunities to investigate the wave-ice interaction over a large distance and different ice types. However, details of ice types and temporal observations are in development. To reach a full understanding and thus a complete waves-inice model requires collaboration from observation and modeling efforts. 455

Appendix
We investigate the decay of the dominant energy component < ) , + , = over a long distance along selected tracks with fixed longitude (145 o W, 146 o W, …, 150 o W). Figure A1 shows < ) , + , = (blue dots) collected within a given longitude interval 150 ± 0.1 o W starting from 74.5 o N towards the north. ( ) ) for ) < 0.019 m -1 is too scattered to show any attenuation trend. While as ) increases, the data become tighter with a clear 460 decay trend. To obtain attenuation coefficient Ñ( ) ), we fit the < ) , + , = by an exponential curve (black line) in each panel based on an exponential wave decay assumption.

465
The resulting Ñ( ) ) against ) is shown in Figure A2, as well as another five curves associated with different longitudes, processed in the same way. Figure A2 shows an increasing trend of Ñ( ) ) against ) as expected, and Ñ( ) <0.019 m -1 ) are mostly negative. While this plain method blinds the variation of ice morphology and peak shifting of the power spectra density, thus it is insufficient to understand the damping mechanism in water-ice interaction. Hence, as our interest is to determine the attenuation in more consistent ice conditions, 470 shorter distance between measuring points is adopted in the main body.

Code and data availability
The code and data necessary to reproduce the results presented in this paper can be obtained from the 475 corresponding author.