# Cross-Spectral Analysis of SAR Altimetry Waveform Tails

Marcel Kleinherenbrink, Frithjof Ehlers, Sergi Hernandez, Frederic Nouguier, Ourania Altiparmaki, Florian Schlembach, Bertrand Chapron

► To cite this version:

Marcel Kleinherenbrink, Frithjof Ehlers, Sergi Hernandez, Frederic Nouguier, Ourania Altiparmaki, et al.. Cross-Spectral Analysis of SAR Altimetry Waveform Tails. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62, 4206615 (15p.). 10.1109/TGRS.2024.3402390. hal-05063282

HAL Id: hal-05063282

https://hal.science/hal-05063282v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cross-Spectral Analysis of SAR Altimetry Waveform Tails

Kleinherenbrink Marcel <fnref n="1" />,*, Ehlers Frithjof <fnref n="1" />, Hernandez Sergi <fnref n="2" />, Nouguier Frederic <fnref n="3" />, Altiparmaki Ourania <fnref n="4" />, Schlembach Florian <fnref n="5" />,<fnref n="6" />, Chapron Bertrand <fnref n="3" />

* Corresponding author : Marcel Kleinherenbrink, email address : m.kleinherenbrink@Tudelft.nl

Abstract :

Until recently, intensity modulations in synthetic aperture radar (SAR) altimetry waveform tails have been considered a nuisance for geophysical-parameter retrieval. These modulations are actually predictable and might be exploited using a spectral analysis of the waveform tails. After Altiparmaki et al. (2022), a more elaborated analysis is performed to improve the interpretation of these SAR altimeter spectra. A fast numerical model is developed to explain the modulation mechanisms in focused SAR altimetry waveform tails. Using numerical solutions, standard analytical closed-form solutions, are demonstrated to be invalid to retrieve ocean-wave-spectra retrievals from nadir altimeters. Although not valid, a closed-form derivation provides intuitive insights about the information contained in an SAR altimetry cross-spectrum. Under moderate environmental conditions (significant wave heights (SWHs) of similar to 2 m), a closed- form solution might still be useful to infer swell-wave spectra from swath-altimetry SAR spectra at incident angles of similar to 4 degrees. Comparable to side-looking SAR ocean processing, the cross-spectral analysis for nadir signals reduces noise and might remove the 180 degrees ambiguity of the wave direction. Since the synthetic aperture length of nadir altimeters is larger than sidelooking imaging SARs (e.g., Sentinel-1, RadarSat, Gaofen-3), sublook processing can be performed to compute multiple cross- spectra for the same scene. With a slightly changing observation geometry, the cross-spectra reveal slightly different parts of the ocean-wave spectrum. The resulting stack of cross-spectra can thus be used to improve the retrieval of ocean-wave parameters. Retrieved ocean-wave parameters shall then enhance the sampling of the global wave field, but also serve to advance more consistent sea-state-bias corrections.

**Keywords** : Cross-spectral analysis, ocean-wave spectra, sentinel-6, synthetic aperture radar (SAR) altimetry, SAR spectra

<footnote n="1">¹ Delft Univ Technol, Dept Geosci & Remote Sensing, NL-2628 CD Delft, Netherlands.</footnote>

<footnote n="2">² IsardSAT, Barcelona 08005, Spain.</footnote>

<footnote n="3">³ Univ Brest, Inst Francais Rech Exploitat Mer, Lab Oceanog Phys & Spatiale, CNRS, F-29280 Plouzane, France.</footnote>

<footnote n="4">⁴ Delft Univ Technol, Dept Astrodynam & SpaceMiss, NL-2628 CD Delft, Netherlands.</footnote>

<footnote n="5">⁵ Tech Univ Munich, Deutsch Geodat Forsch Inst, D-80333 Munich, Germany.</footnote>

<footnote n="6">⁶ Rohde & Schwarz, Munich, Germany.</footnote>
# I. INTRODUCTION

Satellite radar altimetry has already provided observations of sea-surface height (SSH), significant wave height (SWH), and wind speed (WS) for a continuous period of more than three decades. The launch of CryoSat-2 with its delay/Doppler altimeter first largely improves the resolution with respect to standard low-resolution mode (LRM), down to ~300 m in the along-track direction [2]. Delay/Doppler, or unfocused synthetic aperture radar (SAR), processing secondly increased the number of independent looks, to help reduce the noise of

M. Kleinherenbrink and F. Ehlers are with the Department of Geoscience and Remote Sensing, Delft University of Technology, Delft, The Netherlands.

retrieved geophysical parameters. The resolution was further improved to ~0.5 m for static, coherent targets by considering focused SAR processing [3]. Yet, focused SAR data from delay/Doppler altimeters, like CryoSat-2 and Sentinel-3A/B, suffer from strong azimuth ambiguities as a consequence of the closed-burst-mode operation of the onboard instruments [3]. With the launch of the first interleaved SAR altimeter satellite, Sentinel-6 Michael Freilich, focused SAR processing can now be fully exploited [4].

The resulting enhanced resolution makes altimeter measurements sensitive to long-wave modulations of backscatter (figure 1) [5]–[7]. Often regarded as a nuisance, recent studies have shown that the modulations may be exploited to retrieve swell-wave parameters [6]. [1] first attempted to interpret the backscatter modulations in the waveform tail using a spectral analysis, but a full understanding is still needed. Particularly, the non-linear range bunching was not taken into account and left-right modulation differences remained to be interpreted.

This motivates the present study, to extend the work of [1] to provide a more consistent interpretation of nadir and near-nadir SAR altimetry modulations. The new proposed closed-form model is compared to numerical solution and we discuss several of its properties and validity for nadir SAR altimetry and swath altimetry. The numerical model includes non-zero-Doppler geometries, and principles of sublooking and cross-spectral analysis [8] are considered for SAR altimetry. The long aperture length (~2 seconds) indeed enables to compute a 'stack' of cross-spectra, to provide means to improve ocean-wave-spectra retrievals. For validation, the forward model is compared to Sentinel-6 cross-spectral stacks for a few scenarios. Note that the models used in this paper support the description of the SAR altimetry modulation mechanisms, but would require further adjustments (e.g. curved Earth, etc.) if used for ocean-wave-spectra inversions.

# II. A NUMERICAL MODEL

We first describe a numerical model to compute a SAR altimetry spectrum in a zero-Doppler configuration. A flat Earth and a straight, parallel orbit are assumed. The integration time is assumed to be sufficiently short to neglect backscatter changes during the integration time. A right-handed reference system will be used with the x-, y- and z-axis pointing cross-track, along-track and up, respectively with the origin at the satellite nadir location.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Focused SAR radargram from Sentinel-6 acquisition over the Pacific Ocean on the 5th of January 2021. The along-track posting distance is just over 1 m and the range sampling is 0.19 m. The tails of the waveforms exhibit backscatter modulations induced by long ocean waves.</figcaption></figure>

## A. Scene generation

A typical wavy scene is generated from the spectral domain on an $(N_x, N_y)$ regular grid. It consists of the superposition of a Gaussian swell-wave spectrum and a wind-wave spectrum based on [9], such that

$$S(k_x, k_y) = S_{\text{swell}}(k_x, k_y) + S_{\text{wind}}(k_x, k_y), \quad (1)$$

normalized as

$$\text{SWH} = 4\sqrt{\int S(k_x, k_y)dk_ydk_x}. \quad (2)$$

The maximum wavenumbers $(k_x, k_y)$ correspond to wavelengths of about 5 m to ensure most of the velocity variance is captured. Alternatively, buoy or wave-model spectra can also be considered for model verification with Sentinel-6 data. The unresolved wind waves are in that case also modeled with the Elfouhaily spectrum. From the wave spectrum, the elevation $h$ is computed as

$$h = \mathcal{F}^{-1}\left\{\sqrt{0.5S(k_x, k_y)\Delta k_x \Delta k_y N_x^2 N_y^2} \cdot \exp(i2\pi X) + \text{c.c.}\right\}, \quad (3)$$

where $\mathcal{F}^{-1}$ denotes the inverse discrete Fourier transform, $X$ is a uniformly distributed random variable to introduce random phase offsets, $\Delta k_x, \Delta k_y$ is the two-dimensional spectral resolution and c.c. denotes the complex conjugate. Using the deep-water dispersion relation for the pulsation $\omega = \sqrt{gk}$ of long waves and $k = \sqrt{k_x^2 + k_y^2}$, we compute the vertical velocity field as

$$v = \mathcal{F}^{-1}\left\{-i\omega\sqrt{0.5S(k_x, k_y)dk_xdk_yN_x^2N_y^2} \cdot \exp(i2\pi X) + \text{c.c.}\right\}. \quad (4)$$

Finally, the surface slope in the direction of the satellite, is given by

$$s_x = \mathcal{F}^{-1}\left\{-ik_x\sqrt{0.5S(k_x, k_y)dk_xdk_yN_x^2N_y^2} \cdot \exp(i2\pi X) + \text{c.c.}\right\}. \quad (5)$$

