Retrieval of Ocean Wave Heights from Spaceborne SAR over the Arctic Marginal Ice Zone with a Neural Network

,

mode data acquired in the Arctic were collocated with data from four radar altimeters (RA), 23 yielding 126,128 collocated data pairs. These data were separated into training and testing 24 datasets to develop a BPNN model for retrieving SWH. Comparing the S1 retrieved SWH 25 using the testing dataset with the RA SWH yielded a bias of 0.17 m, a root-mean-square error 26 of 0.71 m and a scatter index of 23.05% for SWH less than 10 m. The S1 retrieved SWH 27 were further compared with CFOSAT/SWIM data acquired in the Arctic between August 28 2019 and May 2020 to validate the SWIM performance on wave measurements at different 29 beams. 30

Plain language summary 31
The rapid decline of sea ice in the Arctic creates wider marginal ice zone (a transit from open  32 water to sea ice, MIZ) than ever. Some studies have suggested that interaction between ocean 33 dynamics (e.g., sea surface wind and wave) and sea ice is one possible feedback to retreat of 34 sea ice in the Arctic. Therefore, ocean wave data in the MIZ is highly desirable. Synthetic 35 aperture radar, as an active remote sensing technique, can operate independent on sunlight 36 and weather conditions and image the earth with high spatial resolution. However, due to the 37 complicated imaging process of ocean waves by spaceborne SAR, retrieval of sea state 38 parameters by SAR data has been investigated for decades. Here, we developed an algorithm 39 based on a back propagation neural network to retrieve significant wave height from 40 spaceborne SAR data. This provides a chance of obtaining wave height information in both 41 large coverages and high spatial resolution from satellite observation, and therefore, can 42 contribute to scientific study, offshore operation and shipping in the Arctic. 43

