SAR Altimetry Data as a New Source for Swell Monitoring

This article shows the first spectral analysis of fully‐focused Synthetic Aperture Radar (FFSAR) altimetry data with the objective of studying backscatter modulations caused by swells. Swell waves distort the backscatter in altimetry radargrams by means of velocity and range bunching. These swell signatures are visible in the tail of the waveform. By locally normalizing the backscatter and projecting the waveforms on an along‐/cross‐track grid, satellite altimetry can be exploited to retrieve swell information. The analysis of FFSAR spectra is supported by buoy‐derived swell‐wave spectra of the National Oceanic and Atmospheric Administration network. Using cases with varying wave characteristics, we discuss the altimetry‐derived spectra and relate them to what is known from side‐looking SAR imaging systems. Besides having a vast amount of additional data for swell‐wave analysis, altimeter data can also help us to better understand the side‐looking SAR spectra.

the high-rate sampling of FFSAR data may be suitable for detecting swell-related information and, therefore, for explaining such distortions in the radar backscatter signal.
Since the seventies the SAR and Real Aperture Radar (RAR) imaging mechanisms have been investigated in detail as they have been useful tools for monitoring ocean surface waves (Alpers & Rufenach, 1979;Ardhuin et al., 2017;Larson et al., 1976;Raney, 1971). The main RAR modulations that are responsible for the imaging of ocean waves are the tilt modulation, which is related to local changes to the incident angle by slopes of long waves, range bunching, which is the variation of surface area captured in resolution cells as a consequence of the same slopes, and hydrodynamic modulation, which is triggered by the interaction between the short and long waves (Alpers et al., 1981). The SAR processing introduces one more mechanism, called velocity bunching (Alpers et al., 1981), that seems to dominate in most cases . The velocity bunching is the clustering of the signal due to the vertical motion of the waves. In more detail, the phase-history difference of echoes between nearby along-track points is approximately linear. The motion of a scatterer in the range direction, which corresponds to vertical motion for an altimeter, can replicate this linear behavior causing the projection of the scatterer to be misplaced in the along-track direction. As scatterers ahead and behind the crest of a wave move up and down, respectively, they are both shifted in a SAR image in opposite directions. Depending on the direction of the wave with respect to the radar motion, scatterers are clustered either around the crest or the trough. With a spectral model, containing both RAR and the shift due to velocity bunching, a SAR intensity spectrum can be inverted to obtain part of the wave spectrum (K. Hasselmann & Hasselmann, 1991). The inversion is limited to long waves as the effective resolution of SAR observations over ocean is far less than the theoretical limit (Alpers & Hasselman, 1978;K. Hasselmann et al., 1985), with typical values of 100 m. In contrast to side-looking SAR systems, altimeters are nadir-looking radars. This difference has an impact on the backscatter mechanisms which prevail (Valenzuela, 1978) and also geometrically changes the RAR response. The so-called Bragg, or resonant, scattering dominates for a side-looking system at moderate sea states, while the specular scattering for a nadir-looking system. This paper shows that SAR altimetry data have the potential to infer swell-wave spectra. We introduce a methodology to process altimetry data to estimate a FFSAR spectrum, which is suitable for an inversion in a comparable way as for side-looking SAR data. Moreover, this paper gives a detailed description of the modulation process, shows how this differs from a side-looking system and discusses the adaptions required for a nadir-looking system. Lastly, buoy data are used to support the analysis.

Data and Methods
This section is divided into two parts. The first part introduces the processing steps to obtain a FFSAR spectrum. The second part shortly discusses the buoy data.