For the slopes, the integral is limited to $k_{lim} = 2\pi/5$ rad m$^{-1}$ as the rest of the surface slopes are inherently captured by the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of input elevation, velocity and slope grids for both sides of the track.</figcaption></figure>

quasi-specular scattering model in section II-B. Note, vertical-velocity contribution drops quickly at wavelengths shorter than 10 m, and error estimates for the correlation between vertical velocity and filtered slope are negligible. The integration time is assumed to be short $T \ll 1$ s so that the viewing geometry can be assumed constant. If multiple sublooks are considered, the surface evolves using phasors of the form $\exp(-i\omega t)$, where $t$ is the time with respect to zero-Doppler. Surfaces are generated on both sides of the track corresponding to the part of the waveform tail to be considered. An example of such a realization for both sides of the track is shown in figure 2.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Topview of the sublook geometry where the arrows indicate the direction of sensitivity. Velocity bunching occurs due to shifts in the along-track direction and is maximal near peaks or troughs. Range bunching occurs due to shifts in the cross-track direction and has its maximum on the slope facing the instrument. Tilt modulation is not a shift, but has its maximum on the slope facing the satellite.</figcaption></figure>

## B. SAR altimetry model

A top-view geometry is given in figure 3. At near-nadir angles the dominant mechanism is quasi-specular scattering. For a rough surface without any long waves, specular backscatter is approximated by [10]

$$ \sigma_{0,\text{ref}} = \frac{1}{2 \cos^4(\theta) s_{aw} s_{cw}} \exp\left(-\frac{\tan^2(\theta)}{2s_c^2}\right), \quad (6) $$

where $s_{aw}$, $s_{cw}$ are the maximum (along wind) and minimum (cross wind) mean-square slopes in perpendicular directions, respectively, while $s_c^2$ is the cross-track mean-square slope, computed as

$$ s_c^2 = \iint k_x^2 S(k_x, k_y) dk_x dk_y. \quad (7) $$

The mean-square slopes are computed using an Elfouhaily spectrum between the wavelengths of ~10 m and 1 cm [11]. For a nadir altimeter the incident angle $\theta$ varies rather rapidly as a function of cross-track distance, and backscattered intensities are not constant over the waveform tail. Local long-wave tilt alters the local incident angle

$$ \theta_l = \arctan(\tan(\theta) + s), \quad (8) $$

so that the NRCS becomes

$$ \sigma_0 = \frac{1}{2 \cos^4(\theta_l) s_{aw} s_{cw}} \exp\left(-\frac{\tan^2(\theta_l)}{2s_c^2}\right). \quad (9) $$

For each of the scatters in a narrow Doppler strip, the range as a function of cross-track distance $x$ and elevation, is

$$ R = \sqrt{(H - h)^2 + x^2}, \quad (10) $$

where $H$ the satellite altitude. The scatter velocity causes a different Doppler compared to a static target, to induce an along-track shift. Under the assumption that the along-track distance over which the scatterer travels is small, the along-track shift is computed as

$$ \Delta y = \frac{R}{V} v, \quad (11) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. A realization of the normalized waveform tail and its associated spectrum based on the input of figure 2.</figcaption></figure>

where $V$ is the platform velocity. Using a two-dimensional sinc-function, the unnormalized intensity in a waveform bin $n$ for Doppler strip $m$ is then computed as

$$ I(n, m) = \sum \text{sinc}^2\left(\pi \frac{R - R_n}{R_{res}}\right) \text{sinc}^2\left(\pi \frac{y + \Delta y - y_m}{y_{res}}\right) \sigma_0. \quad (12) $$

The along-track mean waveform tail is then used for normalization [1]. A realization of the normalized waveform tail and its associated SAR spectrum are shown in figure 4.

## III. A CLOSED-FORM SOLUTION

A closed-form mapping between the ocean and SAR-altimetry spectra could provide interpretation leading to fast inversion algorithms. A closed-form solution will unfortunately quickly loses its validity for wave slopes close to, or larger than the incident angle. However, a closed-form derivation provides insights into the behavior of the measurable nadir spectra. This will also benefit a more elaborated explanation of the detected modulations obtained from non-zero-Doppler sublooks, discussed in section IV-A. Still, a reasonably accurate closed-form solution well applies to interpret swath altimeter measurements under moderate environmental conditions. The latter would lead to similar type of inversion as for Sentinel-1 Wave Mode spectra [12].

### A. SAR spectral transform

Consider a surface that is analytically described by

$$ h(x, y) = \sum_k h_k e^{i(k_x x + k_y y)}, \quad (13) $$

where $h_k$ is the amplitude at wavenumber $(k_x, k_y)$. Then we can map any 'linear' image distortion using a simple transfer function $T(k_x, k_y)$, which results in

$$ \tilde{I}(x, y) = \sum_k T(k_x, k_y) h_k e^{i(k_x x + k_y y)}. \quad (14) $$

With that, it is possible to derive a transfer function for the tilt modulation based on the scattering theory discussed in section II-B, given as [13]

$$ T_I = -ik_x \frac{1}{\sigma_0} \frac{\delta\sigma_0}{\delta\theta}, \quad (15) $$
in line with transfer functions used for spectrometers, like SWIM instrument on-board CFOSAT satellite.

Range bunching and velocity bunching do not alter the scattering properties of the surface. In a similar way as the tilt modulations, the shifts are modeled using transfer functions, such that

$$ 
\begin{aligned}
\Delta x &= \sum_k T_x(k_x, k_y) h_k e^{i(k_x x + k_y y)} \\
\Delta y &= \sum_k T_y(k_x, k_y) h_k e^{i(k_x x + k_y y)},
\end{aligned}
\quad (16) $$

with the transfer functions for cross-track migrations [14]

$$ T_x = -\frac{1}{\tan(\theta)} \quad (17) $$

and the along-track migrations [15]

$$ T_y = -\frac{R}{V}(i\omega) \approx -\frac{H}{V}(i\omega). \quad (18) $$

The derivation of the closed-form transform follows previous developments [8], [15]–[17]. Consider a Fourier transform taken of a SAR image $I_v$

$$ \chi(\vec{k}) = \frac{1}{2\pi} \int I_v e^{-i\vec{k}\vec{x}} d\vec{x}. \quad (19) $$