Introduction 44
Prior to the launch of the Chinese French Oceanic Satellite (CFOSAT) with its onboard 45 Surface Waves Investigation and Monitoring (SWIM) sensor, the only sensor capable of 46 imaging ocean waves in two dimensions from space was the spaceborne synthetic aperture 47 radar (SAR), which provides images with high spatial resolution. The SAR imaging 48 mechanism of ocean waves is complex which is generally explained by three modulations: 49 tilt modulation, hydrodynamic modulation and velocity bunching (Valenzuela, 1978;Alpers 50 et al., 1981). While tilt and hydrodynamic modulations are also shared by real-aperture radar 51 as the dominant imaging mechanisms of ocean waves, velocity bunching is unique for SAR 52 to image ocean waves. The moving scatterer of water particles with a velocity either towards 53 or away from a moving SAR sensor, causes an azimuthal shift in SAR images. In addition, 54 velocity bunching in the SAR resolution cell leads to an azimuth cut-off, that is, the minimum 55 SAR-detectable wavelength of ocean waves traveling in the azimuth direction. Therefore, the 56 nonlinearity of SAR ocean wave imaging complicates their retrieval. In the following, we 57 briefly summarize the existing methods used to retrieve ocean wave information in terms of 58 both two-dimensional spectrum and integral wave parameters. 59 The Max Planck Institute (MPI) scheme developed by (Hasselmann & Hasselmann,60 1991; Hasselmann et al., 1996) is the widely used method to retrieve two-dimensional ocean 61 wave spectra from spaceborne SAR data. The MPI method iteratively searches for the 62 minimum of cost function to retrieve wave spectra from SAR by using a numerical ocean 63 wave model (e.g., the WAM model) for the first-guess wave spectra. These first-guess wave 64 spectra provide the wave propagation direction and compensate for the loss of wave 65 information in high-frequency during the SAR imaging process. By this way, nonlinear 66 retrievals can get the complete two-dimensional spectra of ocean waves. Therefore, these 67 methods strongly depend on the first-guess wave spectra as prior information. Alternatively, 68 wind vectors measured by a scatterometer can be utilized to estimate generally missed 69 windsea information by SAR imaging ocean waves, e.g., the semi parametric retrieval 70 algorithm scheme (SPRA) developed by Mastenbroek & De Valk (2000), which also applies 71 full nonlinear mapping relations between ocean waves and SAR imaging. The SPRA 72 combines the observed SAR spectrum with collocated scatterometer wind vectors to estimate 73 the windsea spectrum, while the residual signal in the SAR spectrum is considered as the 74 swell. The SAR image spectra employed by the abovementioned methods are derived from 75 intensity image. Alternatively, the partition rescaling and shift algorithm (PARSA) developed 76 by Schulz-Stellenfleth et al. (2005) inputs the cross spectra derived from single-look-complex 77 SAR data to a nonlinear inversion. This type of nonlinear retrieval methods can generally 78 yield two-dimensional ocean wave spectra, enabling the derivation of integral ocean wave 79 parameters, e.g., the significant wave height (SWH) and mean wave period. Nevertheless, 80 due to their dependency on prior information, these methods inconvenient for wide 81 applications as ocean wave model spectra are generally not publicly available. Moreover, 82 nonlinear retrievals can be degraded to quasi-linear retrievals. By inputting the cross 83 spectrum which resolves the ambiguity of ocean wave propagation in the quasi-linear 84 retrieval (Engen & Johnsen, 1995), one can generally obtain ocean swell spectrum. The 85 advantage of this approach is that prior information is no longer needed. Even though 86 quasi-linear retrievals cannot yield full two-dimensional ocean wave spectra, the obtained 87 swell spectra are particularly important for studying swell propagation and decay (Li, Compared with theoretical-based methods, these empirical type methods can yield full sea  97  state (both wind wave and swell wave) parameters without conducting complicated nonlinear  98 retrievals. Therefore, first-guess spectral information is no longer needed. Furthermore, 99 although they are called empirical algorithms, the input parameters are not chosen randomly. 100 For instance, the SAR normalized radar cross section (NRCS) and the corresponding image 101 spectra, which are often used in empirical algorithms, form the basis of traditional nonlinear 102 SAR ocean wave retrievals. Moreover, CWAVE-type algorithms attempt to incorporate the 103 nonlinear relations among SAR images and ocean wave parameters using 2nd-order 104 polynomials by cross-multiplying the input parameters. However, the nonlinear relationships 105 between SAR image and ocean wave parameters are often too complex to be sufficiently 106 represented by a 2nd-order polynomial. Therefore, the backpropagation neural network 107 (BPNN), which has the ability to fit nonlinear relationships, has been employed to retrieve 108 the SWH from SAR images and to improve the retrieval accuracy. 109 BPNN, a traditional machine learning method proposed in the 1980s (Rumelhart et al., 110 1986), has been shown to be effective at fitting nonlinear problems between input and output 111 parameters. BPNN considers an iteration as the combination between the forward 112 transmission of information and the backward transmission of error. The network is trained 113 iteratively until the global error satisfies the preset accuracy or until the number of training 114 iterations exceeds the specified maximum number of learning iterations. BPNN consists of an 115 input layer, one or multiple hidden layers and an output layer. The input and output layers 116 comprise the input and output data of the model, respectively. The hidden layer, which is not 117 visible to users, is the key to fitting the relationship between the input and the output data. 118 Stopa  RAs and the coarse resolution of CFOSAT/SWIM, limited ocean wave products are available 143 in the Arctic MIZ with both high spatial resolution and large coverage for studying the 144 interaction between ocean waves and sea ice. The Copernicus Sentinel-1A (S1A) and 145 Sentinel-1B (S1B) satellites have been in orbit since April 2014 and April 2016, respectively. 146 This constellation significantly reduces the revisit period, thereby yielding a high temporal 147 resolution, particularly in the polar regions. Approximately 3,000 S1 images are acquired 148 every month in the Arctic, and most of the Arctic can be covered within two days. The twins 149 have extensively acquired data in extra-wide (EW) swath mode and interferometric wide (IW) 150 swath mode in the Arctic. Fig. 1 presents an example of the spatial coverages of the EW data 151 acquired by S1A and S1B within six days in 2019. Additionally, the EW and IW data 152 acquired in the Arctic are generally in polarization combination of co-polarization and 153 cross-polarization, dedicated for sea ice monitoring (e.g. the motivation of this study is to develop an algorithm dedicated for retrieving SWH in the 157 Arctic using S1 data in HH polarization. These data certainly are useful for studying the 158 interaction between sea ice and ocean dynamics, as both sea ice and marine-meteo parameters 159 can be derived from SAR simultaneously. Following the introduction, the datasets used in this study are introduced in Section 2. 170 Section 3 presents the methodology, including the data collocation and the development of 171 the BPNN model to retrieve SWH by S1 EW data in HH polarization. Verification of the 172 BPNN model for retrievals are shown in Section 4. In Section 5, a detailed comparison 173 between the S1 retrieved SWH and the collocated CFOSAT/SWIM data is presented. A 174 summary and the conclusions are given in the last section. 175 2 Datasets 176 2.1 S1A and S1B EW data 177 Most S1A and S1B EW data acquired in the Arctic are in dual-polarization (HH and 178 HV). In this study, S1 Level-1 ground range detected (GRD) data in HH polarization are used 179 to retrieve SWH. S1 EW images have a swath width of 400 km with a spatial resolution of 40 180 m. The radar incidence angle of the EW data ranges from 18.9° in the near range to 47.0° in 181 the far range. Radiometric calibration and thermal noise removal of the EW data are 182 conducted according to the S1 user manual (ESA, 2016). The NRCS 0 is obtained by: 183 where is the digital number read from the tiff data file, is the noise vector, and is 184 the calibration factor. The noise vector and calibration factor are given in the product noise 185 and calibration metadata. 186 The EW GRD data used herein span the period between January 2017 and October 2019, 187 comprising approximately 113,500 images. 188 As most of the S1 EW and IW data acquired over in situ buoys are in VV polarization, 189 we found only 305 pairs of S1 data and National Data Buoy Center (NDBC) buoy data in the 190 period from October 2014 to October 2019 . Therefore, in this study, we used 191 RA measurements of SWH in the Arctic as ground truth to develop the BPNN model. 192 2.2 RA SWH data 193 The RA SWH data are from four missions: CryoSat-2, Jason-2, Jason-3 and SARAL. 194 These RA-measured SWH are screened, and only good quality data are retained. The 195 CryoSat-2 data are provided by the European Space Agency (ESA,196 http://science-pds.cryosat.esa.int/) and can reach latitudes of 88°N. We used pole-to-pole 197 Level-2 CryoSat-2 data with a 1 Hz sampling frequency and extracted the data with values 198 for 'surf_type' of 0 (ocean) and 'flag_instr_op_mode' of 1 (good quality). Jason-2, Jason-3 199 and SARAL are all provided by the European Organization for the Exploitation of 200 Meteorological Satellites (EUMETSAT, https://archive.eumetsat.int/usc/). While the Jason-2 201 and Jason-3 missions can reach latitudes of only 66.15°, SARAL can cover more of the Arctic, 202 up to 81.49°N. We extracted the data of these three RA missions with values for 203 'surface_type' of 0 (ocean) and 'qual_swh' of 0 (good quality). 204 Prior to using the RA data from the four missions above to construct the BPNN model, 205 we conducted cross-comparisons among the four RA missions. The RA missions in each pair 206 were matched with temporal interval less than 1 hour and spatial distance less than 10 km for 207 cross-comparisons in the region above 60°N. The corresponding statistical parameters of 208 these comparisons are listed in Table 1. The twin satellites, Jason-2 and Jason-3, achieve the 209 best agreement with a bias of 0.04 m and a root-mean-square error (RMSE) of 0.00 m. The 210 comparisons between CryoSat-2 and Jason-2/3 also show good compatibility with biases of 211 -0.01/-0.02 m and RMSEs of 0.02/0.01 m. The differences between CryoSat-2 and SARAL 212 and between Jason-3 and SARAL are slightly higher with biases of -0.06 m and 0.06 m, 213 respectively, but the RMSEs of these two comparisons are only 0.01 m and 0.02 m, 214 respectively. Therefore, the discrepancies among the SWH data from these four RA missions 215 are minor, and we did not calibrate the data based on data from a single RA mission. 216 217 3.1 Collocation of the S1 EW and RA data 245 The S1 EW scenes were collocated with the RA data with a temporal window of less 246 than 90 minutes. A total of 4,273 S1 EW scenes were collocated with the data from the four 247 RA missions, among which 1,834 and 2,439 scenes were acquired by S1A and S1B, 248 respectively. The spatial distributions of the collocated S1A and S1B EW images are 249 presented in Fig. 3. Then, the S1 EW sub-images with dimensions of 256 × 256 pixels (i.e., 250 10,240 m×10,240 m) collocated with the RA footprints were collected as matchups. Finally, a 251 total of 153,485 collocation data pairs of S1 and RA data between January 2017 and October 252 2019 were obtained. 253 (a) (b) Figure 3. Spatial distributions of (a) S1A and (b) S1B images collocated with the RA data in 254 the period between January 2017 and October 2019. 255

256
As a large amount of S1 EW data were acquired in the Arctic MIZ, they often present a 257 mixture of sea ice and open water. Therefore, we used the reanalysis daily sea ice cover 258 product (with a grid size of 1 km) of the ice mapping system (IMS) to filter out ice-covered 259 sub-images. 260 In addition, the quality of the S1 sub-images has a significant impact on the SWH 261 retrievals. We used the homogeneity parameter (Schulz-Stellenfleth & Lehner, 2004) to filter 262 out S1 sub-images on presenting some oceanic and atmospheric features not related to ocean 263 surface waves. On the other hand, the IMS data are daily products and have discrepancies 264 with the S1 observations, which are snapshots. Therefore, a homogeneity test can also discard 265 sub-images presenting sea ice features (particularly pancake and icebergs (Lehner & 266 Ocampo-Torres, 2003)) not identified by the IMS data. The homogeneity parameter is 267 defined in (2): 268 where Φ is the power spectral density of each sub-image. Generally, the sea surface is 269 considered homogeneous for <1.05. The statistics of the homogeneity values are shown in 270 Fig. 4. 271  Therefore, we also chose parameters similar to those used in CWAVE-type algorithms to 289 retrieve SWH by the S1 data: the mean NRCS (denoted ̅ 0 ), normalized image variance 290 (cvar), and 20 spectral parameters computed from the variance spectrum of a sub-image. The 291 ̅ 0 and cvar are computed as follows: 292 where 〈 〉 is the mean intensity of an S1 sub-image. 293 The 20 spectral parameters are extracted from the SAR image spectrum using a set of 294 orthonormal functions. The SAR image spectrum is estimated by computing the image 295 periodogram with a two-dimensional fast Fourier transform (FFT) algorithm. These 296 orthonormal functions can extract the features of the image spectrum from 20 different 297 directions. The method for extracting the 20 SAR image spectral parameters is described in 298 detail in the Appendix. 299 The previously developed CWAVE-type algorithms for SAR WV data do not include the 300 parameter of incidence angle, as WV data have fixed incidence angles of approximately 23° 301 for ERS/SAR and ENVISAT/ASAR WM data or angles of 23° and 33° for S1 WV data. 302 However, S1 EW mode data have incidence angles ranging from 19° to 47°, while the NRCS 303 significantly varies with the incidence angle. Therefore, the incidence angle should be 304 included as a key input parameter to the neural network. Previous studies on developing 305 empirical methods for SWH retrieval by SAR data used different forms of incidence angles, 306 such as tan respectively. 320 Figure 5. Structure of the proposed BPNN model for retrieving SWH from S1 data. 322 The function of each node in the network is to calculate the scalar product of the input 323 vector and weight vector using a nonlinear transfer function. This nonlinear transfer 324 function, called the activation function, is the key to improving the approximation ability of a 325 neural network and is expressed as follows: 326 where the activation value of node j is , is the connection weight vector from the 327 nodes of the upper layer to node j of this layer, represents the bias of node j, is the 328 output of node j, and (•) is the activation function of a node. The activation function of the 329 second hidden layer is a sigmoid function (we used logsig), and the activation function of the 330 other hidden layers is the hyperbolic tangent function (tansig); these two functions are given 331 in (8) and (9), respectively. The activation function of the output layer is "purelin", a linear 332 transfer function. Fig. 6 illustrates these activation functions, in which the x-axis and y-axis 333 are the input and output of the nodes, respectively, and the solid line represents their 334 relationship. 335 respectively. 338 After forward-propagating the data in the input layer to the hidden layers, the network 339 computes the result in the output layer. A global error is computed based on the 340 performance function of the mean square error (MSE), which is given as follows: 341 where is the connection weight from the hidden node j to the output node , 342 is the bias of the output node , (•) is the activation function of the output layer node, and 343 is the true value of the training data.
where represents either the input or the output parameters, and are the 366 minimum and maximum values of each parameter, respectively, and represents the 367 normalized input and output data. After normalization, the input and output parameters are 368 between 0 and 1. To use the proposed BPNN model to retrieve SWH from S1 EW data, the 369 output data should be anti-normalized to practical values. 370 Three parameters are assigned as termination conditions. The maximum number of 371 iterations is set to 5,000, and the minimum of MSE is set to 0.001. The maximum failure time 372 is set to 6, where failure is defined when the global error in the current iteration is larger than 373 that in the previous iteration. 374 After training the BPNN model, three statistical parameters, namely, the bias, RMSE 375 and SI, are used to evaluate the comparisons between the S1 retrieved SWH using BPNN and 376 the RA SWH. The three parameters are computed as follows: where is the S1 retrieved SWH and is the RA SWH. 378 Fig. 7 (a) and (b) show comparisons between the S1 retrieved SWH and RA SWH using 379 the training and testing datasets, respectively. With respect to the comparison using the 380 training dataset, the bias of 0.02 m, the RMSE of 0.62 m and the SI of 20.67% show that the 381 S1 retrieved SWH is close to the RA SWH. The comparison using the testing dataset achieves 382 almost identical statistical parameters with a bias of 0.02 m, an RMSE of 0.63 m and an SI of 383 21.07%. This finding indicates that the trained BPNN model has stable performance on the 384 SWH retrieval from S1 EW data. However, these comparisons suggest that the retrieved 385 SWH is lower than the RA SWH when the SWH exceeds 4 m, as indicated by the error bars 386 in Fig. 7. Moreover, the underestimation increases with SWH increasing. 387 (a) (b) Figure 7. Comparisons between the S1 retrieved SWH and RA SWH using (a) the training 388 dataset and the (b) testing dataset. 389 390 A method of duplicating training data in high sea state is used to solve the 391 underestimation afflicting the S1 retrieved SWH. Fig. 8 (a) shows a histogram of the RA 392 SWH in the training dataset suggesting that the amount of data in high sea state is far less 393 than the amount of data in low to moderate sea state, e.g., between 2 and 4 m. This is likely a 394 major cause of the underestimation of S1 retrieved SWH in high sea state. To solve this 395 problem, we arbitrarily changed the distribution of the training dataset to normal distribution 396 (as shown in Fig. 8 (b)) by discarding some training samples with SWH lower than 3.3 m and 397 duplicating samples with SWH higher than 3.3 m, resulting in another training dataset with 398 153,691 data pairs. We retained the original testing data to verify the training of the network, 399 which histogram is shown in Fig. 8 (c). The BPNN was re-trained by using the adjusted training dataset and the retrievals using 404 the new network were compared with the RA data, as shown in Fig. 9. Fig. 9 (a) shows the 405 comparison using the full training dataset (including the duplicated cases in high sea state), 406 and Fig. 9 (b) presents the comparison without including the duplicated cases. Both (a) and (b) 407 suggest that the underestimation of SAR retrievals is effectively resolved using the adjusted 408 training dataset. By comparing Fig. 9 (a) with Fig. 9 (b), one can refer to the effect of those 409 duplicated cases in high sea state in the BPNN training. By excluding the duplicated data 410 from the comparison, all three parameters increase accordingly. The comparison based on the 411 training dataset without duplicating ( Fig. 9 (b)) reveals statistical parameters that are almost 412 identical to those of the comparison using the testing dataset with a bias of 0.17 m, an RMSE 413 of 0.71 m and an SI of 23.05%, as shown in Fig. 9 (c). However, these statistics parameters 414 are higher than those achieved using the original training dataset (Fig. 7 (b)). Therefore, we 415 resolved the underestimation of SAR retrievals from moderate to high sea states but at the 416 cost of increasing the overall statistical parameters. 417 418 (a) (b) (c) Figure 9. Comparisons between the S1 retrieved SWH and RA SWH using (a) the training 419 dataset after duplicating, (b) the original training dataset and (c) the testing dataset. 420

421
In the following, we present two cases to demonstrate the advantages of sea state 422 observations by S1 in the Arctic MIZ based on the proposed method. The first case is in the 423 east of Greenland; the S1A EW data were acquired at 19:15 UTC on 6 December 2018. The 424 retrieved SWH using the developed BPNN model is shown in Fig. 10 (a). Fig. 10 (c) presents 425 the corresponding ERA-5 reanalysis wind field at 19:00 UTC on 6 December 2018, showing 426 a cyclone weather situation with wind speeds above 20 m/s in the northwest of the S1 SWH 427 map leading the SWH to exceed 6 m therein. The overlaid track in Fig. 10 (a) is the 428 collocated CryoSat-2 SWH measurements from 18:45 to 18:46 UTC. The collocated S1 429 retrievals (triangles) with the Cryosat-2 SWH (circles) along the track are shown in Fig. 10  430 (e). From this scatter diagram, the S1 SWH is close to the CryoSat-2 SWH, especially 431 between the latitude of 61°N to 62.5°N, where the difference between the S1 SWH and 432 CryoSat-2 SWH is only 0.10 m. The S1 retrievals are slightly lower than the CryoSat-2 SWH 433 south of 61°N (lower by 0.93 m on average) but are higher than the CryoSat-2 SWH north of 434 62.5°N (where the sea state is generally above 7 m) with significant spatial variation. 435 The second case is also in the east of Greenland but the data were acquired by S1B at 436 08:13 UTC on 28 November 2018. The S1 retrieved SWH is shown in Fig. 10 (b), in which 437 the gray area represents the coverage of sea ice (extracted from the IMS data). Fig. 10 (d)  438 presents the ERA-5 reanalysis sea surface wind field at 08:00 UTC on 28 November 2018. 439 The case shows a strong wind above 15 m/s blowing from the northeast to the southwest, and 440 as a result, the SWH increases from northeast to southwest. The overlaid track represents the 441 measurements of Jason-3 from 09:12 to 09:13 UTC on 28 November 2018, which is 442 approximately 1 hour later than the S1B sensing time. The collocated S1 retrievals (triangles) 443 with the Jason-3 SWH (circles) along the track are shown in Fig. 10 (f). In according with the 444 Jason-3 SWH, the S1 SWH decreases with the increasing of latitude, as shown both in the 445 scatter diagram (Fig. 10 (f)) and in the SWH map ( Fig. 10 (d)). The S1 retrievals are slightly 446 higher than the Jason-3 SWH south of 65 .35°N and between 65.55°N and 65.6°N Figure 10.

Comparison between the S1 retrieved SWH and CFOSAT/SWIM data 464
In this section, we compared the S1 retrieved SWH with the collocated CFOSAT/SWIM 465 data acquired between August 2019 and May 2020. The SWIM measurements at nadir beam 466 and 10° beam were used for a comparison with the S1 retrieved SWH. The S1 retrieved SWH 467 and SWIM SWH were matched up in a temporal interval of 90 minutes. We first extracted the 468 collocated S1 sub-images with dimensions of 70 km × 90 km, which is the same area as the 469 wave cell of the 10° SWIM beam, and then we retrieved the SWH of these sub-images at a 470 2.56 km spatial resolution, finally, we computed the mean values of the retrieved SWH. 471 Finally, we obtained 32,403 collocations between the data from S1 and the SWIM nadir beam 472 and 1,283 collocations between the data from S1 and the SWIM 10° beam. 473 Fig. 11(a) shows the comparison between the S1 retrieved SWH and the collocated 474 SWIM SWH at nadir. The bias is 0.15 m, the RMSE is 0.60 m and the SI is 18.98%, which 475 are similar to the comparison with the RA SWH, as presented in Fig. 9 (c). Fig. 11 (a) also 476 reveals that the S1 retrieved SWH is close to the SWIM nadir SWH when SWH is lower than 477 approximately 5 m, but when SWH is above 6 m, the former is lower than the latter by 478 approximately 0.84 m. The comparisons presented by Hauser et al. (2020) also suggested that 479 the SWIM nadir SWH is at least 0.5 m higher than the ECMWF model SWH when the SWH 480 is above 6 m. Fig.11 (b) Figure 11. Comparisons between the S1 retrieved SWH and SWIM SWH at (a) nadir beam 495 and at the 10° beam on the (b) left track and (c) right track. 496 Due to the nonlinear imaging mechanisms of ocean waves by spaceborne SAR, the 497 retrievals of two-dimensional wave spectra and sea state parameters may suffer problems for 498 short waves or azimuthal-traveling waves. The two-dimensional ocean wave information 499 available from the SWIM sensor provides a unique opportunity to verify whether the SAR 500 retrievals of sea state parameters depend on the wavelength and wave direction. Fig. 12 (a) 501 and (b) show the differences between the S1 retrievals and SWIM SWH at the 10° beam 502 (with a total of 1,283 data pairs by combining the collocations with the left and right SWIM 503 tracks) varying with the dominant wavelength and azimuth wave direction (i.e., the dominant 504 wave direction relative to the S1 azimuth direction) of SWIM. The step sizes of the 505 histograms are 40 m for the dominant wavelength and 10° for the azimuth wave direction, 506 respectively. The overlaid error bars represent the mean absolute bias and its standard 507 deviation. 508 Fig. 12 (a) indicates that the S1 collocations with the 10° SWIM beam data are 509 concentrated on the sea states with dominant wavelength less than 300 m. In this range, the 510 mean absolute SWH bias is less than 0.70 m with limited fluctuations, and moreover, it is not 511 found that the bias increases for the retrievals of waves with a relatively short wavelength, 512 e.g., less than 100 m. With an increasing dominant wavelength, both the bias and the standard 513 deviation increase, which is slightly different from our expectation, as SAR is generally 514 considered suitable for the retrievals of ocean waves with long wavelength. However, the 515 large bias (>0.5 m) obtained for data with a long wavelength (>300 m) may also attribute to 516 quite less amount of collocated data pairs, accounting for only 5.79% of the total number of 517 collocated data pairs. 518 Interestingly, the collocations between the S1 and 10° SWIM beam data are 519 concentrated mainly on the sea states with azimuthal wave direction between 0° and 45° and 520 between 135° and 180°, namely, close to the SAR flight direction. The biases for the 521 collocated data in these two wave directions ranges are generally between 0.50 m and 0.75 m, 522 and moreover, they are quite stable with no dependence on wave traveling directions. For the 523 collocation data pairs with azimuthal wave traveling direction between approximately 60° 524 and 120°, the biases are relatively large, generally larger than 0.75 m, and the fluctuations are 525 quite distinct. These large biases may also be attributed to the smaller amount of collocated 526 data with azimuthal wave directions in this range, accounting for only 4.89% of the total 527 number of collocated data pairs. 528 These two comparisons suggest that the S1 retrieved SWH based on the proposed BPNN 529 model tends to be independent on the wavelength and azimuth wave direction, while more 530 collocations need to be collected in the future to draw a more reliable conclusion. We further presented a case to compare the S1 retrieved SWH with the SWIM data in 536 the Arctic MIZ. Three consecutive S1 EW images were acquired over the east of Greenland 537 from 18:01 to 18:03 UTC on 26 February 2020. Fig. 13 (a) shows the S1 retrieved SWH of 538 this case, in which the gray area represents sea-ice covered area based on the IMS data. The 539 overlaid circles represent the SWIM nadir SWH observations, while the squares to the left 540 and right of the track of circles are the SWH by SWIM at the 10° beam. The black arrows on 541 the squares reflect the dominant wave direction derived from the SWIM data. The SWIM 542 data were acquired at 18:38 UTC, 37 minutes later than the S1 acquisitions. From the 543 northwest to the southeast, the SWH shows a trend of increasing and then decreasing, 544 reaching a peak value of approximately 6 m at 66°N. Fig. 13 (b) is the sea surface wind field 545 provided by the scatterometer onboard CFOSAT obtained at the same acquisition time as the 546 SWIM data presented in Fig. 13 (a). The wind speed rises from 10 m/s to nearly 20 m/s and 547 blows to the southwest, then decreases to 7-8 m/s with wind direction turning from southwest 548 to southeast. 549 550 (a) (b) Figure 13 (a) S1 retrieved SWH map (background) with collocated SWIM SWH at nadir 551 (circles) and at the left and right 10° beam tracks (squares). The arrows on the squares 552 indicate the dominant wave directions. (b) Corresponding sea surface wind field by the 553 scatterometer onboard CFOSAT. The ID of the three SAR images are 554 S1A_EW_GRDM_1SDH_20200226T180128_20200226T180232_031427_039E37_EE97, 555 S1A_EW_GRDM_1SDH_20200226T180232_20200226T180332_031427_039E37_559A, 556 and 557 S1A_EW_GRDM_1SDH_20200226T180332_20200226T180432_031427_039E37_19C3. 558 SWH at nadir, and the solid line represents the collocated S1 retrieved SWH at the SWIM 565 nadir track. The red and blue circle symbols represent the SWIM SWH on the right and left 566 tracks of the 10° beam, respectively. The diamonds with the same color are the collocated S1 567 retrieved SWH. In the region between approximately 64°N and 70°N, the SWH varies from 3 568 m to 6 m, and the SWIM SWH at nadir shows a similar spatial variation as the S1 retrievals. 569 However, the SWIM SWH is systematically higher than the collocated S1 SWH by 570 approximately 0.57 m. 571 The SWIM SWH of the left and right tracks at the 10° beam are also higher than the 572 collocated S1 retrieved SWH by 0.80 m and 0.55 m on average, respectively. Fig. 14  To further investigate this case, we chose three two-dimensional wave spectra provided 580 by SWIM at the 10° beam for demonstration, as shown in Fig. 15. Their integral wave 581 parameters and the collocated S1 retrieved SWH are listed in Table 3. Fig. 15 (a) shows the 582 sea state involving both wind sea (with a peak wavelength of approximately 152 m) and swell 583 (455 m); consequently, the dominant wavelength in region 1 is 145 m. The SWH for this 584 region by S1 and SWIM are similar with values of 3.11 m and 3.56 m, respectively. In region 585 2 ( Fig. 15 (b)), the sea state developed further, with longer dominant wavelength of 184 m 586 compared with the sea state at region 1. With sea state increasing, the difference between the 587 S1 retrieval and SWIM SWH increases to approximately 1.0 m (4.69 m vs. 5.63 m). The 588 two-dimensional wave spectrum presented in Fig. 15 (c) suggests that the sea state in region 3 589 is swell-dominated with a wavelength of 283 m. This swell system should have developed 590 from windsea at previous times, as its wave direction (60.68°, going to) is approximately 45° 591 from the local wind direction (15°). In addition to this dominant swell peak, there is another 592 weak peak with a wavelength of approximately 200 m and a wave direction of approximately 593 15°, consistent with the sea surface wind direction, which may indicate a young swell just 594 leaving the generation area. For region 3, although its wave height is lower than that of region 595 2, the SWIM SWH is still higher than the S1 retrieved SWH by 0.6 m. 596 This case reveals complicated sea state conditions with a mixture of windsea and swell 597 systems. Swells developed in previous times propagated further, and they coexisted with 598 locally generated windsea or young swells, as the high wind field also continuously moved. 599 The SWIM nadir SWH shows better agreement with the S1 retrievals than the data at the 10° 600 beam in this case. Although the SWIM SWH at both the nadir beam and the 10° beam are 601 higher than the S1 retrievals, the differences between the SWIM nadir and S1 retrieved SWH 602 seem to be systematic, while the differences between the 10° beam and S1 retrieved SWH are 603 significantly variable. 604 605 Table 3. Parameters and collocated S1 SWH of the three wave spectra in Fig. 15 integral wave parameters are listed in Table 3. All the wave spectra are oriented with respect 610 to true north (up represents north). 611 612

Summary and conclusions 613
The interaction between ocean dynamics and sea ice in the Arctic starts to draw more 614 attention due to the rapid decrease in the sea ice extent. As the basis of ocean dynamics, 615 accurate measurements of ocean wave parameters by remote sensing data in the MIZ are 616 highly desirable. S1A and S1B have extensively acquired spaceborne SAR data at both high 617 spatial resolution and large spatial coverage in the Arctic, providing unique advantages in the 618 acquisition of sea state information in the Arctic MIZ. Therefore, in this study, we focus on 619 developing a practical method to derive sea state parameter of the SWH from S1 SAR data, 620 which can be used further to study the interaction between ocean waves and sea ice. 621 Previous studies have demonstrated that empirical algorithms are practical to derive 622 integral wave parameters, e.g., the SWH and mean wave period, from SAR data than 623 traditional nonlinear inversions, as these algorithms do not need a priori information. In this 624 study, we adopted the idea of previous SAR-ocean wave empirical algorithms, but 625 incorporated these ideas into a BPNN model. BPNN  team, which indicates that the S1 retrievals should be of relatively good quality. However, a 650 comparison of the same dataset of S1 retrievals with the SWIM SWH at the 10° beam (on 651 either the left track or the right track) shows that the SWIM SWH is much higher than S1 652 retrievals with a bias of approximately 0.4 m and an RMSE of 0.90 m. Moreover, both the 653 statistical analysis and the case study indicate that the differences between the 10° SWIM 654 beam SWH and the S1 retrievals vary considerably. Although the difference between the S1 655 retrieved SWH and the SWIM SWH at the 10° beam is rather large, the S1 retrieved SWH is 656 independent of the dominant wavelength and azimuthal wave direction, indicating that the 657 proposed BPNN model can yield stable retrievals of SWH by S1 data. 658 These comparisons suggest that the quality of SWIM wave data should be improved in 659 the future. In October 2020, the SWIM development team announced that the current 660 modulation transfer function (MTF) has been adjusted and the reprocessing of all SWIM data 661 since the beginning of the mission will be triggered. We expect better SWIM data for further 662 research. 663 On the other hand, the S1 retrievals based on the proposed BPNN also have room for 664 improvement. One issue that remains unresolved is that it is difficult to retrieve correctly 665 SWH less than 1.5 m due to the insensitivity of SAR signals to low sea states. SAR 666 cross-polarization data are less saturated with high wind compared with co-polarization data 667 (Monaldo, et al., 2017). We recently developed a robust method for denoising S1 668 HV-polarized data , therefore, we expect to obtain better results for high sea 669 states by combining data in both HH and HV polarization. Furthermore, to date, only integral 670 wave parameter of SWH has been retrieved based on a neural network; hence, it might be 671 possible to retrieve two-dimensional wave spectra based on deep learning methods without 672 through the complicated nonlinear inversions. 673 The MATLAB code to retrieve SWH by the S1 data in HH polarization using the 674 proposed method was published in Zenodo for public sharing (Wu, 2020). 675 Acknowledgments: The S1 SAR data are downloaded from the Copernicus data hub 676 (https://scihub.copernicus.eu/), the CryoSat-2 data are downloaded from the CryoSat-2 677 Science Server (https://science-pds.cryosat.esa.int/), the Jason-2/3 data and SARAL data are 678 downloaded from EUMETSAT (https://archive.eumetsat.int/), and the CFOSAT SWIM data 679 are downloaded from AVISO (https://aviso-data-center.cnes.fr). Use of the reference data of 680 IMS data (https://www.natice.noaa.gov/ims/), the ERA-5 data 681 (https://cds.climate.copernicus.eu/cdsapp#!/home) and the CFOSAT scatterometer data 682 (https://osdds.nsoas.org.cn) is also acknowledged. This work was supported in part by the 683 National Key Research and Development Project (2018YFC1407100) China. 684

Appendix: Estimation of the SAR Image Spectrum 685
The SAR image spectrum is estimated by computing the image periodogram with a 2-D 686 FFT algorithm. The idea is to divide an image with 256 × 256 samples into 2 × 2 687 sub-images with 128 × 128 samples and then to compute the FFT of each sub-image and 688 obtain the power spectral density. Finally, the SAR image spectrum is acquired by computing 689 the average of 4 power spectral densities. 690 The 2-D FFT is applied to every normalized sub-image G: 691 where 128 represents the size of every sub-image. The power density spectrum is 692 denoted by : 693 = ( ) 2 #( 2) Then, summing the four power density spectra and averaging them, the entire SAR 694 image spectrum P is given by: 695 The power density spectrum needs to be normalized to ensure that the integral of the 696 image in the frequency domain is equal to that in the spatial domain. The normalized image 697 spectrum is denoted as ̅ : respectively. 724