SAR Altimetry Data Processing
A flowchart of the main processing steps is given at the top of Figure 1, while an example accompanies it to better understand the behavior of the backscatter altimetry signal. The latter concerns data from a descending orbit of CryoSat-2 over the Channel Islands of California on 2 January 2020. The presence of swell waves, with a period of 20 s and direction 281° with respect to North, is confirmed by a nearby buoy of the National Oceanographic and Atmospheric Administration-National Data Buoy Center (NOAA-NDBC) network (https://www.ndbc.noaa. gov/). The FFSAR L1b multilooked waveforms used in this study were obtained from the ESA Grid Processing-on-Demand (GPOD) FFSAR service. This service is based on the Omega-Kappa SAR focusing algorithm described in Guccione et al. (2018). For the demonstration, we use radargrams of 500 waveforms of 200 Hz multilook posting rate, which corresponds to an along-track sampling of about 34 m, with a total sample length of about 17 km. Note, that if single-look waveforms are available, we may opt a weighted multilooking, because it acts as a filter and this, in turn, may have an impact on the FFSAR spectra. The first panel, i.e., Figure 1a, depicts a radargram. A closer look at the tail of it reveals undulations that are caused by swells.
For an accurate application of the method, the waveforms in the radargram should be aligned. A change of the leading edge of one bin leads to a misprojection of tens of meters across track. Small variations of the leading edge location can be dealt with in the ground-projection step, but larger variations can lead to normalization problems. A suitable realignment strategy is a two-step process: (a) the waveforms are averaged over a distance larger than the along-track-projected swell wavelength and retracked, and (b) consecutively a phasor is applied to the complex values in FFSAR processing after the range-residual-phase correction to realign the retracked bin to a pre-defined reference. As for the cases considered in this study the waveforms are aligned to less than a bin and no wave-spectrum inversion is performed, realignment is not applied.
The first step of our method is the normalization of the radargram intensity ( Figure 1b). This is required to compute the normalized FFSAR spectrum at a later stage so that it can be used in a wave spectrum retrieval algorithm comparable to the one used for side-looking SAR systems, such as Sentinel-1. In the case of Sentinel-1, the intensity of a scene of 20 × 20 km (wave mode) is normalized by first subtracting the mean intensity and then a division by the mean intensity. However, due to the drop of the power in the trailing edge of the SAR waveforms, such a normalization is not suitable. We overcome this by first computing the expected intensity ̂ for each bin in the radargram. The expected intensity is used to compute the normalized intensity I N , that is, where I indicates the radargram intensity. Note that, I should resemble the scaling of a SAR image. The expected intensity ̂ is computed from the data by first applying a Gaussian filter in along-track direction with a width (defined as two times the standard deviation) of 25 waveforms. This corresponds to about 850 m along track. The filter length is cut at 51 waveforms (a sensitivity analysis is done supportive material in Supporting Information S1). Then we fit a polynomial of degree 4 to the tail of each waveform. The tail includes all bins between 4 and 7 km across track. Finally, ̂ is obtained by evaluating the polynomial at the waveform bins, and subsequently used to compute I N .
In the second step, the normalized waveforms are projected on the ground. The distance between the cross-track points is determined as where l cross indicates the horizontal distance between the range bins projected on the ground, dR is the range sampling interval of 0.234 m, n is the range bin, n ref the leading edge reference bin and h is the mean satellite altitude for the area of interest. In along-track direction, the distance between successive points is determined by the along-track sampling. To relate each range bin to a unique cross-track location the radargram is projected as if all backscatter is received from the right side. Consequently, the projected normalized radargram contains a correctly projected backscatter signal of the right side and a left-right flipped backscatter signal from the left side. Backscatter distortions, caused by a monochromatic swell wave, are therefore expected to exhibit a crossing pattern. In the last step, we compute the FFSAR spectrum as the square of the Discrete Fourier Transform (DFT) amplitude. Prior to the DFT operation, the projected data are interpolated to an equidistant grid with a resolution of 10 m. Note that the DFT is only applied to the trailing edge (corresponding to 4-7 km across track) of the waveforms (Figure 1c). To suppress the speckle noise in the FFSAR spectra, a Gaussian filter with a width of four pixels (two standard deviations) is applied to it. The resulting spectrum for the example track is given in Figure 1d. Although the common way to suppress noise is to use cross-spectra computed by subswath processing (Engen & Johnsen, 1995), this is currently not feasible with the GPOD FFSAR processor.
We have to note that with the transmission of pulses by CryoSat-2 and Sentinel-3A/B, ghosts appear at approximately 90 m along track (Egido & Smith, 2017). The spectral response under swell conditions is normally more than one order larger than the spectral response of the ghosts. In case of low sea states, the spectral response of the ghosts cannot be ignored and should be appropriately handled.

Buoy Data
Buoy data from the NOAA-NDBC network are used to verify and discuss the FFSAR spectra. From these publicly available data we collected the swell-wave spectra along with the peak period (i.e., the period with the maximum wave energy) and the direction from which the waves with the corresponding peak period are coming. The data have an accuracy of ±1 s and ±10°, for the peak period and direction, respectively.
We demonstrate the method for an area in the Northeast Pacific Ocean. This area is chosen because there is often swell present and a buoy is located in offshore waters to limit the effects of coastlines and bathymetry (depth ranges from 100 to 1,100 m). In this area CryoSat-2 operates in SAR mode. CryoSat-2 data were considered within a radius of 50 km from the buoy. Four cases, identified as having swell systems are discussed here. The upper panel of Figure 3 illustrates a map with four CryoSat-2 tracks, two descending and two ascending, as well as the buoy location. To limit effects of spatial and temporal variations between the buoy and FFSAR spectra, the time needed for the swell system to cross both measurement locations should be taken into account. This time offset, Δt, was calculated as where Δx indicates the distance calculated from the middle point of each overpass to the buoy, c g is the group velocity calculated using linear wave propagation theory and considering intermediate conditions, ϕ w is the wave direction with respect to the North, and ϕ m the relative direction of the satellite track and the buoy with respect to the North. Then, based on the relative position of the satellite and the buoy with respect to the wave direction, i.e., whether the waves cross first the satellite observation location or the buoy, the time offset was accordingly added or subtracted from the acquisition time of the satellite data in order to select the buoy epoch.

Hasselmann's SAR Spectrum
To support the interpretation of the obtained FFSAR spectra, we introduce the analytical model relating the SAR spectrum to a wave spectrum presented by K. Hasselmann and Hasselmann (1991).
which is a function of the cross-track wave number k x and the along-track wave number k y . This spectral mapping depends on the characteristic function (Krogstad, 1992) with ρ yy being the correlation of azimuthal shifts as a function of wave spectrum S, and ρ II the correlation function of linear RAR modulations containing range bunching, tilt and hydrodynamic modulation also as a function of wave spectrum S. Higher-order terms, H, are ignored as it simplifies the equation, and as their contribution is an order smaller than ρ II (Engen & Johnsen, 1995;Krogstad et al., 1994). The spectral transform and its correlations are further elaborated in Engen and Johnsen (1995) Schulz-Stellenfleth and Lehner (2002), and others. Note that, as the characteristic function depends on wavenumbers k x and k y , the integral Equation 4 is not a Fourier Transform. In fact, the integral should be evaluated for each wavenumber separately. In this paper this simple form of equation (H. Li et al., 2019) is dissected into separate terms to support the interpretation of FFSAR spectra.

SAR Modulation
First, the − 2 ( (0,0)) term (referred to as "term A") can be taken out of the integral. The term ρ yy (0, 0) represents a scaled velocity variance (Kerbaol et al., 1998;Schulz-Stellenfleth & Lehner, 2002). At higher sea states, the velocity variance increases. Term A causes a decrease of SAR spectral power with increasing k y , which represents the reduction of along-track resolution as a consequence of random linear motions of scatterers. Further resolution loss can occur due to higher-order motions of scatterers, which decrease the decorrelation time (Alpers & Hasselman, 1978;Kerbaol et al., 1998).
Second, the term 2 (referred to as "term B") represents, together with the first term, the velocity bunching. The second term, which is non-linear, increases with an increasing angle of wave propagation (here the angle is defined with respect to the cross-track direction). Under some circumstances, for example, weak wind and strong swell, secondary responses appear in the SAR spectrum (Schulz-Stellenfleth & Lehner, 2002). The exponential term becomes equal to 1 in case k y = 0, which implies that velocity bunching does not occur with waves traveling in the cross-track direction.

RAR Modulation
The last term, ρ II , (referred to as "term C") represents the first-order RAR response and can be considered as the response of the ocean surface if it was not moving during the time of overpass. It is affected by range bunching, tilt, and hydrodynamic modulations (Alpers et al., 1981;K. Hasselmann et al., 1985;Schulz-Stellenfleth et al., 2005). It is different for side-looking and nadir-looking systems as described in the introduction. In contrast to the velocity bunching, RAR modulation also occurs with cross-track traveling waves (B. Chapron et al., 2001). The discussion on the RAR modulation is supported by simulation results, provided in supportive material S3 in Supporting Information S1.
In Figure 2 a geometric representation of a cross-track propagating swell wave is shown. Hydrodynamic interaction of long waves with shorter ones causes an increase of roughness just after the crest, which reduces the specular reflection. The maximum specular reflection caused by hydrodynamic interaction is therefore expected just behind the trough. Tilt causes an increase of returned power where incident signal is close to normal to the surface, which is also close to the point where range bunching is maximal as at that point the largest surface area is captured between two range isolines. Tilt and range bunching are maximal on the slope facing the satellite. For small significant wave heights there is one maximum near the slope maximum (see supportive material figures S3.2A/B in Supporting Information S1, but if the wave slopes are larger than the incident angle, two maxima appear that move toward the crest and the trough (see supportive material figures S3.3A/B in Supporting Information S1). On one side of the satellite, the range bunching and tilt occur on same slope as hydrodynamic modulation, but on the other side not (Figure 2b).
For a near-nadir looking radar system the hydrodynamic modulation is considered small compared to range bunching and tilt (H. Li et al., 2021) and is therefore not considered. For swell waves considered in this paper, tilt modulation is much smaller than range bunching (see supportive material S3 in Supporting Information S1) and will therefore also be ignored in the further analysis. Note that the range bunching is very non-linear as the slopes of the swell waves are typically larger than the incident angle. Even for a monochromatic wave the spectral response will be smeared as we might have two maxima on one satellite-facing slope of the wave with a distance between them that depends on wavelength, amplitude, and direction (see supportive material figures S3.2A-S3.3B in Supporting Information S1). When waves propagate in the along-track direction the range bunching is nearly zero. A small signal, though, remains as crests and troughs are captured in different range bins, which have a different projected cross-track resolution. This complexity is further enhanced by overlay (Figure 2a), which makes it difficult to derive a closed-form expression for the spectral response of range bunching. At slopes 10.1029/2021GL096224 6 of 11 smaller than the incident angle, the range bunching might be approximated by the equations used for the Chinese-French Oceanography Satellite (CFOSat) (H. Li et al., 2021;Jackson, 1987).
Due to radial velocities the RAR modulated signal are also moved in the along-track direction and get bunched (Figure 2b). If only range bunching and velocity bunching are considered and the waves are symmetric with Figure 2. Geometric representation of the Real Aperture Radar (RAR) imaging modulations in nadir-looking systems as a function of the altimeter cross-track footprint from 3,000 to 4,000 m for a 450 m swell wave (a). Light and dark blue dashed lines represent the troughs and crests, respectively. The maxima of the hydrodynamic modulation which results from reduced roughness near the troughs are indicated with red dots. The middle panels of (a) show the approximate sinusoidal behavior of the hydrodynamic modulation. The maxima of the range bunching (green dots) and tilt modulation (black dots) coincide with minimum local incident angle (bottom panels of a). Range bunching is maximal where the surface is aligned with the range isolines (light gray lines, upper panels of a). Top view of the corresponding RAR and Synthetic Aperture Radar (SAR) modulations (b). The SAR modulation, i.e., velocity bunching, is maximal either near troughs or crests. The change of focal location due to surface motion is indicated with the black arrows. Therefore for this specific case, the velocity bunching is maximal near the crests.
respect to the crest, the expected modulation is out-of-phase left and right of the track whereas the expected spectral power is the same. However, it is a statistical process, which causes a deviation from the expected value (see supportive material figures S3.4A/B and S3.7A/B in Supporting Information S1) and a difference between the spectral amplitudes left and right. Simulations indicate that velocity bunching and range bunching have the same order of magnitude, with range bunching often dominating (see supportive material S3.5A-S3.6B and Figure 3. A map of the study area is given on the top. The red and blue lines represent the descending and ascending CryoSat-2 tracks, respectively, while the black triangle indicates the buoy with World Meteorological Organization number: 46219). Below, each column of panels show the projected radargrams, the corresponding fully-focused Synthetic Aperture Radar (FFSAR) spectra and the swell-wave spectra obtained from the buoy for two descending (cases 1 and 2) and two ascending (cases 3 and 4) overpasses. The red crosses, plotted on top of the SAR spectra, represent the wavenumber vector with the maximum wave energy based on the buoy measurements. The acquisition date along with the in situ observations (SWH, significant wave height; DWP, dominant wave period; DWD, dominant wave direction) and the mean roll angle of the satellite are given at the top. The SWH is calculated as the average of the highest one-third of all of the wave heights.

Results and Discussion
The panels of the second row in Figure 3 show the projected radargrams used for the computation of the FFSAR spectra, which are presented in the panels of the third row. Note that the axes of the latter have been reversed to easily compare the FFSAR spectra with the swell-wave spectra given in the bottom panel. Four peaks are observed, one per quadrant. This differs from a SAR spectrum computed by a side-looking SAR system, where only a signal of (k x , k y ) is mirrored at (−k x , −k y ). The ambiguity of 180° in a side-looking system is normally accounted for by using the imaginary part of a cross spectrum (Engen & Johnsen, 1995), which can also be done by subswath processing of altimetry data. In FFSAR spectra two additional peaks, albeit with less power, are present in quadrants that do not exhibit substantial wave energy. As the altimeter observes both sides of the ground track at the same time, the expected value will get approximately mirrored around the ky-axis if velocity bunching and range bunching are dominant. In that case, the (double-sided) FFSAR spectrum, P 2s (k x , k y ), is approximated as where P l and P r the approximated spectra from the left-and right-hand sides of the ground track, and P(k x , k y ) as in Equation 4. γ is a weighting factor close to 0.5 depending on platform roll, the antenna pattern, but can be also influenced by cross-track waves and the skewness in slope probability distribution function (Munk, 2009).
In three of the four shown cases the signals in two quadrants are substantially weaker than expected. This can partly be explained by a statistical deviation from the expected value (see differences between the left-and rightside spectra in supportive material figures S3.7A/B in Supporting Information S1). There are also other possible deterministic causes for the discrepancy. First, there might be geophysical effects, such as the interaction of waves with currents and bathymetry that cause a cross-track difference in the backscatter modulation between both sides of the track. Second, there are instrumental effects such as the roll of the satellite that cause the gain to vary between both sides of the satellite, which is indicated with the term γ in Equation 6. However, from the four cases shown here and in supportive material S2 in Supporting Information S1, it appears that the return power in the quadrants the waves are actually moving is higher than in the other two quadrants.
The latter indicates an asymmetric upwave-downwave response. The model for velocity bunching and range bunching alone (see supportive material S3 in Supporting Information S1) is not able to reproduce a deterministic asymmetric response. Therefore either the wave properties are not symmetric and/or mirrored around the crest (i.e., slope and radial velocity) or a substantial hydrodynamic modulation would be required. At the moment of writing this asymmetry is not fully understood. Note that, H. Li et al. (2021) also described an up-to-down wave asymmetry in the fluctuation spectrum computed from CFOSAT data with a larger spectral response in the upwave direction.
By comparing the FFSAR spectra with the buoy-derived swell-wave spectra it is observed that the maximum energy closely coincides for all four cases, but the wavenumber difference between ocean-wave spectra and SAR spectra peaks is generally non-zero. Discrepancies can occur due to the non-linear nature of the SAR and RAR responses. A reduction of signal might occur in the cross-track direction, as velocity bunching (term B) goes to zero and the spectral response caused by range bunching (term C) is typically a bit smaller in the cross-track direction as the response smears over more wavenumbers (see supportive material S3 in Supporting Information S1). An approach of the cut-off wavelength (related to term A) in the azimuth direction, decreases the response at increasing k y . Waves with a relatively high along-track wavenumber will therefore be dampened. The (along-track) resolution depends on the sea state, or more precisely the velocity variance (term A) of the ocean (see section about the SAR spectrum), but also on the coherence time of scatterers (Alpers & Rufenach, 1979;Buchhaupt et al., 2021). The typical coherence time of a C-band SAR imaging system (e.g., Sentinel-1, Envisat and others) in open water is 50 ms (Carande, 1994). At 50 ms the velocity variance is the dominant driver for resolution loss. The coherence time for a Ku-band altimeter, with its shorter wavelength, is smaller and therefore the resolution loss is likely more severe. Lastly, small differences between the location of peak energy in the buoy-derived swell-wave and FFSAR spectra may occur because of the non-linear responses of SAR, which are captured by higher order terms in the G-function (Equation 5) for which RAR and SAR are not independent.

Conclusions and Recommendations
This paper presents a method to compute SAR spectra from FFSAR altimetry data and discusses the spectral response with the help of analytical and numerical modeling. FFSAR radargrams of CryoSat-2 data show modulations of power in the trailing edge of the waveforms in the presence of swells. A normalization of the intensity and a re-projection of the range bins as a cross-track location allow for a spectral analysis in a similar way as for side-looking SAR systems. Features of the FFSAR spectra are compared to a model developed for side-looking SAR data.
Side-looking SAR spectra only have an ambiguity of 180°. FFSAR altimetry spectra, however, exhibit power in four quadrants due to signals received from both sides of the ground track. In addition, the RAR response differs from that in a side-looking SAR system in terms of geometry and scattering mechanism. The largest term of RAR, range bunching, appears to be at the same order of magnitude and often even larger than the velocity bunching spectral response, depending on the wave height and direction. The inversion of an FFSAR altimetry spectrum into an ocean-wave spectrum requires a better understanding of the RAR response.
SAR altimeters show potential as an additional means to retrieve long-wave spectra. Having currently four SAR altimeters in orbit, i.e., CryoSat-2, Sentinel-3A/B and Sentinel-6, the data set available for studying swell will be greatly expanded. As the swell waves are observed from two different instruments operating at different frequencies the SAR and RAR responses can be studied in more detail by looking at cross-overs between side-looking SAR and altimeters. Besides that, the sampling is greatly enhanced, which is especially beneficial for the study of waves radiating from tropical cyclones. Finally, the scaled integral of FFSAR derived long-wave spectra provides an estimate of significant wave height for swell only. In combination with the significant wave height obtained from low-resolution mode altimetry, we can discriminate between swell-and wind-wave heights.

Data Availability Statement
The buoy data are publicly accessible through the National Oceanographic and Atmospheric Administration-National Data Buoy Center from https://www.ndbc.noaa.gov/. The fully-focused altimetry data that support the findings of this study are openly available in the 4TU.ResearchData repository of TU Delft (https://data.4tu.nl/) at https://doi.org/10.4121/17198435.