The signal at SAR image location $\vec{x}$ actually originates from a set of scatters, shifted there from other locations $\vec{x'}$. It can be written as [15]

$$ \chi(\vec{k}) = \frac{1}{2\pi} \int I_r e^{-i\vec{k}(\vec{x} - i\Delta\vec{x'})} d\vec{x'}, \quad (20) $$

where $\Delta\vec{x'}$ represent the cross-track and along-track shifts and $I_r$ is the intensity image after tilt bunching only. The (cross-)spectrum is computed as $\langle \chi_1(\vec{k}) \chi_2^*(\vec{k}) \rangle$. The SAR image spectrum can be written

$$ P(k_x, k_y) = \frac{1}{(2\pi)^2} \iint G(x, y, k_x, k_y) e^{-i(k_x x + k_y y)} dxdy, \quad (21) $$

which is a slightly modified version of equation 6 in [8], where the integral should be computed for each $(k_x, k_y)$ separately. Following Gaussian statistics, the G-function is described with a set of cross-correlations as

$$ 
\begin{aligned}
G(k_x, k_y, x, y) = e^{k_x^2 \mu_{xx} + k_y^2 \mu_{yy} + k_x k_y (\mu_{xy} + \mu_{yx})} \\
& \quad (1 + \rho_{II} + ik_x(\rho_{Ix} - \rho_{xI}) + ik_y(\rho_{Iy} - \rho_{yI}) \\
& \quad + ik_x^2 \mu_{Ix} \mu_{xI} + ik_y^2 \mu_{Iy} \mu_{yI} + ik_x k_y (\mu_{Ix} \mu_{yI} \\
& \quad + \mu_{Iy} \mu_{xI})),
\end{aligned}
\quad (22) $$

with

$$ \mu_{ab}(x, y) = \rho_{ab}(x, y) - \rho_{ab}(0, 0). \quad (23) $$

For clarity, the $(x, y)$ dependence of the cross-correlations is left out. The complete derivation of the G-function can be found in [8]. The cross-correlations are computed in the spectral domain as

$$ 
\begin{aligned}
\rho_{ab} = & \frac{1}{(2\pi)^2} \iint \left[ \frac{1}{2} T_a(k_x, k_y) T_b^*(k_x, k_y) S(k_x, k_y) + \right. \\
& \left. \frac{1}{2} (T_a(-k_x, -k_y) T_b^*(-k_x, -k_y))^* S(-k_x, -k_y) \right] e^{i(k_x x + k_y y)} dk_x dk_y,
\end{aligned}
\quad (24) $$

with $S$ the wave spectrum and $T_{a,b}$ are transfer functions. Higher-order terms are ignored leading to

$$ G(k_x, k_y, x, y) = e^{k_x^2 \mu_{xx} + k_y^2 \mu_{yy} + k_x k_y (\mu_{xy} + \mu_{yx})} (1 + \rho_{II}), \quad (25) $$

and the SAR altimetry spectrum

$$ 
\begin{aligned}
P(k_x, k_y) = & \frac{1}{(2\pi)^2} \\
& e^{-k_x^2 \rho_{xx}(0,0) - k_y^2 \rho_{yy}(0,0) - k_x k_y (\rho_{xy}(0,0) + \rho_{yx}(0,0))} \\
& \iint e^{k_x^2 \rho_{xx} + k_y^2 \rho_{yy} + k_x k_y (\rho_{xy} + \rho_{yx})} \\
& Le^{-i(k_x x + k_y y)} dxdy,
\end{aligned}
\quad (26) $$

with $L = (1 + \rho_{II})$. Compared to previous efforts [1], this expression explicitly includes non negligible non-linear range bunching effects as cross-track shifts (exponential terms involving $\rho_{xx}$).

## B. Discussion of the closed-form solution

At near-nadir angles, the closed-form solution has certainly more limitations than at higher incident angles, but can still be useful to derive information regarding along- and across-track resolution.

1) **Cut-off wavelengths:** The terms outside the integral of Eq. 26 encode the cross- and along-track resolution of the resulting SAR spectra. The term with zero-lag azimuth auto-correlation $\rho_{yy}(0,0)$ corresponds to the azimuth fall-off function [18]. Using the transfer function for the along-track migrations, it is related to the along-track cut-off wavelength

$$ \lambda_c \propto \pi \sqrt{\rho_{yy}(0,0)} = \pi \frac{R}{V} \sqrt{\sigma_v^2} \quad (27) $$

and depends on the vertical velocity variance $\sigma_v^2$ of detected wave systems [19]. Note, the equation omits correlations between velocities and slopes, acting to reduce the velocity variance detected by the altimeter [20]. $\sigma_v^2$ shall thus be interpreted as a weighted velocity variance. Besides that, short-time life of scatters, for example, wave breakers, is ignored [13], possibly leading to an underestimation of the cutoff at high sea states. The relatively high altitude of Sentinel-6 (1300 km) causes an increased resolution loss with respect to that of other SAR altimeters (CryoSat-2, Sentinel-3). At a typical velocity variance of 0.4 m$^2$s$^{-1}$, the along-track cut-off wavelength for Sentinel-6 is > 200 m, nearly twice larger than for CryoSat-2 and Sentinel-3.

The term with zero-lag auto-correlation $\rho_{xx}(0,0)$ describes the cross-track fall-off function. With

$$ \rho_{xx}(0,0) = \frac{\sigma_e^2}{\tan^2(\theta)} \quad (28) $$

a function of the elevation variance $\sigma_e^2$. The cross-track resolution is therefore proportional to the SWH, such that

$$ \lambda_{ct} \propto \pi \sqrt{\rho_{xx}(0,0)} = \pi \sqrt{\frac{\sigma_e^2}{\tan^2(\theta)}} = \pi \frac{\text{SWH}}{4 \tan(\theta)}. \quad (29) $$

Note that near nadir incident angle $\theta$ quickly changes in the waveform tail, so it cannot be considered constant. At
a incident angle of 0.4° and a SWH of 2 m, the cross-track cut-off is already > 200 m. Swell-wave spectral retrieval will benefit from a large measurement window that covers the full beam-limited footprint of Sentinel-6.

2) Nadir altimeter: We consider an altimeter altitude of 800 km, a chirp bandwidth of 320 MHz and a synthetic aperture length of 500 m. For both the numerical model and the closed-form solution, we estimate the expected value for a spectrum without noise generated from two square scenes between 4000 m and 7000 m on both sides of the ground track. For the closed-form solution, we use an incident angle corresponding to 5500 m across track in a flat-Earth geometry. It is assumed that the satellite does not roll and is located at zero-Doppler. For the numerical model, an average of ten realizations is taken.

Figure 5 shows a comparison of the closed-form and numerical model for swell waves propagating in three different directions. For clarity, we only model the right side of the altimeter ground track, which yields only two ambiguities. There is a moderate wind system present which reduces the resolution in the cross-track and along-track directions. The reduced resolution eliminates any signals from the wind waves and therefore only swell waves modulate the intensity with signals visible in the spectra. The power spectral density of the swell signal changes with direction because of the bunching is dominant across track while velocity bunching is dominant along track [1]. The discrepancies between the closed-form and numerical models also depend on the wave direction. In the cross-track direction, the power spectral density for this particular scenario differs by almost an order of magnitude, while in the along-track direction, both models agree quite well.

Among the processes responsible for these discrepancies, two are related to the non-linear range migration. To discuss this in detail, we first expand the range-migration term into a Taylor series

$$ \exp(k_x^2 \rho_{xx}) = 1 + \frac{k_x^2 \rho_{xx}}{1!} + \frac{k_x^4 \rho_{xx}^2}{2!} + \dots \quad (30) $$

Consider a monochromatic swell system $\eta = A \sin(k'_x y)$, it results in a range-migration cross-covariance function

$$ \rho'_{xx} = \frac{A^2 \cos(k'_x x)}{2 \tan^2(\theta)}. \quad (31) $$

Inserting it into the Taylor expansion,

$$ \exp(k_x^2 \frac{A^2 \cos(k'_x x)}{2 \tan^2(\theta)}) \approx 1 + \frac{k_x^2}{\tan^2(\theta)} \frac{A^2 \cos(k'_x x)}{2} + \frac{k_x^4}{\tan^4(\theta)} \frac{A^4 \cos^2(k'_x x)}{8} + \dots \quad (32) $$

the monochromatic swell system appears to cause responses at $k_x = Nk'_x$, where N is a (positive) integer. This results in peaky modulations, the orange lines in figure 6. For a constant incident angle and low waves, a numerical one-dimensional model well reproduces the closed-form solution. Yet, the incident angle rapidly changes with increasing cross-track distance and the closed-form solution is not accurate.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Numerical (top) and closed-form (bottom) nadir-altimeter SAR spectra simulations for three swell-wave spectra with mean propagation directions of 0°, 45° and 90° from left to right, respectively. The unresolved wind waves are based on a 10 m/s stress-equivalent wind speed with a fetch of 200 km.</figcaption></figure>

To illustrate this effect, the auto-correlation function for a monochromatic wave is approximated

$$ \rho'_{xx} = \frac{H^2 A^2 \cos(k'_x x)}{2x^2}, \quad (33) $$

where H is the satellite altitude and x the cross-track distance. This auto-correlation function corresponds to the normalized range bunching intensity modulation for a monochromatic wave [15]

$$ I_{rb} = \exp(k_x' \frac{HA \cos(k'_x x)}{x}) \quad (34) $$

as shown by the pink line in panel-a of figure 6. The closed-form model (orange line) does not explicitly includes the incident-angle dependence, e.g. the ratio $(H/x)$ is set constant. Therefore the closed-form understandingly deviates from the numerical model, while the model from equation 34 more closely matches. A computationally efficient implementation of the incident-angle dependency in a closed-form solution is not trivial. The varying incident angle does not only cause a trend in modulation directly, but also indirectly a varying cut-off with the incident angle.

While the incident angle variation introduces a rather small discrepancy, the second issue is much more problematic. When $k_x A > \tan(\theta)$ the waves become supercritical (the slope becomes larger than the incident angle), which results in two maxima (panel b). The crest and trough switch range location at one fourth of the wavelength. These maxima are normally not separable as the limited cross-track resolution smears the response. For (near-)supercritical waves, the spectral response changes from the analytical model and the closed-form solution are no longer valid. In the case of a 400 m cross-track swell, supercriticality is already reached at a fraction of a meter. Supercritical velocities can also occur for sea ice [21], which results in comparable patterns in the along-track direction. For the Doppler migration term a similar Taylor expansion can be used to illustrate that supercritical velocities correspond to $k_y A > \frac{V}{\omega R}$ and that upwind and
<figure><img src="image_7.png" /><figcaption>Fig. 6. A comparison of two one-dimensional analytical models with a numerical model at steep incident angles for very low swell-wave heights (top). A comparison between an analytical model with a numerical model for moderate wave height (bottom). Oversampling is applied for interpretation.</figcaption></figure>

downwind slope switch position if the Doppler shift becomes one-fourth of the swell wavelength [22]. Note that the swell radial velocities in the examples of figure 5 are not large enough to be supercritical.

More subtle, other differences are present between the numerical model and the closed-form solution. Given the instrument sampling rate, the normalization and the resampling of the waveforms can induce artefacts that should be carefully considered. However, the non-linear range (and velocity) bunching terms already demonstrate that, unless the swell-wave height is very low, the closed-form solution does not provide an accurate representation of the SAR altimetry spectrum.

