Spatiotemporal observations of wave dispersion within sea ice using Sentinel1 SAR TOPS mode

In this paper, we derive the dispersion relation for waves within ice ‐ covered waters from the Sentinel ‐ 1 Interferometric Wide swath mode by using an innovative new implementation of the method proposed by Johnsen and Collard (2009). We apply this method to the spatial burst overlap area that is present due to the terrain observations with progressive scans (TOPS) technique. The advantage of this implementation is the use of a larger time separation between subsequent images that results in a less noisy and higher ‐ quality imaginary spectrum, which allows us to obtain spatiotemporal information of the dispersion relation. We studied seven wave events in the Barents Sea, where we have accompanying in situ data of sea currents and sea ice draft, available from the Barents Sea Metocean and Ice Network measurement program. Having simultaneous data on the currents is a signi ﬁ cant improvement over previous studies and allows us to quantify these data's in ﬂ uence on the dispersion relation. The derived dispersion relation from Sentinel ‐ 1 is derived for long waves propagating through icy waters (peak wavelengths 100 – 350 m), and it does not deviate from the theoretical open ‐ water dispersion relation. At present, however, the spatial resolution of synthetic aperture radar (SAR) data is too coarse and does not enable the study of the dispersion relation for short waves within sea ice. method. The new implementation bene ﬁ ts from the steering of the satellite antenna in a certain mode to allocate overlapping images of the ocean surface, which have not been previously used. Further, our study area is located in the Barents Sea, where we have data on, among others, the ocean currents. Having these data is an improvement in and of itself, as many previous studies lacked data that were collected on site. Our results show that we can study the wave dispersion relationship for long waves within the sea ice and, for some cases, in open water. At present, short waves cannot be studied with satellite imagery because the amount of details the satellite can capture is insuf ﬁ cient. modify the open ‐ water wave dispersion considerably. To demonstrate this, we computed the MAPE for a theoretical case using the maximum observed current speed of 1.06 m/s. The MAPE is then up to 11% higher than the case with the actual observed current speeds, which clearly shows it is vital to have accompanying data on the currents.