3) *Swath altimeter*: Based on the preliminary design of Sentinel-3 Next Generation, a swath altimeter flying at an altitude of 800 km with a chirp bandwidth of 200 MHz and a synthetic aperture length of 500 m is considered. We estimate the expected values of the signal spectrum for a small scene between 52000 m and 56000 m across track with incident angles close to approximately 4 degrees. The incident angle for the closed-form model is kept constant and corresponds to 54000 m across track in a flat-Earth geometry.

The closed-form and numerical models for swell and wind waves traveling in three different directions are shown in figure 7. Compared to the SAR spectra of the nadir altimeter, wind-wave system starts to be visible with the high resolution in the cross-track direction. Although the cross-track resolution of the swath altimeter is not as high as a side-looking SAR system, like Sentinel-1, the modulation due to range bunching is much stronger which allows to have a better constraint on the wind-wave system. A full inversion of the two-dimensional wind-wave system, without external information, is still impossible with the clearly visible cut-off in the along-track direction. The discrepancies of the swell signals between the closed-form and numerical spectra are smaller than for the nadir altimeter, even though the slopes considered here are higher. Under moderate swell heights, typically lower than

<figure><img src="image_8.png" /><figcaption>Fig. 7. Numerical (top) and closed-form (bottom) swath-altimeter SAR spectra simulations for three swell-wave spectra with mean propagation directions of 0°, 45° and 90° from left to right, respectively. The unresolved wind waves are based on a 10 m/s stress-equivalent wind speed with a fetch of 200 km.</figcaption></figure>

4 m, the supercriticality threshold $k_x A > \tan(\theta)$ is not exceeded. Neither does the change of incident angle over the considered 3 km wide scene play a significant role (figure 8).

For low (swell-)wave heights, simplifications might be made using the first-order Taylor expansion (see figure 8) of the cross-track migration term

$$e^{k_x^2 \rho_{xx}} \approx 1 + \frac{k_x^2 \rho_{xx}}{1!}, \qquad (35)$$

with which we simplify equation 26 to

$$P(k_x, k_y) = \frac{1}{(2\pi)^2} e^{-k_y^2 \rho_{yy}(0,0)} \iint e^{k_y^2 \rho_{yy}} L' e^{-i(k_x x + k_y y)} dxdy, \qquad (36)$$

where $L' = 1 + \rho_{II} + \rho'_{xx} + ... \approx 1 + \rho_{II} + \rho'_{xx} = 1 + \rho'_{II}$ and $\rho'_{xx} = k_x^2 \rho_{xx}$. This form of the closed-form equation matches the approximation used in side-looking SAR [17], [23], and the one given by [1]. It inherently assumes that the cross-track resolution loss is negligible and/or that secondary peaks are not present. Considering a linear cross-track migration term substantially reduces the computation time and simplifies the inversion (see for the derivation in [8]) at the cost of limited applicability and reduced performance. Alternatively, a second- or third-order approach might be performed. In principle, this would allow the retrieval to follow a comparable path as the Sentinel-1 ocean-wave retrieval algorithm [12].

The processing of a swath-altimeter SAR spectrum involves either a Welch or Bartlett periodogram averaging. Averaging of periodograms for Sentinel-3 Next Generation Topography is done in the outer swath. If we assume this swath covers 48 km to 60 km across track, it corresponds to incident angles between 3.4 to 4.3 degrees in a flat-Earth geometry. A forward modeled spectrum using the mean incident angle will therefore differ from the observed average as $\tan(\theta)$ changes nonlinearly. Improved performance is expected when the same
<figure><img src="image_9.png" /><figcaption>Fig. 8. A comparison of two one-dimensional analytical models with a numerical model at steep incident angles for a swell system only and a swell system in the presence of a wind-wave system.</figcaption></figure>

averaging is taken into account in the model, which comes at the cost of a higher computation time.

## IV. SUBLOOK PROCESSING

We now focus on the use of numerical simulations for the nadir altimeter for which closed-form solutions are normally not valid. Nadir radar altimeter overpasses allow for a large integration time compared to side-looking SAR systems and swath altimeters. Instead of integrating over the full aperture, it is possible to split the aperture into sublooks and compute multiple independent single-look waveforms corresponding to the same location. By combining the DFTs of the normalized waveform tails of different sublooks, cross-spectra are computed. In side-looking SAR, sublooking and its associated cross-spectral analysis are exploited to both remove the speckle noise pedestal [8], [24] and a 180-degree wave propagation ambiguity. It can be anticipated that the phase is not only a function of wave propagation, but also related to changes of the observation geometry. For the nadir altimeter, it will be shown that not only the cross-spectral phase is a function of the observation geometry, but also the cross-spectral density. Furthermore, because of the larger geometrical diversity (due to the long aperture length) of nadir altimeters, a stack of cross-spectra can be computed, which can be exploited to further constrain the ocean-wave spectral inversion procedure.

### A. The non-zero-Doppler geometry model

Suppose we focus on a narrow Doppler strip at an along-track position $y_{ref} = Vt$ with respect to the platform position (figure 9). The local incident angle with a flat surface becomes a function of cross-track distance and the distance to the along-track strip

$$ \theta = \arctan\left(\frac{\sqrt{y_{ref}^2 + x^2}}{H}\right). \qquad (37) $$

The ground-range direction is not perpendicular to the along-track direction anymore. We therefore write the equation for the surface slope along the direction of the satellite as

$$ s = s_x \cos(\phi_r) + s_y \sin(\phi_r), \qquad (38) $$

where $\phi_r$ is the angle between the ground-range direction and the cross-track direction, given as

$$ \phi_r = \arctan\left(-\frac{y_{ref}}{x}\right). \qquad (39) $$

The cross-track and along-track slopes are computed with

$$ \begin{aligned} s_x &= \mathcal{F}^{-1}\{-ik_x \sqrt{0.5S(k_x, k_y)dk_xdk_y N_x^2 N_y^2} \\ & \qquad \exp(i2\pi X) \exp(-i\omega t) + \text{c.c.}\} \\ s_y &= \mathcal{F}^{-1}\{-ik_y \sqrt{0.5S(k_x, k_y)dk_xdk_y N_x^2 N_y^2} \\ & \qquad \exp(i2\pi X) \exp(-i\omega t) + \text{c.c.}\}. \end{aligned} \qquad (40) $$

From these equations, it can easily be shown that $s$ is highest when the waves are traveling in the ground-range direction (perpendicular to the range-isolines in figure 9). Consider a Doppler strip ahead of the nadir location and the altimeter only illuminates the surface on the right side of the ground track. A wave traveling at 0-90 degrees or 180-270 degrees with respect to cross-track, will therefore experience a larger tilt modulation than waves traveling in the other two quadrants. If the altimeter would only illuminate the left side of the ground track this would be reversed. As nadir altimeters illuminate both sides of the ground track, the tilt modulation is expected to be stronger in two quadrants in case a non-zero-Doppler geometry is considered, unless we have along-track travelling waves. For the numerical model, we also adjust the NRCS equation

$$ \sigma_0 = \frac{1}{2 \cos^4(\theta_l) s_{aw} s_{cw}} \exp\left(-\frac{\tan^2(\theta_l)}{2s_r^2}\right), \qquad (41) $$

where $s_r^2$ is the mean-square slope in the radar look direction given as [25]

$$ s_r^2 = \frac{s_{aw}^2 s_{cw}^2}{s_{cw}^2 \sin(\phi_r')^2 + s_{aw}^2 \cos(\phi_r')^2} \qquad (42) $$

and therefore $\phi_r'$ is the ground-projected angle between the radar look direction and the semi-major axis of the surface roughness ellipse. Note that $\sigma_0$ and its derivative depend on both the magnitudes and directions of the long-wave slope (i.e. $\theta_l$) and the mean-square slope.

Range bunching will also enter in a non-zero-Doppler geometry. The along-track distance $y_{ref}$ between the satellite and the target should be considered when projecting the waveforms onto the surface. For a narrow Doppler strip, a rather simple adjustment to the range gives

$$ R = \sqrt{(H-h)^2 + y_{ref}^2 + x^2}. \qquad (43) $$

For a wider Doppler strip, the range varies in the along-track direction which cannot be ignored. This results in

$$ R = \sqrt{(H-h)^2 + (-y_{ref} + dy)^2 + x^2}, \qquad (44) $$
<figure><img src="image_10.png" /><figcaption>Fig. 9. Topview of a forward-looking geometry (sublook 1) and a backward looking geometry (sublook N). Elevation and Doppler shifts occur in the direction of the arrows. The peak of the range bunching modulation is therefore located on the wave slope facing the ground track. Doppler shifts are not perpendicular to the elevation shifts anymore. Velocity bunching still occurs near the crests and troughs of the wave, but is maximal when waves are propagating perpendicular to the ground-range direction. Tilt modulation is maximal on the slopes facing the satellite. The blue squares indicate the illuminated area corresponding to the waveform tail.</figcaption></figure>

where *dy* is the difference in along-track direction between the scatter and the reference range. Note that *dy* is considered to be positive in the flight direction and *y*<sub>ref</sub> is considered negative for a forward-looking geometry. By first not considering any velocity bunching, it is possible to infer the range bunching effects. Consider a rather wide Doppler strip and a monochromatic wave with low amplitude. A crest at the furthest along-track location in the Doppler strip can compensate for the increased range of a positive *dy*. Therefore, scatters from troughs and crests might be mapped near each other, which would reduce the range bunching modulation. As with the tilt, this causes a left-right difference. Ignoring secondary effects, the maximum range bunching occurs when *φ*<sub>k</sub> − *φ*<sub>r</sub> ≈ 0. A wide Doppler strip therefore acts as a filter that depends on the wavelength and direction of the waves, the width of the strip, and the distance between zero-Doppler and the strip. Note also that a scatter at the edge of two Doppler strips gets projected at a different range in both waveforms, as in one Doppler strip and has along-track location *y*<sub>ref</sub> + *dy* and in the other *y*<sub>ref</sub> − *dy*. For the numerical simulation the ranges therefore also have to be recomputed for each Doppler strip.

The most prominent change of modulation occurs due to Doppler shifts and the associated velocity bunching. Doppler shifts do not cause scatters to move in the along-track di-

<caption>TABLE I</caption>

RELATIVE SIZE OF THE MODULATIONS FOR A MONOCHROMATIC WAVE
TRAVELING BETWEEN 0-90 DEGREES WITH RESPECT TO CROSS TRACK.
RANGE AND VELOCITY BUNCHING ARE ABBREVIATED AS RB AND VB,
RESPECTIVELY.

<table>
  <thead>
    <tr>
      <th>Look</th>
      <th>Side</th>
      <th>Tilt</th>
      <th>RB</th>
      <th>VB</th>
      <th>Cut-off rot.</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Forward</td>
      <td>Right</td>
      <td>+</td>
      <td>+</td>
      <td>-</td>
      <td>+</td>
    </tr>
    <tr>
      <td>Forward</td>
      <td>Left</td>
      <td>-</td>
      <td>-</td>
      <td>+</td>
      <td>-</td>
    </tr>
    <tr>
      <td>Backward</td>
      <td>Right</td>
      <td>-</td>
      <td>-</td>
      <td>+</td>
      <td>-</td>
    </tr>
    <tr>
      <td>Backward</td>
      <td>Left</td>
      <td>+</td>
      <td>+</td>
      <td>-</td>
      <td>+</td>
    </tr>
  </tbody>
</table>

rection alone in a non-zero-Doppler geometry. As shown in Kleinherenbrink et al. (to be submitted), line-of-sight-moving scatters experience a shift approximately along range isolines. For nadir altimeters, the range isolines are arcs on the surface (figure 9). At a sufficient distance from nadir, these arcs can locally be approximated by a straight line with the direction

$$
\phi_{sh} = \phi_r + \frac{\pi}{2} \qquad (45)
$$

for the right side of the ground track. The scatterer still
requires the same along-track Doppler shift Δ*y*<sub>D</sub> as a con-
sequence of its own geophysical Doppler and therefore the
total shift will be larger than for a zero-Doppler geometry.
The cross-track shift is related to the along-track shift as

$$
\Delta x_D \approx \frac{\Delta y_D}{\tan(\phi_{sh})} = \frac{R}{V \tan(\phi_{sh})} v. \quad (46)
$$

In fact, the negative cross-track shift ensures the increased
range due to the positive along-track shift to be compensated.
The inclined shift also introduces a left-right difference. Con-
sider again a monochromatic wave propagating in the ground-
range direction φ<sub>r</sub> on the right side of the ground track.
Scatters with vertical velocities are moved in the direction φ<sub>sh</sub>,
which is aligned with the wave slopes, so no velocity bunching
will occur. On the left side of the track, scatters are moved in
the direction π − φ<sub>sh</sub> = π/2 − φ<sub>r</sub>. This is perpendicular to the
wave slopes and therefore maximum velocity bunching will
occur. The size of the modulations for the left and right sides
of the track are summarized in table I, where the plus and
minus signs indicate relatively strong and weak modulations.

As the Doppler shift becomes two-dimensional in a non-
zero-Doppler geometry, it will also contribute to the cross-
track cut-off wavelength. This will introduce a rotation of the
fall-off ellipsoid. Using the function for the cross-track shift,
we can adapt the cross-track transfer functions to include a
velocity bunching term

$$
T_x = - \frac{1}{\tan(\theta)} - \frac{H}{V \tan(\phi_{sh})} (i\omega). \quad (47)
$$

Inserting this into the cross-correlation functions
leads to a non-zero contribution of the cross-term
e⁻⁻ᵏₓkᵧ(ρₓᵧ(0,0) + ρᵧₓ(0,0)) as ρₓᵧ(0,0) + ρᵧₓ(0,0)) does
not cancel (see appendix). In the presence of wind waves, the
variation of the PSD in the spectral stack in dominated by
the cut-off rotation as we will show in the following section.

B. Cross-spectral stack model results

Sublooks are generated by assuming a quasi-static geometry,
i.e. within the sublook linear Doppler shifts are considered and
the geometry is kept constant during the 500 m integration length. The realizations of the intensities from both sides of the nadir track are linearly added together and consecutively normalized and resampled. A discrete Fourier transform is applied to the resampled normalized waveform tails of the sublooks after which a cross-spectrum is computed from two sublooks separated by 1000 m between their respective centers. Eventually, nine cross-spectra are computed and stored in a cross-spectral stack. For the discussion in this section, we assume a swell spectra with a mean propagation directly of 45° as a reference.

1) *Swell-wave modulations*: To isolate the various effects discussed in the previous section, we first consider a surface with only a swell-wave system and homogeneous roughness. In this section, we limit the discussion to swell traveling in the quadrants and omit the discussion of along-track and cross-track travelling waves. Figure 10 shows the absolute value of spectral stacks computed from a time-varying ocean surface and observation geometry for the modulations due to tilt, range shifts, and velocity shifts separately. Any wind or wave direction mentioned in the coming section will follow the Cartesian convention: positive counterclockwise from the x-axis, which is across-track. To ease the discussion, we will designate quadrants 1, 2, 3, and 4 to be the upper-left, upper-right, lower-left, and lower-right quadrants, respectively.

## Tilt modulation

Modulations due to tilt are typically orders of magnitude weaker than those from the elevation and Doppler shifts. For the results in panel a, we assumed a very low surface roughness, which causes the Jacobian of the specular scattering to be relatively large. Under mean sea states (SWH = 2 m, $U_{10}$ = 8 m/s), the tilt modulation is even weaker than shown in the top panels of figure 10. Note that the relative change in modulation is rather strong. As expected, in the zero-Doppler geometry (center column) four ambiguities are present with equal power-spectral density (PSD), a consequence of the homogeneous roughness. The signal in Q1 stems from the left side of the mirroring of the signals on the left side of the track, while the signals in Q2 are the signal arriving from the right side of the track. The signals in Q3 and Q4 are their respective 180-degree ambiguities. In the forward-looking geometry (negative satellite y-locations) the slopes of swell propagating at a direction of 45° are facing the satellite if they are on the right side of the track. On the left side of the ground track the swell crests are nearly aligned with the ground range. Tilt modulations from the right side of the track are therefore stronger and only a signal in Q2 is visible. At a backward-looking geometry (positive satellite y-locations) the modulations originating from the left side of the ground track are stronger. Figure 11 shows how the peak PSD evolves in slow time (during the overpass). As visible, the histories are mirrored around zero-Doppler. The signals from the right side of the track (Q2 and its 180-degree ambiguity Q4) decrease in slow time, while the other two increase. Due to the relatively weak tilt modulation, its effect on the overall modulation history can be ignored under mean sea states.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Absolute values of numerically modeled cross-spectra for five geometries considering modulation mechanisms separately. The input is a Gaussian swell-wave system with a mean propagation direction of 45°. Top: spectral response of tilt modulation at 17 m along-track resolution. Second row: spectral response of velocity bunching at 17 m 'along-track resolution'. Third row: spectral response of range bunching at 17 m along-track resolution. Last row: spectral response of range bunching at 88 m along-track resolution.</figcaption></figure>

## Velocity bunching

Far stronger modulations are observed as a consequence of velocity bunching (figure 10). Its history is opposite of that of tilt modulation, i.e. for the considered swell-wave direction the velocity bunching is strongest in Q1 and Q4 for the forwards looks. At -6000 m, the mean ground-range direction in the considered ocean surface is nearly 45°, which corresponds to the swell propagation direction. Remember that the maximum velocity bunching occurs when swell is propagating perpendicular to the ground-range direction, which is the case for the left side of the track. As such, in the forward-looking geometry, the maximum PSDs are found in Q1 and Q4.

## Range bunching