Introduction
The rapid change in sea ice extent (e.g., Comiso et al., 2017) and sea ice thickness (e.g., Renner et al., 2014) in the Arctic region provides a clear sign of climate change. The shrinkage in sea ice cover has led to longer open-water seasons, increased wave intensity , new shipping routes (Melia et al., 2016), and the growth of offshore operations. To reduce the risks and hazards of coastal and offshore operations, it is of the utmost importance to accurately forecast the wave climate. However, in the polar regions, wave prediction is problematic due to the presence of sea ice. Not only do the waves influence the ice cover through, for instance, ice breakup (Collins et al., 2015;Kohout et al., 2016) and by compacting ice edges (Stopa, Sutherland, et al., 2018), but the ice also impacts waves by refraction, shoaling, reflection, scattering, dissipation, and changing the dispersion relation (e.g., Squire, 2007;Squire et al., 1995). Therefore, waves propagating through icy waters show amplitude attenuation and have a different propagation speeds compared to waves traveling in open water. Both the change in the wave speed and the attenuation of the wave amplitude can be described by a dispersion relation, for which various mathematical formulae exist throughout the literature (see section 2). It is common to express the wave number in these formulae using complex notation, which allows the description of the amplitude attenuation (the imaginary part) together with changes in the wavelength (the real part). In this paper, we obtain the real part of such dispersion relations from Sentinel-1 (S1) Interferometric Wide (IW) swath mode in the marginal ice zone (MIZ), where waves have been observed by synthetic aperture radar (SAR) hundreds of kilometers into the ice pack (e.g., Ardhuin et al., 2015;Carsey et al., 1989).
A major challenge in validating the theoretical formulae of the wave dispersion in icy waters and for determining the model parameters is the need for data. The data can be collected from either in situ measurements or remote sensing. The complication of retrieving observations of the dispersion relation is that both the temporal and spatial measurements are simultaneously required. In situ measurements from buoys provide time series, that is, in the temporal domain, while observations from remote sensing provide information in the spatial domain. Furthermore, the data collected from in situ measurements are sparse, are typically available for a limited amount of time, and have proved to be challenging and expensive to collect in the harsh environment of the polar regions. With its large and continuous day and night coverage, in addition to its outstanding performance in all weather conditions, SAR onboard the S1 constellation is a valuable source of wave information in ice-covered oceans (e.g., Ardhuin et al., 2017;. The study by Collins et al. (2017) gives an overview of available observations (visual, in situ, remote sensing, and laboratory) of wave dispersion in the polar regions. From these observations, both lengthening and shortening of waves in pancake and brash ice are observed. All previous studies on remote sensing observations of wave dispersion use airborne or satellite SAR imagery to investigate the change in the twodirectional spectrum for waves entering icy waters (Liu et al., 1991;Shuchman et al., 1994;Wadhams & Holt, 1991). As noted in these studies, the analysis of SAR imagery to retrieve the wave number spectrum is nontrivial. A modulation transfer function (MTF) is required, which in open water consists of tilt modulation, hydrodynamic modulation, and velocity bunching (e.g., Alpers et al., 1981;Hasselmann et al., 1985). The method presented in this study does not require the use of the MTF because we assume a linear ocean-to-SAR transform and we directly retrieve the dispersion relation from the SAR image cross-spectrum (see section 4.5). Furthermore, in our analysis, we include the effects of sea currents on the wave dispersion relation by using in situ measurements. This addition is a significant improvement over previous studies since a major drawback of all the previous studies on the change in wave number in the MIZ is the lack of simultaneous current measurements (Collins et al., 2017).
The Arctic Sea State field experiment is a measurement campaign that was recently conducted in the western Arctic region . Wave buoys and shipboard X-band marine radar were used to derive the dispersion relation for waves in icy waters. From these observations, it was concluded that the dispersion of long waves in sea ice does not deviate significantly from the theoretical open-water dispersion relation (Cheng et al., 2017, their Figure 2). Collins et al. (2018) performed the most comprehensive study on the dispersion relation for waves in sea ice to date by using the Arctic Sea State field experiment data. For wave frequencies lower than 0.3 Hz, the measurements were close to the theoretical open-water dispersion relation. For wave frequencies in the range of 0.3 to 0.5 Hz, a deviation between measurements and the theoretical open-water dispersion relation was found. Collins et al. (2018) claimed that this deviation is consistent with the mass-loading model (see section 2) but did not follow the contour lines exactly.
While the latest study by Collins et al. (2018) presented detailed observations of the dispersion relation from in situ measurements, the aim of this study is to present and validate an approach to retrieve the wave dispersion relation from satellite SAR imagery. We use an innovative implementation of the method proposed by Johnsen and Collard (2009) on the spatial burst overlap area of S1 IW data. We analyze seven S1 IW images taken over the Barents Sea during the period of the Barents Sea Metocean and Ice network (BaSMIN) measurement program. This measurement campaign, together with the studied wave events, is introduced in section 3. Our method to derive the wave dispersion relation is described in section 4. Section 5 presents the results and is divided into three parts: (1) a comparison between the derived wave dispersion relation in open water and the one obtained in ice-covered seas, (2) the validation of the dispersion relation for relatively long waves, and (3) the derived dispersion relation for the shortest observed wavelengths within the sea ice. Finally, the results are summarized in section 6 before we draw the main conclusions.

Wave Dispersion
Wave dispersion (or frequency dispersion) refers to the difference in the wave speeds of wave components with different frequencies: Waves with different wavelengths propagate at dissimilar speeds and therefore disperse (Snodgrass et al., 1966). The derivation of the theoretical open-water wave dispersion relation is based on the assumptions of irrotational fluid motion, inviscid water, and incompressibility. It relates the spatial and temporal wave kinematics and is defined by where ω and k denote the angular frequency and the wave number expressed in a fixed frame of reference, respectively, g is the gravitational acceleration, d is the water depth, U is the current velocity, and α is the angle between the direction of the mean current and the wave propagation. When waves propagate through ice-covered oceans, equation (1) will be altered in such a way that causes changes in the wavelength, the phase and group velocities, the wave direction, and the wave height while keeping the frequency and action unchanged (Collins et al., 2017). Greenhill (1887) was the first to propose a modification to the theoretical open-water dispersion relation to make it applicable to ice-covered waters; over the past few decades, substantial effort has been put into the development of theoretical models to derive dispersion relations for waves in icy waters. These models include the mass-loading model (Weitz & Keller, 1950), the viscous layer model (Keller, 1998), and the viscoelastic model (Wang & Shen, 2010). All these models are based on different assumptions of the ice and water layers. Therefore, it is unlikely that one formulation can be used for all types of ice. Following the notation of Meylan et al. (2018), the general dispersion relation for waves traveling in deep water in the presence of sea ice and the absence of sea currents may be written as where the factor Q(k, ω) is the modification of the theoretical open-water dispersion relation. If no ice is present, Q(k, ω) is zero, and the theoretical open-water dispersion relation is restored. As an example, the simplest dispersion relationship for waves traveling in ice waters is the mass-loading model, for which the factor becomes where ρ ice and ρ w are the ice and the water densities, respectively, c is the ice concentration, and h is the ice thickness. For other wave dispersion formulae and the corresponding Q(k, ω), see, for instance, Meylan et al. (2018).
One useful application of wave dispersion relations valid for icy waters is the estimation of the sea ice thickness. A method was first developed by Wadhams and Holt (1991), who used the mass-loading model to estimate the ice thickness in the MIZ. By observing the ice concentration, wave number and angular frequency, the ice thickness can be computed from equations (2) and (3). Wadhams et al. (2002) extensively tested this method and obtained the wave number in sea ice, k i , by deriving the directional wave number spectrum in ice with the Hasselmann inversion scheme (Hasselmann & Hasselmann, 1991) and the cross-spectral inversion scheme (Engen & Johnsen, 1995). The wave number in open water, k ow , was acquired from the directional wave number spectrum in open water, and by using the theoretical open-water dispersion relation (ω 2 = gk ow ), the angular frequency, ω, was found. Wadhams et al. (2002) found that the derived ice thickness was overestimated. According to Collins et al. (2017), a possible explanation is that the spatial resolution of the SAR is too poor, so a shift of one wave number bin already results in a large change in the wavelength. In later studies by Wadhams et al. (2004) and Wadhams et al. (2018) , the viscous layer model by Keller (1998) was used, and good agreement with in situ measurements of pancake ice was found.

Site Description, BaSMIN Measurement Program, and Studied Wave Events
To assess the results from our method to retrieve the wave dispersion relation from satellite SAR imagery, we used in situ data collected during the BaSMIN measurement program. This 3-year metocean and ice measurement campaign took place from October 2015 until October 2018. It is a joint industry project led by Equinor, and the measurement campaign was conducted by Fugro GEOS Ltd. The data were kindly made available for academic research, and this is the first time the data are being published. During this field campaign, five wave-scan moorings and five ice profilers were deployed, and their locations are indicated in Figure 1. The wave-scan moorings measured the wave parameters, near-surface currents, seawater temperature, seawater salinity, and meteorological parameters. The ice profilers measured the ice draft, seawater temperature, ice drift speed, and current velocity.
The bathymetry of the Barents Sea is shown in Figure 1c. The Ice Moring Stations IC3, IC4, and IC5 are located at water depths of 157, 139, and 173 m, respectively. Therefore, when using the deep-water limit (H > 0.5 · L, where H is the water depth and L is the wavelength), waves with wavelength smaller than approximately 300 m can be considered deep-water waves. The Ice Mooring Stations IC1 and IC2 are placed at water depths of 52 and 68 m, respectively. For most wave events, this depth range can be classified as intermediate water. For this reason, we only looked at satellite passes close to IC3-IC5, as they are in deep water.
In Figure 1a, the ice thickness is provided for the date 4 April 2017 (Tian-Kunze et al., 2018) as taken from the Soil Moisture Ocean Salinity Earth Explorer mission (Kaleschke et al., 2012). The ice thickness in Figure 1a is between 20 and 40 cm around the ice mooring stations, which is consistent with the measured ice draft at locations IC3-IC5 but is slightly underestimated at IC1 and IC2 (not shown). The ice concentration ( Figure 1b) is from Advanced Microwave Scanning Radiometer 2 (Spreen et al., 2008) for the same date. From both the Soil Moisture Ocean Salinity ice thickness and ASMR2 ice concentration, there is an ice tongue extending south, encompassing Ice Mooring Stations IC1 and IC2. This ice tongue follows the Spitsbergen Bank because of the limited water depth, where ice forms more rapidly than in deeper water, and because of the sea ice drift. The sea ice drift has a directional trend toward the southwest, and its magnitude is on average larger for IC1 and IC2 than for IC3-IC5, which is in agreement with the general surface circulation in the Barents Sea (e.g., Oziel et al., 2016).
The relevant data collected during the measurement campaign for this study are the current speed and the ice draft. The current speed is necessary to calculate the wave dispersion (see equation (1)). At IC1 and IC2, they generally vary between 0.05 and 0.3 m/s and, at IC3-IC5, between 0.2 and 0.8 m/s. The maximum current speed observed during the entire campaign was 1.06 m/s, observed at IC1. For the 3 years of the measurement campaign, ice was mostly recorded from March to May. The ice was typically thin with a mean ice draft varying between 0.1 and 0.5 m. According to the sea ice nomenclature (World Metrological Organization, 1970), which relates the ice thickness to ice categories, young ice and first year thin ice are found. Furthermore, by studying some Sentinel-2 optical images, we found that the ice close to the ice edge is usually inhomogeneous, with ice floes and some open-water patches present.
Seven wave events were selected to test and validate our proposed method to derive the wave dispersion relation from S1. An overview is provided in Table 1, and these events were selected, so they cover a wide range of peak wavelengths. For the events with long waves (#1-4), the dispersion relation in sea ice is expected to be identical to the open-water dispersion relation (Collins et al., 2018), and hence, these events are used to validate our algorithm. The Events 5-7 have peak wavelengths in the range of 55-70 m and are in the remainder of this paper referred to as short waves. These waves were the shortest observed peak wavelength where the waves were still visible within the sea ice.

Method Description
The method described here results in spatiotemporal observations of the wave dispersion, that is, simultaneous measurements of the wave number and the angular frequency. Our method is based on the S1 WV Level-2 algorithm proposed by Johnsen and Collard (2009), which is innovatively adapted to S1 TOPS mode Single Look Complex (SLC) data. The traditional processing technique to estimate the cross-spectra (Engen & Johnsen, 1995) from strip-map data is performed on single-look data. This methodology usually solves the 180°wave ambiguity and removes the speckle noise bias. However, it results in a noisy and poor-quality imaginary cross-spectrum when applied to TOPS SLC data. Therefore, we propose to apply the cross-spectra estimation method proposed by Johnsen and Collard (2009) on the burst overlap area, which leads to a higher-quality imaginary cross-spectrum. This method can be used to more reliably determine the wave direction.
The burst overlap area is present because the S1 constellation uses the Terrain Observation with Progressive Scans SAR (TOPSAR) technique (Zan & Guarnieri, 2006) in its IW and Extra Wide swath modes. In addition to beam steering across the range direction, as with the traditional scanning SAR acquisition, the antenna beam is also steered in the azimuth direction, which results in a single-look area and a burst overlap area for each subswath. The scanning geometry is sketched in Figure 2, and the single-look and burst overlap areas are indicated.
The well-established method by Engen and Johnsen (1995) is commonly used to produce multiple independent looks from the SLC image of the single-look area and to consequently compute the cross-spectrum. The different sublooks that are produced have a maximum separation time of approximately half the integration time (≈0.4 s) . Applying this technique to TOPS SLC data results in a poor-quality imaginary part of the cross-spectrum because the progressive sweep of the antenna results in a dwell time of only 0.15 s. Moreover, the phase change (introduced by wave propagation) between the resulting subsequent looks is small and, therefore, noisy. In this study, we overcome this challenge by using two intensity images from the burst overlap area, in which we have two looks of the same area that are separated in time by approximately 1.8-2.1 s, depending on the subswath and the range location. The phase change between subsequent images is now significantly larger and less noisy. Therefore, the quality of the computed imaginary part of the crossspectrum via our method is much better than that calculated using the Note. The wave and wind parameters are taken from wavescan mooring WS3. Hs = significant wave height; L ws = wavelength; W ad = wave direction "coming from," relative to north and positive clockwise; W is = wind speed; W id = wind direction "coming from," relative to north and positive clockwise; h ice = daily mean ice thickness computed as the average from Ice Mooring Stations IC3-IC5; U = current speed averaged over Ice Mooring Stations IC3-IC5.
Figure 2. Acquisition geometry of the S1 TOPSAR technique (adapted from Park et al., 2018). There are three subswaths for the Interferometric Wide swath mode, and the burst cycle time is approximately 2.6 s. During each burst, the antenna beam is steered forward in the azimuth direction (orange dotted line to the purple dotted line) for each subswath. Overlap areas occur between the subswaths and between the bursts for each subswath. The overlap area between the bursts of a subswath, used in the processing, is indicated in the figure. standard processing technique, which is necessary for the accurate estimation of the wave dispersion relation. In addition, the wave direction can be resolved more accurately than the standard processing of the single-look area.
In Figure 3, we present a flow chart that shows the steps of the proposed method. In the next few sections, the different steps are explained in detail. Finally, we present an example case to illustrate the results.

SAR Input Data and Preprocessing
We use data acquired from the S1 constellation operating in IW swath mode. S1 carries a C-band radar instrument and has an incidence angle that ranges from 29°to 46°. Multiple polarization options are available, but throughout this study, we use the HH polarization. The total swath width is approximately 250 km, and it consists of three subswaths. We use the raw Level 0 data that can freely be downloaded from the European Space Agency's Sentinel Data Hub. These data are processed to Level 1 SLC data product with the same method that the European Space Agency (2013) uses but without merging the bursts in the azimuth direction and the subswaths in the range direction. This preprocessing step results in standard SLC images with a maximum spatial overlap between the bursts and swaths.

Step 1: Find the Extent of the Sub-images
After the preprocessing of the data, the next step is to define multiple subimages with the spatial overlap area of consecutive bursts. We first find the area where the two bursts overlap at the azimuth time. However, we cannot use the entire overlap area because the targets located at the beginning and at the end of each burst are not illuminated by the entire antenna pattern. Therefore, we find the length of the subscene in the azimuth direction by plotting the normalized intensity averaged over the range for both overlapping bursts. The length of the sub-image in the azimuth direction (y direction) is then found using the following criterion: b I nþ1 y ð Þ− b I n y ð Þ <0:5; where b I nþ1 y ð Þ and b I n y ð Þ are the normalized intensities averaged over the range for burst n + 1 and burst n, respectively. The limit of 0.5 is a compromise between using as many points in the azimuth direction as possible while excluding the part where the intensities fade out at the beginning and at the end of each burst. For all the images, this step results in approximately 320 points in the azimuth direction. In the range direction, we used 1,000 points for each sub-image. The resolution of the SLC data in the azimuth and range directions is approximately 14 and 5 m, respectively, resulting in a total sub-image area of approximately 4.5 × 5 km. The number of bursts used for each wave event differs. Bursts are manually selected for processing when waves are visible. Therefore, for events with long waves, more bursts are used as the waves propagate further into the ice pack. For each burst, we used multiple sub-images in the range direction, with 50% overlap, which typically results in 30 sub-images for each burst. A sub-image is manually processed when ocean waves are visible within the sea ice.

Step 2: Compute the Cross-Spectra
The image cross-spectrum between two overlapping images from two different bursts is computed following the method proposed by Johnsen and Collard (2009). The cross-spectrum (P (n,n+1) ) is defined by where I (n) and I (n+1) are the two overlapping intensity sub-images of burst n and burst n + 1, respectively. The two images are separated by the look separation time τ, and k ¼ k x ; k y À Á is the wave number domain (k x in the range direction and k y in the azimuth direction). The asterisk (*) denotes the complex conjugate, and < > represents ensemble averaging over the sub-images. Next, the inverse of the crossspectrum is computed, which is referred to as the cross-variance function, ρ (n,n+1) and is given by where x is the target pixel (with range and azimuth values). The cross-variance function is smoothed with a Hanning window, h, raised to the power of β. Throughout this study, we used a value of 2 for β. The final cross-spectrum is computed as the Fourier transform of the smoothed cross-variance function according to Finally, wavelengths larger than 500 m are filtered out by masking the middle of the cross-spectrum.

4.4.
Step 3: Convert to a 1-D Spectra To obtain the 1-D cross-spectra, we first compute the real part, P real k ð Þ; and imaginary part, P imag k ð Þ; of the cross-spectrum (equations (8) and (9)).
We then convert from a Cartesian wave number grid into polar coordinates: The dominant wave direction is identified where the maximum value in the real part of the cross-spectrum matches the maximum value of the imaginary part of the cross-spectrum (e.g., Bao & Alpers, 1998). The 1-D spectra are found by integrating over the angular coordinate of the maximum value of the imaginary cross-spectrum θ max , with a range of ±5°, using trapezoidal numerical integration: We integrate over a range of 10°angular coordinates, and the resolution of the wave number, k, is set to 0.0015. The choice of both the wave number resolution and the range in the angular coordinate should be as small as possible, so we only look at one set of waves. At the same time, the range should be large enough to smooth the noise as much as possible. After some trial and error, we found that the chosen values provided acceptable results.

Step 4: Compute the Wave Dispersion
The nature of the ocean-to-SAR transform depends on the imaging regime, which can be linear, quasi-linear, and nonlinear (Krogstad, 1992). In the linear imaging regime, the 2-D SAR image cross-spectrum is related to the ocean wave spectrum through (Engen et al., 2000) P k; θ; τ ð Þ≈T 2 k; θ ð Þ·S k; θ ð Þ·e iφ k;τ ð Þ ; where T(k, θ) is the MTF and S(k, θ) the ocean wave spectrum. The linear regime does not occur very often in the open ocean as it requires a low significant wave height (and hence, a low azimuth cutoff) and relatively long waves. Within the sea ice, however, the azimuth cutoff is much smaller because the short wavelengths of wind-driven seas, which are the dominant contribution to the azimuth cutoff effect, are mainly absent because the ice cover acts as a low-pass filter (Collins et al., 2015). Therefore, we argue that SAR images taken over the sea ice are in the linear imaging regime. With this assumption, we can compute the phase 10.1029/2019JC015311 spectrum, φ(k, τ), from the imaginary part and real part of the 1-D crossspectrum given in equations (10) and (11). In this case, the MTF from equation (12) cancels out, and we get Both the imaginary part and the real part are smoothed using the moving average filter. The relation between the phase spectrum and the angular wave frequency, ω(k), is The look separation time between the images, τ, depends on the subswath and the range coordinate (see Table 2). For a certain range coordinate, the look separation time is found by a linear interpolation between the near and far range coordinates.

Test Case
To demonstrate the abovementioned steps of our proposed method, we give an example within the sea ice of the wave event that occurred on 4 April 2017. This event is characterized by relatively long waves (peak wavelength of 190 m; see Table 1). Two subsequent bursts for the same subswath are presented in Figures 4a and 4b. The red line in Figures 4a and 4b is found from the azimuth time and indicates where the next burst starts and where the previous burst ends, respectively. In Figure 4c, the normalized intensity averaged over the range is plotted for two overlapping bursts. It can be clearly seen that the intensity of a burst fades at the beginning and at the end of each burst. The length between the two thick black lines satisfies the criterion given in equation (4) and is the number of pixels in the azimuth direction. A typical obtained sub-image is presented in Figure 4d, which is the result of Step 2 in Figure 3.
The real part and the imaginary part of the cross-spectrum from two sub-images (the green rectangles in Figures 4a and 4b) are shown Figures 4e and 4f, respectively. The resolution of the cross-spectrum in the range direction dk x ¼ 2π 1000·5 À Á is 0.0013 and 0.0015 [rad/m] in the azimuth direction dk y ¼ 2π 300·14 À Á . The resolution in the azimuth direction is slightly different from burst to burst, as the number of points used in azimuth differ slightly for each sub-image. In Figure 4g, the phase spectrum is given, computed with equation (13). The resulting 1-D real and imaginary image spectra are presented in Figure 4h, which is the result of Step 3. In addition, the computed angular frequencies from Step 4 are plotted and compared with the theoretical open-water dispersion relation, which uses data on the currents collected from Ice Mooring Stations IC3-IC5 at a depth of 8 m (Table 1). For this example, the computed angular frequency from S1 is almost identical to the theoretical open-water angular frequency in the range where the wave energy is located (between approximately k = 0.025 and 0.04).
The portion of the imaginary and the real parts of the cross-spectrum that is used to calculate the phase, equation (13), does include noise. Therefore, we can rewrite equation (13) We estimated the noise of the imaginary spectrum (ImNoise) and the noise of the real spectrum (ReNoise) by calculating the variance in a small area of the cross-spectrum, far away from the area of interest. For example, in the event shown in Figure 4, we estimated the variance in the area: k x between 0.6 and 0.65 and k y between 0.2 and 0.25. For all the cases considered in this paper, the estimated noise was less than 1% of the peak values of the cross-spectrum, and hence, noise did not greatly influence the results of the wave dispersion.
To illustrate the advantage of our proposed method in which the burst overlap area is used compared to the single-look area, we compute the imaginary spectrum for the case within the sea ice as presented in Figure 4 using the single-look area. The result is presented in Figures 5a and 5b. The quality of the spectrum derived Journal of Geophysical Research: Oceans from the burst overlap area is much better and less noisy. For this example, the wave direction can still be seen from the spectrum of the single-look area. However, in the open ocean, the imaginary spectrum becomes even more noisy, making it sometimes impossible to derive the wave direction accurately from the single-look area. This result can be greatly improved by using the burst overlap area instead, where we have a larger time separation that results in a better-quality imaginary spectrum. An example in the open water is presented in Figures 5c and 5d, which shows the imaginary spectrum at the location of Wavescan Mooring WS3. The measured 2-D wave spectrum at this location is given in Figure 5e. The imaginary spectrum of the single-look area (Figure 5c) is in this case too noisy to extract the wave direction and wavelength. However, these wave parameters can be retrieved from the imaginary spectrum from the burst overlap area (Figure 5d) and have values of 67°and 217 m for the mean wave direction and the peak wavelength, respectively. These values are close to the ones obtained from the measured 2-D wave spectrum (Figure 5e), which has a mean wave direction of 64°and a peak wavelength of 193 m.

Applicability of the Method-Linear and Nonlinear Imaging Regime
Among the wave events considered in this paper, the angular frequency as a function of the wave number, that is, the wave dispersion relation, is underestimated in open water, except for Wave Event 1. We believe that this underestimation may be attributed to nonlinear imaging effects caused by the random motions of short waves. These random motions blur the image and lead to the azimuth cutoff effect (e.g., Kerbaol et al., 1998). Linearity of ocean-to-SAR transform is assumed when computing the phase from equation (13) in our method, and therefore, we do not expect our method to work in the nonlinear imaging regime. The black arrow indicates the north, and the circles show the wavelength. Note that the imaginary spectrum of the burst overlap area within the sea ice (b) is identical to the one presented in Figure 4. In (e), the measured 2-D wave spectrum measured at wavescan mooring WS3 is presented.

Journal of Geophysical Research: Oceans
To demonstrate this, we apply our proposed method to waves propagating from open water into sea ice, which allows us to derive and compare the wave dispersion relations for linear and nonlinear cases. We present the derived wave dispersion relations both in the open water and within the sea ice for two wave events in Figure 6. These two wave events are characterized by long waves and occurred on 4 April 2017 and on 16 February 2018 (Events 1 and 2 in Table 1). The derived values of the angular frequency are compared with the theoretical open-water dispersion relation, which is calculated using in situ data of the currents. Note that the comparison by means of the theoretical open-water dispersion relation is also valid for images taken over the sea ice since the dispersion of long waves within the sea ice is nearly identical to the theoretical open-water dispersion relation (Cheng et al., 2017;Collins et al., 2018). In Figure 6, the wave dispersion obtained from S1 is in good agreement with the theoretical open-water dispersion relation for images taken over the sea ice for both events. However, the obtained angular frequencies in the open water show an underestimation for the 4 April 2017 event, while again, a good agreement is found for the 16 February 2018 case.
Factors that could lead to differences between the results from S1 and the theoretical open-water dispersion relation could be a too noisy imaginary spectrum. This is, for instance, the case when instead of the burst overlap area, the single-look area is used. Another factor could be that our assumption of linearity of the ocean-to-SAR transform is violated, which is the case for the 4 April 2017 event in the open water. To illustrate this, we estimated the linearity for both cases by calculating the nonlinearity coefficient, equation (16), proposed by Alpers et al. (1981), which was used in the study by Shen et al. (2018) to show linear conditions. This is the case when the value is equal or smaller than 0.3.
with R, V, and ϑ the radar parameters that are the target range, platform velocity, and the radar incidence angle, respectively, a is the wave amplitude, k p the peak wave number, and φ p the peak wave direction relative to the azimuth direction. The computed nonlinearity coefficient is 0.  minimizing the error of a fitted Gaussian function to the normalized azimuth profile of the cross-covariance function. This method minimizes the function (Johnsen & Collard, 2009): where λ c is the estimated azimuth cutoff wavelength. This method has been indicated to give reasonable results over the open ocean (Kerbaol et al., 1998;Stopa et al., 2015). The computed azimuth cutoff wavelength is approximately 306 m for the 4 April 2017 event and approximately 260 m for the 16 February 2018 event. On 4 April 2017, there was a higher wind speed and a larger significant wave height, resulting in a higher value. For the 16 February 2018 event, the peak wavelength (350 m) is well above the azimuth cutoff wavelength. In contrast, for the 4 April 2017 event, the peak wavelength (160 m) is smaller than the azimuth cutoff wavelength, and the peak wavelength is distorted by the cutoff. As stated in section 4.5, we argue that the images taken over the sea ice are in the linear regime. Most of the short waves are not present in images taken over the sea ice because the sea ice scatters and dissipates the highfrequency waves. For that reason, waves appear clearer in images taken over the sea ice than in images taken over the open ocean  and that the azimuth cutoff wavelength is much smaller in the sea ice than in the open water. As a result, only the images taken over the open water for the 4 April 2017 event, presented in Figure 6a, are in the nonlinear imaging regime, and the peak wavelength is distorted by the azimuth cutoff. Therefore, we conclude that our method does not work when the assumption of linearity is violated. The good results of the open water case for the 16 February 2018 event ( Figure 6c) and the cases within the sea ice, Figures 6b and 6d, give us confidence in the validity of our method when we are in the linear regime.

Analysis of the Derived Wave Dispersion Relations for Long Waves Within Sea Ice
We validate our method by comparing the obtained angular frequencies with the values from the theoretical open-water dispersion relation because the wave dispersion of long waves within sea ice is nearly identical to the open-water wave dispersion relation. We use three different wave events with long waves, that is, Wave Events 2-4 in Table 1. The results are shown in Figure 7. The results are in good agreement with the theoretical open-water dispersion relation in the area where the wave energy is located. We evaluated the performance by calculating the mean absolute error and the mean absolute percentage error (MAPE). These errors are defined as follows: where ωs1 is the angular frequency computed from S1 and ωow is the angular frequency calculated from the theoretical open-water dispersion relation. The MAPE for the three wave events, based on the average of the observations, is approximately 1%. However, quite a large spread in the observations can be observed. This spread is roughly ±7% when looking at the right panels of Figure 7. Nonetheless, the good agreement with the theoretical open-water dispersion relation gives us confidence in the validity of our method.
An advantage of this study is the availability of in situ data on the currents. We investigated the influence of the currents by comparing the three wave events presented in Figure 7 with the theoretical open-water dispersion relation, without using the currents. The MAPE, when including currents, is approximately 1-2% lower compared to the case without currents. For these specific cases, the error does not change substantially because the magnitude of the currents is quite weak, with values between 0.1 and 0.2 m/s. Having stronger currents will modify the open-water wave dispersion considerably. To demonstrate this, we computed the MAPE for a theoretical case using the maximum observed current speed of 1.06 m/s. The MAPE is then up to 11% higher than the case with the actual observed current speeds, which clearly shows it is vital to have accompanying data on the currents.

Analysis of Derived Wave Dispersion Relations for Short Waves Within Sea Ice
The obtained dispersion relations of the three wave events with the shortest waves observed, that is, Wave Events 5-7 in Table 1, are plotted in Figure 8. The different wave events shown here include fewer observations than the ones with long waves because there are significantly fewer sub-images in which waves are visible within the sea ice. Furthermore, the spread in observations is larger compared to the events with long waves (Figure 7) because all the sub-images of the events with short waves are close to the ice edge. The ice here is very inhomogeneous; for example, broken-up ice and resulting ice floes in addition to patches of open water are found here. Nevertheless, the average of the observations does show a small deviation from the open-water dispersion relation. Recall that the derived dispersion relationship for long waves within the thin ice did not deviate from the open-water dispersion relation, which is consistent with findings in literature.
In addition to comparing the derived dispersion relation with the theoretical open-water dispersion relation, we attempted to back-calculate the ice thickness by looking at the simplest dispersion relation model, including sea ice, that is, the mass-loading model. To do so, we plotted this mass-loading model as a function of the ice thickness in addition to the derived dispersion relation from S1 in Figure 8. The ice thickness at Mooring Stations IC3-IC5 is approximately 0.2-0.5 m. For such a small ice thicknesses, we can see that the deviation from the theoretical open-water dispersion relation is minor according to the mass-loading model. It is evident from Figure 8 that the desired small deviation from the theoretical open-water dispersion cannot be sensed by the satellite, and we cannot use these observations to quantitatively extract the ice thickness.
To obtain a deviation sufficiently large from the theoretical open-water dispersion relation that can be sensed by S1, we would need to study even shorter waves. Collins et al. (2018) found that the deviation from the open-water dispersion relation for frequencies larger than 0.3 Hz (corresponding to a wavelength of approximately 14 m) was significant. Therefore, we would need to measure in this range of wavelengths to obtain a substantial deviation that can be sensed by S1. However, these measurements are not currently possible because the spatial resolution of S1 is too coarse. Having higher-resolution data will not only allow the study of shorter waves but will also lead to a better accuracy in determining the wave dispersion. The reason for this is that we can average more over the same spatial area, which will lead to a less noisy imaginary spectrum and therefore a better estimate of the phase and the wave dispersion.

Conclusions
In this study, we present a new method for retrieving spatiotemporal observations of ocean waves from the S1 constellation. We propose a new and innovative implementation of the method proposed by Johnsen and Collard (2009) on the spatial burst overlap area, which is present due to the S1 TOPS acquisition mode. The advantage of this method is that the time separation between two individual looks (images) is significantly larger, and, therefore, a larger phase change is present. This method leads to a less noisy and better-quality imaginary spectrum, which can potentially be used for a better determination of the wave direction. Furthermore, the proposed method does not require the MTF because we assume a linear ocean-to-SAR transform within the sea ice. The assumption is based on the fact that the ice cover acts as a low-pass filter, attenuating high-frequency waves. Short wavelengths of wind-driven seas are the dominant contributors to the azimuth cutoff; hence, this cutoff is much smaller within the sea ice.
We used seven wave events in the Barents Sea to test and validate our method. For this study site, we have accompanying in situ measurements of the sea currents and the ice draft. The use of the current measurements is a major advantage, as many previous studies that focused on the change in the wave number in the MIZ lacked simultaneous data on the currents. For our studied cases, the use of the currents generally improved the results by reducing the MAPE by 1-2%. This relatively small contribution of the currents is because the current speed during our studied events are low. However, it is evident that strong currents significantly modify the wave dispersion relation and can lead to considerable errors if not included. The computed wave dispersion within the sea ice for long waves (peak wavelengths larger than approximately 100 m) is in good agreement with the theoretical open-water dispersion relation. This agreement gives us confidence in the validity of our assumption of linearity and confirms that the dispersion relation is not altered for these waves, which agrees with previous findings. In the open water, the estimated angular frequency, as function of the wave number, is underestimated for almost all the studied wave events. We show that the reason for this is that the assumption of linearity is violated in these cases.
Three wave events were studied with peak wavelengths in the range of 50-70 m, which were the smallest wavelengths we found where waves were visible within the sea ice. We investigated whether it is possible to back-calculate the ice thickness using the mass-loading model. However, the combination of the small ice thickness observed at the ice mooring stations and the range of wavelengths studied causes the deviation from the theoretical open-water dispersion relation to be minor. This deviation is too small to be able to be sensed by S1. Moreover, the observations from S1 show quite a large spread due to the inhomogeneous ice cover close to the ice edge. Broken-up ice and resulting ice floes, as well as open-water patches, are present.
For the deviation to be sensible by the satellite, we would have to study waves with a wavelength of approximately 10-20 m, which requires higher-resolution data. Besides being able to study shorter waves, having higher-resolution data allows more averaging over the same spatial area, which eventually will lead to an accuracy sufficient to detect the effects of thin ice on the wave dispersion of short waves.