Range bunching is the strongest modulational effect (figure 10), but has a relatively weak variation in slow time when processed at high resolution (figure 11). As discussed in the section IV-A, the slow-time variation of range bunching is enhanced when the data is processed at a lower resolution. In the presence of wind waves, Doppler shifts also reduce the (along-track) resolution, which can result in similar behavior. However, as will be discussed later, wind waves introduce additional complexity so that the interpretation is not that trivial. Range bunching is largest when the alignment of the wave direction and the line-of-sight is closest and therefore signals from the right side of the track (Q2 and its 180-degree ambiguity Q3) are stronger in the forward-looking geometry. The dominance of range bunching is contested by velocity bunching [1] for particular wave directions and observation geometries.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Peak PSD variation in each of the quadrants as a function of the along-track position for each of the modulation mechanisms separately. The input is a Gaussian swell-wave system with a mean propagation direction of 45°. A 17 m along-track resolution is considered for the first three panels. The last panels show range bunching at an along-track resolution of 88 m.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13. Numerical model of cross-spectra (absolute values) for five different geometries and four different swell-wave directions in the presence of a fully-developed wind-wave system at a stress-equivalent wind of 8 m s⁻¹.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Top rows: absolute values of numerically modeled cross-spectra for five geometries considering all modulation mechanisms for swell-wave systems traveling at 45° and 135°. In the bottom two panels their associated peak PSD variations on the left and right, respectively.</figcaption></figure>

## Full response

The spectral response of the combined modulations is shown in figure 12 for two swell-wave directions. As the modulations are not in-phase with each other and their (relative) amplitudes and phases (near the crest, trough, or maximum slope) depend on the observation geometry, the interpretation of a combined signal is not always straightforward. However, because the tilt modulation is an order of magnitude smaller it can safely

be ignored. The integral effect of all processes leads to a modulation that resembles the pattern of range bunching, independent of the observations geometry (first row). Its associated peak PSD variation (bottom left panel) is even weaker than that of range bunching at higher resolution. This suggests the opposing behavior of velocity and range bunching counter each other. The combined effects lead to an along-track variation of the peak modulation of only 5% for this particular set of conditions. Depending on the direction of the swell it might be larger as the modulation variations of either velocity or range bunching will start to dominate.

The second row shows the spectral response of a similar swell system propagating at 135°. Its associated peak PSD variations are shown in the bottom-right panel. There is nearly no difference between the spectral response of a swell-wave system propagating at 45° or 135° as the signals get mirrored over the ground-track. Any minor deviations arise from differences in the realizations of the surface. However, the signals in Q2/Q3 originate from the left side of the track instead of the right. In terms of the absolute spectral response, there is also no difference between a swell system propagating at 45° and 225° or 315°. The peak velocity bunching will occur shift 180-degree in phase for linear waves, but the magnitude and phase difference with respect to range bunching will remain the same. So even if modulations from the left and the right side of the track have different peak PSD variations, the absolute values of a cross-spectral stack will not solve ambiguities.

2) *The effect of wind waves:* In theory, the line-of-sight mean-square slope depends on the wind direction. The mean-square slope determines the backscatter and can therefore increase or reduce the relative 'weight' between the left and right side of the track. It can easily be shown that this effect is rather small for the considered incident angles. The Jacobian
of the backscatter is also altered by wind direction, which would introduce a tilt modulation that is dependent on the alignment of swell and the short wind waves. If the wind direction is known, it would allow to removal of the left-right ambiguity in the cross spectra. However, as the tilt modulation is orders smaller than the velocity bunching and range bunching this effect will not be detectable.

In figure 13 cross-spectra are shown for the same swell-wave system as discussed before, but superpositioned with a fully-developed wind-wave system. Wind waves are responsible for a large fraction of the wave height and the velocity variance, which reduces the cross-track and along-track resolution (see section III-B1). As discussed in section IV-A, the two-dimensional falloff ellipse changes orientation as a result of the cross-track shifts introduced by the radial velocity variance. The first row of cross-spectra shows that the cutoff indeed appears to rotate. The swell-system propagating at 45° has strong signals in Q2 and Q3 (originating from the right side of the ground track) for the forward-looking geometry, while the signals are strongest in quadrant Q1 and Q4 (originating from the left side of the ground track) is the backward-looking geometry. The imbalance in spectral energy between the four quadrant, described in [1], are therefore likely the result of a non-zero-Doppler processing geometry. Note that the spectral density can be even enhanced with respect to the zero-Doppler geometry. In part, the induced range-bunching variation can be interpreted as coming from the loss of the (along-track) resolution due to velocity bunching (discussed in section IV-B1). Whether the swell is propagating at 45°, 135°, 225° or 315°, the modulation behavior is nearly identical. Also with wind waves present, the four ambiguities cannot be resolved using the absolute or real part of the cross-spectra. A discussion on the imaginary part of the cross-spectra is provided in the supplementary material. Due to a combination of model incompleteness and noise, the imaginary parts of the cross-spectra remain difficult to interpret.

### C. Comparison with Sentinel-6

Our model is compared with SAR spectra computed from the tails of focused Sentinel-6 waveforms. Input for our Sentinel-6 data processing are level-1a I/Q samples. Using the Delft AltimeteR Toolbox (DART) the waveforms are focused [26]–[28]. We ensure that the waveforms align (figure 1) to a fraction of the (oversampled) range bin size of 19 cm using a range-migration correction based on a smooth surface, which is in this case the geoid. An iterative approach might be applied, where the focused waveforms are retracked and the resulting sea-surface height is along-track low-pass-filtered to serve as a reference surface, although only marginal differences are observed in open-ocean observations at moderate conditions.

To avoid storing multiple waveforms for each ground-track location, we apply a procedure like in SAR wave mode data, where the Doppler spectrum is computed from a single-look complex image in order to generate the sublooks [29]. The (phase) alignment of the leading edge is critical for this procedure and therefore we avoid large tracker jumps. It is also

not possible to keep the RCMC constant over more than ~10 waveforms as it will introduce discontinuities in the sublooked radagrams. The Doppler spectrum is computed with

$$D(x, ky) = \mathcal{F}_y\{I(x, y)\}, \qquad (48)$$

where $\mathcal{F}_y$ is the Fourier transform in the along-track direction and $Ix, y$ is the single-look complex radagram. The Doppler spectrum is then split using a windowing function $W_i(k_y)$. The window has width $N_w/N_s$, where $N_w$ is the number of waveforms to be considered and $N_s$ is the number of sublooks we intend to create. Eventually, the radagram for sublook $i$ is computed as

$$I_i(x, y) = \mathcal{F}_y^{-1}\{D(x, k_y)W_i(k_y)\}. \qquad (49)$$

Three sublooked waveforms are shown in figure 14. If the range history is nearly constant over the considered region, this is comparable to splitting the aperture before focusing. As visible in the figure, the sensitivity to backscatter modulations is look-dependent. In the forward look the curves bend with increasing along-track distance towards the leading edge and for the backward look vice versa. After that, the normalization is performed consistent with the numerical model. Note that we take the curvature of the Earth into account for the projection of the waveform samples onto the surface.

The forward model is estimated using colocated ERA5 WAM spectra supplemented by a parametrised short-wave spectrum [9] based on ERA5 winds. Note that ERA5 winds and WAM spectra are based on two-hour average outputs from the model. The WAM-based spectra will deviate from the true ocean-wave spectra, which will prohibit validation in terms of absolute PSD. The cross-spectral stacks for three different scenarios are compared with the forward modeled spectra based on WAM.

Figure 15 shows a scenario where the swell system is propagating in one of the quadrants, which closely compares to the previously discussed scenarios. There is an offset in terms of peak PSD, but the behavior in slow time is well captured. The peak intensities in all four quadrants closely equal to the zero-Doppler geometry. In the forward-looking geometry, the PSD in Q1 and Q4 is decreased while the PSD in Q2 and Q3 is increased, while for the backward-looking geometry, this behavior reverses. In the case of near-cross-track (appendix B, figure 16) and near-along-track (appendix B, figure 17) the spectral signals from left and right start to overlap, yielding only two spectral peaks. The overlap of signals from both sides of the tracks leads to interference. Signals from either side of the quadrants might be suppressed slightly by the looking geometry, albeit with limited sensitivity. The retrieval of swell-spectral parameters in the latter case will therefore anyhow be prone to additional noise from interference.

## V. CONCLUSIONS

This paper provides an in-depth analysis of SAR altimeter waveform tail modulation and its associated cross-spectral analysis. It describes the implementation of a fast numerical model. Using this numerical model, the closed-form solution, generally used for side-looking SAR spectra into ocean-wave
<figure><img src="image_15.png" /><figcaption>Fig. 14. Three radargrams obtained from splitting the Doppler spectrum of a Sentinel-6 aperture into seven sublooks. The top radargram is rendered from the first 1/7th part of 1.7 second total integration time and corresponds to a forward look, where the satellite is located 3.4 km before zero-Doppler. The bottom radargram corresponds to a backward look at 3.4 km from zero-Doppler. The along-track posting rate is kept at ~1 m.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Real part of five cross-spectra estimated from Sentinel-6 data compared to a model. The five cross-spectra are computed from seven sublooks of a 1.7 s aperture, so that the separation between the centers is two sublooks, or about 0.4 s.</figcaption></figure>

spectra, is demonstrated to be inaccurate at steep incident angles. In contrast, SAR spectra from swath altimeters might be well interpreted by closed-form solutions for low to moderate sea states. Even though not fully accurate, closed-form solutions can still be used to derive first-order properties. In particular, it is confirmed that along-track and cross-track resolutions are proportional to the velocity variance and the

SWH, respectively.

This paper also introduces the concept of the cross-spectral stack, in which cross-spectra are formed from multiple sublooks of the seconds-long aperture of the altimeters. The cross-spectral stack enables improved analysis of the swell spectrum, thanks to a varying observation geometry. These changes indeed help infer the directional cut-off wavelength and the modulation, and to also increase the signal-to-noise of the long-wave signals. The non-zero Doppler left-right differences also provide an explanation for the uneven energy distribution in the four quadrants observed in [1]. The imaginary part, or phase, of the cross-spectra may then provide additional information on the wave propagation. The imaginary parts of the SAR altimetry cross-spectral stack, however, do not fully represent the ones computed from Sentinel-6 data. Still, the present model is useful to guide the interpretation and analysis of SAR altimetry waveforms under moderate conditions. It also provides a framework for modeling spectra from SAR altimetry waveform tails. It can be modified and extended to involve additional processes, like non-linear waves [30] or a more realistic geometry.

Quite importantly, the present development opens to more systematically consider cross-spectral stacks to help more accurately retrieve ocean-wave spectral information for certain combinations of wavelength and propagation direction. Re-
trieved ocean-wave parameters shall then contribute to better monitor global wave field. It shall also serve to advance more consistent sea-state-bias corrections.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Real part of five cross-spectra estimated from Sentinel-6 data compared to a model. The mean swell propagation direction is cross-track. The five cross-spectra are computed from seven sublooks of a 1.7 s synthetic aperture, so that the separation between the centers is two sublooks, or about 0.4 s.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Imaginary part of five cross-spectra estimated from Sentinel-6 data compared to a model. The mean swell propagation direction is cross-track. The five cross-spectra are computed from seven sublooks of a 1.7 s synthetic aperture, so that the separation between the centers is two sublooks, or about 0.4 s.</figcaption></figure>

## APPENDIX A

The cross-correlation terms in equation 26 for a zero-Doppler observation are written as a function of transfer functions (see equation 24), such that

$$ 
\begin{aligned}
\rho_{xy} &= f(T_x(k_x, k_y)T_y^*(k_x, k_y)) = f\left(-\frac{H}{V}\frac{i\omega}{\tan(\theta)}\right) \\
\rho_{yx} &= f(T_y(k_x, k_y)T_x^*(k_x, k_y)) = f\left(\frac{H}{V}\frac{i\omega}{\tan(\theta)}\right),
\end{aligned}
\quad (50) $$

where we have inserted the transfer functions from section III-A. This shows that the sum of the two cross-correlation terms becomes

$$ \rho_{xy}(0,0) + \rho_{yx}(0,0) = 0 \quad (51) $$

and therefore cancel each other.

For the non-zero-Doppler geometry the transfer functions from section IV-A, which result in the cross-correlations

$$ 
\begin{aligned}
\rho_{xy} &= f(T_x(k_x, k_y)T_y^*(k_x, k_y)) \\
&= f\left(\frac{H^2}{V^2}\frac{\omega^2}{\tan(\phi_{sh})} - \frac{H}{V}\frac{i\omega}{\tan(\theta)}\right) \\
\rho_{yx} &= f(T_y(k_x, k_y)T_x^*(k_x, k_y)) \\
&= f\left(\frac{H^2}{V^2}\frac{\omega^2}{\tan(\phi_{sh})} + \frac{H}{V}\frac{i\omega}{\tan(\theta)}\right).
\end{aligned}
\quad (52) $$

By separating the Fourier Transform integrals for the imaginary and the real part, it can be shown that the resulting sum of the cross-correlations results in

$$ \rho_{xy}(0,0) + \rho_{yx}(0,0) = 2\frac{H^2}{V^2}\frac{\sigma_v^2}{\tan(\phi_{sh})}, \quad (53) $$

where $\sigma_v^2$ is the velocity variance and $\phi_{sh}$ is the angular dependence as defined in section IV-A.

## REFERENCES

[1] O. Altparmaki, M. Kleinherenbrink, M. Naeije, C. Slobbe, and P. Visser, "SAR Altimetry Data as a New Source for Swell Monitoring," *Geophysical Research Letters*, vol. 49, no. 7, p. e2021GL096224, Apr. 2022. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021GL096224

[2] R. Raney, "The delay/Doppler radar altimeter," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 36, no. 5, pp. 1578-1588, Sep. 1998. [Online]. Available: http://ieeexplore.ieee.org/document/718861/

[3] A. Egido and W. H. F. Smith, "Fully Focused SAR Altimetry: Theory and Applications," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 55, no. 1, pp. 392-406, Jan. 2017. [Online]. Available: http://ieeexplore.ieee.org/document/7579570/

[4] C. J. Donlon, R. Cullen, L. Giulicchi, P. Vuilleumier, C. R. Francis, M. Kuschnerus, W. Simpson, A. Bouridah, M. Caleno, R. Bertoni, J. Rancaño, E. Pourier, A. Hyslop, J. Mulcahy, R. Knockaert, C. Hunter, A. Webb, M. Fornari, P. Vaze, S. Brown, J. Willis, S. Desai, J.-D. Desjonqueres, R. Scharroo, C. Martin-Puig, E. Leuliette, A. Egido, W. H. Smith, P. Bonnefond, S. Le Gac, N. Picot, and G. Tavernier, "The Copernicus Sentinel-6 mission: Enhanced continuity of satellite sea level measurements from space," *Remote Sensing of Environment*, vol. 258, p. 112395, Jun. 2021. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0034425721001139

[5] F. Reale, E. Pugliese Carratelli, A. Di Leo, and F. Dentale, "Wave Orbital Velocity Effects on Radar Doppler Altimeter for Sea Monitoring," *Journal of Marine Science and Engineering*, vol. 8, no. 6, p. 447, Jun. 2020. [Online]. Available: https://www.mdpi.com/2077-1312/8/6/447

[6] P. Rieu, T. Moreau, E. Cadier, M. Raynal, S. Clerc, C. Donlon, F. Borde, F. Boy, and C. Maraldi, "Exploiting the Sentinel-3 tandem phase dataset and azimuth oversampling to better characterize the sensitivity of SAR altimeter sea surface height to long ocean waves," *Advances in Space Research*, vol. 67, no. 1, pp. 253-265, Jan. 2021. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0273117720306840

[7] T. Moreau, E. Cadier, F. Boy, J. Aublanc, P. Rieu, M. Raynal, S. Labroue, P. Thibaut, G. Dibarboure, N. Picot, L. Phalippou, F. Demeestere, F. Borde, and C. Mavrocordatos, "High-performance altimeter Doppler processing for measuring sea level height under varying sea state conditions," *Advances in Space Research*, vol. 67, no. 6, pp. 1870-1886, Mar. 2021. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S027311772030911X

[8] G. Engen and H. Johnsen, "SAR-ocean wave inversion using image cross spectra," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 33, no. 4, pp. 1047-1056, Jul. 1995. [Online]. Available: http://ieeexplore.ieee.org/document/406690/

[9] T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, "A unified directional spectrum for long and short wind-driven waves," *Journal of Geophysical Research: Oceans*, vol. 102, no. C7, pp. 15781-15796, Jul. 1997. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/97JC00467

[10] R. Valenzuela, "AND OCEANIC WAVES - A REVIEW," 1978.
[11] F. Nouguier, A. Mouche, N. Rascle, B. Chapron, and D. Vandemark, "Analysis of Dual-Frequency Ocean Backscatter Measurements at Ku- and Ka-Bands Using Near-Nadir Incidence GPM Radar Data," *IEEE Geoscience and Remote Sensing Letters*, vol. 13, no. 9, pp. 1310-1314, Sep. 2016. [Online]. Available: http://ieeexplore.ieee.org/document/7508981/

[12] B. Chapron, H. Johnsen, and R. Garello, "Wave and wind retrieval from sar images of the ocean," *Annales Des Télécommunications*, vol. 56, no. 11-12, pp. 682-699, Nov. 2001. [Online]. Available: https://link.springer.com/10.1007/BF02995562

[13] H. Li, A. Mouche, H. Wang, J. E. Stopa, and B. Chapron, "Polarization Dependence of Azimuth Cutoff From Quad-Pol SAR Images," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 57, no. 12, pp. 9878-9887, Dec. 2019. [Online]. Available: https://ieeexplore.ieee.org/document/8798973/

[14] F. C. Jackson, "THE PHYSICAL BASIS FOR ESTIMATING WAVE-ENERGY SPECTRA WITH THE RADAR OCEAN-WAVE SPECTROMETER," 1987.

[15] H. E. Krogstad, "A simple derivation of Hasselmann's nonlinear ocean-synthetic aperture radar transform," *Journal of Geophysical Research: Oceans*, vol. 97, no. C2, pp. 2421-2425, Feb. 1992. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/91JC03010

[16] K. Hasselmann and S. Hasselmann, "On the nonlinear mapping of an ocean wave spectrum into a synthetic aperture radar image spectrum and its inversion," *Journal of Geophysical Research: Oceans*, vol. 96, no. C6, pp. 10713-10729, Jun. 1991. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/91JC00302

[17] H. E. Krogstad, O. Samset, and P. W. Vachon, "Generalizations of the non-linear ocean-SAR transform and a simplified SAR inversion algorithm," *Atmosphere-Ocean*, vol. 32, no. 1, pp. 61-82, Mar. 1994. [Online]. Available: http://www.tandfonline.com/doi/abs/10.1080/07055900.1994.9649490

[18] D. Lyzenga, "Numerical Simulation of Synthetic Aperture Radar Image Spectra for Ocean Waves," *IEEE Transactions on Geoscience and Remote Sensing*, vol. GE-24, no. 6, pp. 863-872, Nov. 1986. [Online]. Available: http://ieeexplore.ieee.org/document/4072557/

[19] V. Kerbaol, B. Chapron, and P. W. Vachon, "Analysis of ERS-1/2 synthetic aperture radar wave mode imagettes," *Journal of Geophysical Research: Oceans*, vol. 103, no. C4, pp. 7833-7846, Apr. 1998. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/97JC01579

[20] C. Buchhaupt, A. Egido, W. H. Smith, and L. Fenoglio, "Conditional sea surface statistics and their impact on geophysical sea surface parameters retrieved from SAR altimetry signals," *Advances in Space Research*, vol. 71, no. 5, pp. 2332-2347, Mar. 2023. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0273117722011425

[21] F. Ardhuin, J. Stopa, B. Chapron, F. Collard, M. Smith, J. Thomson, M. Doble, B. Blomquist, O. Persson, C. O. Collins, and P. Wadhams, "Measuring ocean waves in sea ice using SAR imagery: A quasi-deterministic approach evaluated with Sentinel-1 and in situ data," *Remote Sensing of Environment*, vol. 189, pp. 211-222, Feb. 2017. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0034425716304710

[22] T. Yoshida, "Numerical research on clear imaging of azimuth-traveling ocean waves in SAR images," *Radio Science*, vol. 51, no. 7, pp. 989-998, Jul. 2016. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016RS006000

[23] J. Schulz-Stellenfleth, "Spaceborne synthetic aperture radar observations of ocean waves traveling into sea ice," *Journal of Geophysical Research*, vol. 107, no. C8, 2002. [Online]. Available: http://doi.wiley.com/10.1029/2001JC000837

[24] A. Goldfinger, "Estimation of Spectra from Speckled Images," *IEEE Transactions on Aerospace and Electronic Systems*, vol. AES-18, no. 5, pp. 675-681, Sep. 1982. [Online]. Available: http://ieeexplore.ieee.org/document/4102711/

[25] V. Kudryavtsev, "On radar imaging of current features: 1. Model and comparison with observations," *Journal of Geophysical Research*, vol. 110, no. C7, p. C07016, 2005. [Online]. Available: http://doi.wiley.com/10.1029/2004JC002505

[26] M. Kleinherenbrink, M. Naeije, C. Slobbe, A. Egido, and W. Smith, "The performance of CryoSat-2 fully-focussed SAR for inland water-level estimation," *Remote Sensing of Environment*, vol. 237, p. 111589, Feb. 2020. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0034425719306091

[27] F. Ehlers, F. Schlembach, M. Kleinherenbrink, and C. Slobbe, "Validity assessment of SAMOSA retracking for fully-focused SAR altimeter waveforms," *Advances in Space Research*, vol. 71, no. 3, pp. 1377-1396, Feb. 2023. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0273117722010651

[28] F. Ehlers, C. Slobbe, M. Verlaan, and M. Kleinherenbrink, "A noise autocovariance model for SAR altimeter measurements with implications for optimal sampling," *Advances in Space Research*, vol. 71, no. 10, pp. 3951-3967, May 2023. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0273117723001734

[29] R. Husson and P. Vincent, "Sentinel-1 Ocean Swell Wave Spectra (OSW) Algorithm Definition," 2021.

[30] F. Nouguier, C. Guérin, and B. Chapron, ""Choppy wave" model for nonlinear gravity waves," *Journal of Geophysical Research: Oceans*, vol. 114, no. C9, p. 2008JC004984, Sep. 2009. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2008JC004984

Marcel Kleinherenbrink obtained his master's degree at the faculty of Aerospace Engineering of Delft University of Technology, Delft, The Netherlands in 2013. He obtained his Ph.D. degree, on the observation of sea-level change and vertical land motion from space, at the faculty of Civil Engineering and Geoscience of Delft University of Technology in 2018. He specialises in the observation of the ocean surface using active radar instruments. In 2019 he became involved in the Earth Explorer 10 mission Harmony, for which he implements forward models, investigates inversion strategies and provides performance estimates. Since 2018, he has been involved in numerous (ESA) projects related to future missions and novel retrieval algorithms for ocean observations.

Frithjof Ehlers received his Master of Science in Physics of Complex Systems from the University of Goettingen, Germany, in 2019. After that, he gained experience in processing of Remote Sensing Data while working on the Aeolus mission during a Traineeship at ESTEC, Netherlands. Currently, he is pursuing his Ph.D. degree at the department of Geoscience and Remote Sensing, Faculty of Civil Engineering, at the Delft University of Technology. His work is focused particularly on new processing techniques for SAR altimeter data, regarding oceans, coastal zones and inland waters.

Sergi Hernández received the M.S. degree in Telecommunications Engineering from the Universitat Politècnica de Catalunya, Barcelona, in 2020. He is currently pursuing his Ph.D. degree in signal theory and communications at the Universitat Politècnica de Catalunya, Barcelona, co-directed by isardSAT-Cat, Barcelona. In 2019, he joined isardSAT as an Research and Development Engineer. His main activities are the research, development, and implementation of altimetry data processing algorithms.

Frédéric Nouguier received the Agrégation and M.S. degrees in applied physics from the Ecole Normale Supérieure de Cachan, Cachan, France, the M.S. degree in physical methods for remote sensing from the University of Paris-Diderot, Paris, France, and the Ph.D. degree in physics from Aix-Marseille University, Marseille, France, in 2009. He is currently a Research Scientist with the Laboratoire d'Océanographie Physique et Spatiale, Institut Français de Recherche pour l'Exploitation de la Mer, Plouzané, France. He has experience in applied mathematics, physical oceanography, and electromagnetic wave theory and its application to ocean remote sensing.

Ourania Altiparmaki obtained her M.Sc. degree in Geoinformatics at the faculty of Engineering of Aristotle University of Thessaloniki, Greece in 2017. Her master's thesis focused on height system unification applications using satellite altimetry and GNSS data. Currently she is pursuing a Ph.D., specializing in developing advanced image processing techniques to monitor ocean waves from Synthetic Aperture Radar (SAR) altimetry. Her research focuses on satellite missions such as CryoSat-2, Sentinel-3A/B, and Sentinel-6A.
**Florian Schlembach** works as a guest researcher at the Deutsches Geodätisches Forschungsinstitut at the Technical University of Munich, Germany, pursuing a PhD, and is employed as a development engineer at Rohde&Schwarz in signal and spectrum analysis. His research interest is in satellite altimetry and digital signal processing. His academic background lies in the field of electrical engineering. He did his Bachelor's degree in media technology and graduated with his Master degree in 'Communication and Signal Processing' at the Technical University of Ilmenau, Germany. After that, he worked in the R&D industry at Rohde&Schwarz in Munich, Germany, where he worked for five years as a development engineer in mobile communication (LTE 4G/5G).

**Bertrand Chapron** was born in Paris, France, in 1962. He received the B.Eng. degree from the Institut National Polytechnique de Grenoble, Grenoble, France, in 1984, and the Doctorate National (Ph.D.) degree in fluid mechanics from the University of Aix-Marseille II, Marseille, France, in 1988. He was a Post-Doctoral Research Associate with the NASA/GSFC/Wallops Flight Facility, Wallops Island, VA, USA, for three years. He has experience in applied mathematics, physical oceanography, electromagnetic waves theory, and its application to ocean remote sensing. He is currently responsible for oceanography at the Space Laboratory, IFREMER, Plouzané, France.