# On the Effects of Ocean Surface Motion on Delay-Doppler Altimetry

Louis Marié, Frederic Nouguier, Douglas C. Vandemark, Fabrice Ardhuin, Bertrand Chapron

► To cite this version:

Louis Marié, Frederic Nouguier, Douglas C. Vandemark, Fabrice Ardhuin, Bertrand Chapron. On the Effects of Ocean Surface Motion on Delay-Doppler Altimetry. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62, 4207725 (25p.). 10.1109/TGRS.2024.3387667. hal-04731638

HAL Id: hal-04731638

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

Submitted on 16 Jan 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
# On the Effects of Ocean Surface Motion on Delay-Doppler Altimetry

Marié Louis <fnref n="1" />*, Nouguier Frederic <fnref n="1" />, Vandemark Douglas C. <fnref n="2" />, Ardhuin Fabrice <fnref n="3" />, Chapron Bertrand <fnref n="1" />

## Abstract :

The Poseidon-4 radar altimeter on board Sentinel-6 "Michael Freilich" (S6-MF) offers unique opportunities to assess the impact of ocean surface motion on Delay-Doppler altimetry. In this paper, earlier "frozen-sea" studies of the instrument response to an isolated sea surface facet are extended to include the effect of surface motions in its Delay-Doppler Map signature. Integrating this elementary signature over the instrument field of view, an analytical stacked echo waveform model, the IASCO waveform model, is then derived. This waveform is validated against the well-established SAMOSA waveform model for the special case of a frozen sea. Model sensitivity to changes in surface significant wave height, vertical velocity standard deviation, and the "Geophysical Doppler" vector U GD projection along the satellite ground-track velocity are discussed. These developments provide theoretical and analytical means to jointly exploit the S6-MF conventional and Delay-Doppler radar waveforms to improve estimates of Essential Climate Variables (Sea Level, Significant Wave Height), and to retrieve and map two new observables, the along-track projection of the "Geophysical Doppler" vector and the ocean waves vertical velocity standard deviation. These new variables, being sensitive to higher-order spectral moments of the wave directional spectrum, may help to mitigate sea state range bias impacts on altimeter sea level measurements.

**Keywords** : Sentinel-6 MF, Delay-Doppler Altimetry, tracking waveform, surface displacement, surface motion

<footnote n="1">¹ Laboratoire d'Océanographie Physique et Spatiale, Université de Brest, CNRS, Institut Francais de Recherche pour l'Exploitation de la Mer, IRD, Plouzané, France</footnote>

<footnote n="2">² Ocean Process Analysis Laboratory, University of New Hampshire, Durham, NH, USA</footnote>

<footnote n="3">³ Laboratoire d'Océanographie Physique et Spatiale, Université de Brest, CNRS, Institut Francais de Recherche pour l'Exploitation de la Mer, IRD, Plouzané, France</footnote>

<footnote>* Corresponding author : Louis Marié, email address : Louis.Marie@lfremer.fr</footnote>
# I. INTRODUCTION

Satellite radar altimetry provides unique data on sea level, surface wave significant height and surface wind speed. In practice, these observations are obtained by fitting a theoretical function to the satellite-measured waveform, where the waveform is the average of many observed radar echoes as a function of time delay between the transmission and reception of radar pulses. This fitting process is known as retracking, and differing methods involve the choice of both a theoretical waveform shape and a fitting cost function.

Improving the quality of the retrieved altimeter observations thus involves improving not only instrument characteristics, but also the appropriateness of the parametric waveform model used in the retracking process.

Prior to the SIRAL instrument on Cryosat-2, the information contained in the phase of the return signal from successive radar pulses was discarded: only the echo signal intensity was considered, and multi-pulse burst averaging was used to improve the Signal-to-Noise Ratio (SNR) and to mitigate speckle

noise before retracking [1]. Echo waveforms produced in this way were satisfactorily understood in the framework set by the original papers of [2] and [3]. In this “Conventional Altimetry” (CA) framework, the overall waveform is considered to be the sum, on an intensity basis, of the echos returned by the individual scattering elements of the sea surface present in each range bin. It is obtained as the “triple convolution” [3], [4] of a function characterizing the intensity weighting at the range bin level caused by the instrument range resolution mechanism (the range “Point Target Response” (PTR) of the instrument), the probability distribution function of the vertical displacement due to surface waves of scattering elements with respect to their rest positions, and a function describing the radiated power distribution on the sea surface by the system antenna, the “Flat Surface Impulse Response” (FSIR). This does not account for the fact that scattering elements might be displaced from their Doppler frequency bins by surface wave orbital velocities, and assumes that the distribution of surface elevations is independent of the range bin [5].

Retrackers based on this waveform model typically retrieve three geophysical parameters from the radar signal [6]: the normalized radar cross-section $σ_0$, which is related to the surface slope statistics [7]–[9], the standard deviation of the vertical surface displacement, from which the significant wave height $H_s$ is derived, and the mean range from satellite to sea surface over the instrument footprint. The altimeter range is further processed to provide along-track Sea Level and Sea Level Anomaly estimates with cm-level precision.

The SIRAL instrument onboard CryoSat-2 was the first altimeter instrument to allow coherent (phase-preserving) processing of successive radar echos. This permits implementation of the Delay-Doppler Altimetry (DDA) concept proposed by [10], involving Synthetic Aperture Radar (SAR) processing applied to successive echos. This approach allows separation of the signal contributions according to their respective along-track surface-return bins. This yields a two-dimensional Delay-Doppler Map (DDM) of the received power distribution across time delay and Doppler shift. A “Range cell Migration Correction” (RCMC or RMC) is typically applied to waveforms from each Doppler bin to realign them. This removes the along-track distance contribution from the instrument-to-surface slant range. At this stage, signal accumulation (“stacking” or “multi-looking”) is performed before retracking of the resulting 1-D waveform is performed.

Several approaches have been proposed for the signal accumulation stage:

- In the “SAR altimetry” approach [11], SAR-processed radar waveforms viewing the same groundpoint during successive
bursts of pulses are gathered together before summation. Fully-Focussed SAR, in which the final summation itself is phase-preserving, has also been successfully implemented [12]. This can provide observations of the surface elevation of very compact water bodies including lakes, canals [13], and sea-ice leads [14].

— In the alternative “LR-RMC” approach recently proposed by [15], SAR-processing and RMC are applied on bursts of pulses, and this is followed by an incoherent stacking of the intensity. This approach does not yield a finer along-track resolution of the observations, but it does increase the SNR prior to retracking.

The effect of the RMC step, already clearly explained in [10], is to redistribute the total energy contained in the echo waveform, shifting the weight towards the instant of encounter of the radar pulse with the sea surface at nadir. The resulting echo waveform is strongly affected, featuring a sharp peak, instead of a step, at the instant of encounter. This redistribution of energy increases the SNR of the useful portion of the signal, thereby increasing the accuracy of the tracking process.

Despite this potential for improved precision measurement, adoption of DDA was at first hampered by the misconception that it had fundamentally different pulse repetition frequency (PRF) requirements from CA, and that the two modes of operation, as they were implemented for instance in SIRAL, were mutually exclusive. This precluded inter-calibration of the results, and incurred an unacceptably high risk of putting in jeopardy the continuity of the vitally-important sea level observations timeseries started with the 1992 Topex/Poseidon mission [16].

The feasibility of an “open-burst” mode of operation, which could allow simultaneous processing of the same radar echos using the two methodologies, was finally recognized by [16], [17]. This opened the way for the selection of a Delay-Doppler capable instrument for the Sentinel-6 “Michael Freilich” (S6-MF, [18]) mission.

Following this adoption of DDA, there is now a clear need to optimize a new generation of altimeter waveform retracking models that may ultimately serve to improve satellite mission ground segment retracking.

Starting from [19], a number of DDA waveform retrackers have been derived, with varying degrees of sophistication. The studies by [20] and [21] were the first to benefit from and address actual Cryosat-2 data. The authors of [20] propose a lucid exposition of the problem, analytically deriving a number of steps, before resorting to the numerical evaluation of a number of integrals. Results show convincing comparisons with Cryosat-2 ocean return waveforms. The authors of [21], in the context of geodesy, go further analytically, but overlook the issue posed by the finite (and effectively coarse) resolution of the along-track SAR processing. The proposed waveform peak is too narrow for given significant wave height $H_s$, and would have led to consistent overestimation of $H_s$ if the method had been used operationally. In the framework of the ESA-funded SAMOSA project, [22] derived a semi-analytical model that has since been adopted by a large user community. Independently, [23] proposed another semi-analytical model, placing specific emphasis on the retracking

of waveforms when the sensor data are affected by instrument antenna mispointing.

These contributions essentially follow from the original [3] framework and are based on a triple convolution in the range direction only. They do not take into account the effectively-finite azimuthal resolution of the actual SAR processing nor the fact that sea surface movements (and not just height displacements) affect the signature of scattering facets in the Delay-Doppler Map. These factors introduce a PTR and a probability distribution function in the Doppler direction as well. The analysis framework needs to involve a triple convolution in the azimuth direction also.

Several steps in this direction were made by [24]–[27]. These studies propose a computational approach involving a convolution with an azimuth PTR [24], and inclusion of a probability distribution function for the instantaneous radial velocity of the surface scattering elements [25]. They discuss statistics of these elements in relation with the overall sea state characteristics in [26]. Finally, all these elements are consolidated in [27].

The initial motivation for the present study comes from a recent conference presentation [28]. But results to follow have been derived independently from that work and the previously cited studies. Developments are performed within the framework of the ESA-funded “Interpretation and Analysis of SKIM Campaign Observations” (IASCO) project, focused on the analysis of airborne Doppler radar observations of the sea surface collected as part of the phase-A work for the Surface Kinematics Multiscale (SKIM) concept that was a candidate satellite mission for ESA Earth Explorer 9 [29], [30].

The paper is organized as follows. The measurement problem context, including key instrument characteristics and acquisition geometry, is summarized in section II. A derivation of the instrument response to a single scattering facet, the basic building block of the overall instrument response, is discussed in section III. Detailed consideration is given to the problem of statistical averaging over realizations of the random sea surface process, in terms of one facet’s instantaneous elevation and velocity. The idea here is to reproduce for a moving sea surface the initial “frozen-sea” analysis of [3]. A probabilistic theory of these fluctuations is detailed in the Gaussian sea state framework in Appendix A. The integration over the instrument Field of View (FoV), required to upscale the results obtained at the facet level to expected full surface ocean altimeter waveforms, is discussed in section IV. This section culminates with the analytical IASCO SAR altimeter waveform model Eq. (20), which is compared with the Open-Source PySAMOSA [31] implementation of the accepted SAMOSA model [22]. Implications and example results using this model are then discussed in section V followed by a study summary.

## II. PROBLEM CONTEXT

### A. Main instrument characteristics

The S6-MF mission is thoroughly described in [18], but several key characteristics are revisited for the sake of completeness. Relevant S6-MF orbit and Poseidon 4 instrument characteristics are summarized in table I.
The primary mission payload is Poseidon 4, a dual-frequency Ku-/C-band altimeter. A notable difference with respect to the preceding Poseidon systems is that the chirped pulse deramping process is performed numerically after digitization instead of analogically in the receiver front-end. It is also capable to work in the “open-burst” mode, where it continuously alternates between radar transmit and receive operations. This is a significant improvement compared to Cryosat-2 and Sentinel-3A instruments. Those sensors had to alternate between much longer burst transmission and burst reception periods (“closed-burst” mode, see e.g. [32] for an in-depth discussion). This means that the Poseidon 4 signal can be processed using both the CA and DDA approaches at all times. One downside of the open burst mode is that the pulse travel-time sets strong constraints on the maximum transmit pulse repetition frequency $f_p$, which has to be roughly halved with respect to the previous instruments (the implications of this design trade-off are discussed in [17]), and to vary continuously, to adjust for travel-time variations along the platform orbit. Finally, the instrument features a significantly improved onboard processor [18], which is capable of performing onboard the range-resolution, along-track Fourier transformation, and RMC steps of DDA processing. This capability is however not in permanent use, and the DDA processing is currently performed on ground.

<table><caption>TABLE I
SENTINEL 6-MF / POSEIDON-4 DESIGN PARAMETERS.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Mean flight altitude</td><td><i>h</i></td><td>1347 km</td></tr><tr><td>Mean flight velocity</td><td><i>v</i><sub>s</sub></td><td>6967 m.s<sup>-1</sup></td></tr><tr><td>Mean Earth Radius</td><td><i>R</i><sub>E</sub></td><td>6371 km</td></tr><tr><td>Orbital factor</td><td>κ = 1 + <i>h</i>/<i>R</i><sub>E</sub></td><td>1.21</td></tr><tr><td>Average pulse repetition frequency</td><td><i>f</i><sub>p</sub></td><td>9178 Hz</td></tr><tr><td>Carrier frequency (Ku-band)</td><td><i>f</i><sub>c</sub></td><td>13.575 GHz</td></tr><tr><td>Carrier wavelength (Ku-band)</td><td>λ</td><td>2.21 cm</td></tr><tr><td>Pulse duration</td><td><i>T</i><sub>p</sub></td><td>32 µs</td></tr><tr><td>Chirp bandwidth</td><td><i>B</i></td><td>-320 MHz</td></tr><tr><td>ADC sampling frequency</td><td><i>f</i><sub>s</sub></td><td>395 MHz</td></tr><tr><td>I/Q samples per pulse</td><td><i>N</i><sub>s</sub></td><td>20480</td></tr><tr><td>Chirp rate</td><td><i>s</i></td><td>-1.0 10<sup>13</sup> Hz.s<sup>-1</sup></td></tr><tr><td>Range-Doppler ambiguity effective time shift</td><td>δ<i>t</i><sub>rr</sub> = <sup><i>h</i><sub>tr,k</sub></sup>/<sub><i>c</i></sub> + <sup><i>f</i><sub>c</sub></sup>/<sub><i>s</i></sub></td><td>3.13 ms</td></tr><tr><td>Pulses per burst</td><td><i>N</i><sub>p</sub></td><td>64</td></tr><tr><td>Burst duration</td><td><i>N</i><sub>p</sub>/<i>f</i><sub>p</sub></td><td>6.97 ms</td></tr><tr><td>Antenna 3 dB half-beamwidth</td><td><sup>θ<sub>3dB</sub></sup>/<sub>2</sub></td><td>0.665°</td></tr><tr><td>2-way ambiguity velocity</td><td><sup>λ<i>f</i><sub>p</sub></sup>/<sub>4</sub></td><td>~ ±50.7 m.s<sup>-1</sup></td></tr><tr><td>2-way ambiguity <i>y</i></td><td><sup>λ<i>f</i><sub>p</sub> <i>h</i><sub>0</sub></sup>/<sub>4 <i>v</i><sub>s</sub></sub></td><td>~ ±9.8 km</td></tr><tr><td>1-way range diversity over FoV</td><td><i>r</i> - <i>h</i><sub>0</sub> ≃ κ <sup><i>y</i><sup>2</sup></sup>/<sub>2<i>h</i><sub>0</sub></sub></td><td>~ [0; 43.2 m]</td></tr><tr><td>Elevation angle θ at end of unambiguous zone</td><td><sup>λ<i>f</i><sub>p</sub></sup>/<sub>4<i>v</i><sub>s</sub></sub></td><td>0.417°</td></tr></tbody></table>

One other consequence of the reduced $f_p$ of Poseidon 4 is the fact that the Nyquist frequency $f_p/2$ of the SAR processing is also lower, with one ramification that follows. The 1-way power radiation gain diagram $G$ of the instrument antenna can be conveniently written as a Gaussian function depending on the elevation angle $\theta$ [3], [20]:

$$G(\theta) = G_0 \exp\left(-\frac{2\theta^2}{\gamma}\right), \quad \gamma = \frac{\sin^2(\theta_{3dB})}{2\log(2)} \quad (1)$$

where $G_0$ is the boresight gain of the instrument and $\theta_{3dB}$ is the half-power beamwidth. As shown in table I, the elevation angle at the along-track position where the Doppler frequency shift due to the platform flight velocity is equal to the Nyquist frequency is of the order of 0.42°, i.e. actually smaller than the antenna 3 dB half-beamwidth, 0.66°. A significant portion of the instrument Field of View (FoV) is thus aliased in the SAR processing, and appears as conspicuous sidelobes in the DDM.

## B. Measurement geometry

<figure><img src="image_2.png" /><figcaption>Fig. 1. Side-view of the nadir altimeter acquisition geometry. The Earth center, satellite antenna center-of-phase and nadir are respectively denoted $O$, $S$ and $N$. The satellite is at altitude $h$ above the reference spheroid of radius $R_E$. The satellite flight velocity with respect to the Earth surface is $\vec{v}_s$, with norm $v_s$, tangential component $v_t$, radial component $v_r$ (positive upwards). At the current observation point $P$, the Earth surface is at height $\eta$ above the reference spheroid of radius $R_E$, the elevation angle is denoted $\theta$, the satellite zenith angle is denoted $\psi$, the geocentric angular distance from nadir is $\alpha$. The vertical distance between the reference spheroid and its tangent plane at nadir is denoted $\Delta$ (it is always negative). The vector from $S$ to $P$ is denoted $\vec{r}$, and its norm, the radar slant range to $P$, is denoted $r$.</figcaption></figure>

1) *From the side:* Figure 1 presents a side-view of the satellite trajectory and acquisition geometry in a vertical plane that is locally tangent to the trajectory in the Earth-bound frame [33]. A number of useful coordinate systems can be defined (see also Fig. 2). The current observation point $P$ can be referred to by its geocentric angular distance from nadir $\alpha$, azimuth with respect to the satellite ground track $\varphi$ and distance from the Earth center $R_E+\eta$, or by its cylindrical $(\rho, \varphi, z = \Delta+\eta)$ or Cartesian $(x, y, z = \Delta+\eta)$ coordinates in the frame of reference locally tangent to the sphere and centered
on the nadir. $\Delta$ denotes the distance in the z direction between the reference spheroid and the plane tangent to the sphere at nadir. It is always negative. In this work the y-direction of the Cartesian coordinate system is parallel to the satellite velocity vector.

Given the typical beam apertures used in nadir altimetry, $\rho$ is of the order of at most $10^4$ m. At such scales, it is easy to see that the observation point distances from nadir following the spherical surface, $\zeta$, and along the tangent plane ($\rho$, equal to $\sqrt{x^2 + y^2}$), differ by less than 5 mm, and can be used interchangeably. From this one gets that $\theta \simeq \tan(\theta) = \rho/h$.

It is well known [4] that Earth sphericity must be taken into account for its vertical deflection of the surface, with $\Delta$ varying with $\rho$ as $\Delta = -\rho^2/(2R_E)$ and reaching values of the order of 8 m at the edge of the instrument beam. This effect is also found in the expression of the slant range from $S$ to $P$,

$$r = \sqrt{\left(h - \eta + \frac{\rho^2}{2R_E}\right)^2 + \rho^2} \approx h - \eta + \frac{\kappa}{2h}\rho^2,$$

which increases faster with $\rho^2$ than in a flat surface approximation by a factor $\kappa = 1 + h/R_E$, equal to 1.21 for S6-MF flight parameters. And although $\theta$ is well approximated by a flat-Earth expression $\rho/h$, the satellite zenith angle must also be corrected for sphericity and is given by

$$\psi = \kappa \frac{\rho}{h}.$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Top-view of the nadir altimeter acquisition geometry. The nadir is denoted by $N$. The y-direction of the Cartesian coordinate system is parallel to the ground projection of the satellite flight velocity vector, $\bar{v}_s$. The azimuth angle $\varphi$ is counted clockwise from the y-direction. The “Geophysical Doppler” vector is $\bar{U}_{GD}$. The current point $P$ can be referred to by its cylindrical coordinates ($\rho, \varphi, \Delta + \eta$), its Cartesian coordinates ($x, y, \Delta + \eta$), or by its Delay-Doppler coordinates $r = \sqrt{(h - \Delta - \eta)^2 + \rho^2}$ and $f$. The isorange lines are represented by the circular grid centered on $N$, the isodops in absence of mean surface motion by the dashed lines grid, and the actual isodops showing the effect of mean surface motion by the thin lines grid.</figcaption></figure>

2) *From above:* Figure 2 represents the acquisition geometry seen from above. The satellite flight velocity aligns with

the positive y-direction. The current point $P$ on the surface can be labeled using different sets of coordinates, Cartesian ($x, y, z = \Delta + \eta$), cylindrical ($\rho, \varphi, z = \Delta + \eta$), or, more interestingly, Delay-Doppler coordinates that are constructed from its slant range from the satellite, $r$, equal to $h - \eta + \frac{\kappa}{2h}\rho^2$, and the Doppler shift at its location, $f$.

At a given location, the Doppler shift due to a dynamic ocean is composed of a “frozen-sea” contribution and a wave-induced contribution emerging from the correlation of the Lagrangian velocity of the specular facets and their backscattering cross-sections. This contribution is discussed in section III and Appendix A, where it is shown to represent the projection on the radar line-of-sight of an additional “Geophysical Doppler” vector $U_{GD}$. In [30], this vector is further decomposed into a “Current Doppler” vector $U_{CD}$, i.e. the ocean “Total Surface Current Vector” and a “Wave Doppler” vector due to the intrinsic motion of the sea surface specular facets. Note, the “Wave Doppler” term is the dominant component (see Eq. 2 in [34]), with magnitude of the order of 2.5 to 3.5 m.s$^{-1}$ at Ku-band, and a weak dependence on sea state and wind speed.

Denoting the magnitude of this vector, $U_{GD}$, and the azimuth it points to, $\varphi_{GD}$, the expression for the total Doppler frequency shift can be read from Eq. (13) as

$$f = \frac{2}{\lambda} \left[ -v_r + \frac{yv_t}{h} - \tan(\psi) U_{GD} \cos(\varphi - \varphi_{GD}) \right],$$

with $\lambda$ the wavelength of the carrier radio waves. The radial velocity term is discussed in section III. This can be readily removed as a “Doppler Centroid Anomaly”.

Turning to the Geophysical Doppler contribution, $U_{GD}$ is small with respect to $v_t$. Introducing $\varepsilon$, the fractional Doppler shift change observed due to the surface motion when looking along the satellite flight direction, equal to

$$\varepsilon = -\kappa \frac{U_{GD}}{v_t} \cos(\varphi_{GD}),$$

and the deflection angle

$$\varphi_{\varepsilon} = -\kappa \frac{U_{GD}}{v_t} \sin(\varphi_{GD}),$$

$f$ can be expressed to first order in $U_{GD}/v_t$ as

$$f = -\frac{2v_r}{\lambda} + \frac{2v_t}{\lambda h} \rho(1 + \varepsilon) \cos(\varphi - \varphi_{\varepsilon}). \quad (2)$$

When the satellite travels “against the wind”, $\varphi_{GD} \simeq 180^{\circ}$, and $\varepsilon$ is positive, indicating an excess Doppler shift with respect to the “frozen sea” situation. For the Sentinel-6 MF flight parameters, $\varepsilon$ is of order $4.3 \ 10^{-4}$, and $\varphi_{\varepsilon}$ is of order $0.03^{\circ}$ at most. The frozen-sea isodops are represented as dashed lines in figure 2, while isodops accounting for a (much exaggerated) Geophysical Doppler contribution are represented in thin lines.

### III. NEAR-NADIR DELAY/DOPPLER RADAR: AT THE FACET LEVEL

Near-nadir radar backscattering is well represented in the so-called Physical Optics, Kirchhoff, approximation. The electromagnetic field backscattered to the instrument from the sea
surface is dominated by the weighted composition of individual returns from distributed stationary-phase (i.e. “specular”) surface facets [35], [36]. Generic life history traits of such specular points were addressed in a Gaussian sea surface context by [37]. In short, such “facets” appear and disappear at well-defined instants, and have well-defined trajectories in between. They thus constitute meaningful objects, and can be ascribed Lagrangian trajectories and velocities during their lifetime.

The phase of each contribution to the overall field is given by the two-way facet distance to the phase center of the instrument antenna, modulo the wavelength. The phases are considered random, uniformly distributed and statistically independent. Upon taking products of the field with itself at the same or different instants or locations and performing ensemble averages, only terms involving the product of the contribution of one facet with its complex-conjugate can survive. Practically, this means that interference between the different scattering elements within the scene do not contribute to the average total received power.

These assumptions lead to Rayleigh statistics for the field, and to the fact that the average backscattered power is proportional to the illuminated surface even though the computed field expectation value remains zero, in analogy with the phenomenology of the Wiener process. Two other consequences are that the instrument received power PTR is equal to the squared modulus of its amplitude PTR, and that applying the RMC step of DDA to complex amplitudes or detected power yields the same result for the average power Delay/Doppler Map (DDM).

What follows next is a derivation dedicated to modeling the Delay/Doppler instrument power response to the return from a single facet. This is then used as a basic building block to model and more fully understand near-nadir DDM measurements from a multi-faceted moving sea surface.

### A. Facet kinematics

Consider an isolated facet whose rest position is $(x_0, y_0, z_0 = \Delta + \eta)$. It is affected by surface wave motion, and at the middle of the observation interval it is located at $(x_0 + \tilde{x}, y_0 + \tilde{y}, z_0 + \tilde{z})$ and moves with respect to the Earth with a Lagrangian velocity $(\tilde{u}, \tilde{v}, \tilde{w})$. Assume this Lagrangian velocity to be constant, and small with respect to both the satellite velocity and the speed of light. If we denote by $\tau_n$ the instant of transmission of the n-th radar pulse with respect to the middle point of the observation interval, and by $t'$ the flight time of the pulse since its transmission, the distance from the radar to the facet evolves with time $\tau_n + t'$ as:

$$r(\tau_n + t') = \sqrt{\begin{array}{l} (x_0 + \tilde{x} + \tilde{u}(\tau_n + t'))^2 \\ + (y_0 + \tilde{y} + (\tilde{v} - v_t)(\tau_n + t'))^2 \\ + (h - (z_0 + \tilde{z}) + (v_r - \tilde{w})(\tau_n + t'))^2 \end{array}}$$

In this expression, the variations of all terms are small with respect to the instrument range resolution over the ranges of variation of both $t'$ and $\tau_n$: all variations can be neglected in the analysis of CA processing. Some terms, however, can

vary by a sizable fraction of a wavelength over these time scales. These have an impact on the phase of the return signal from the facet, and must be taken into account in the analysis of coherent processing algorithms. A robust approximate expression, neglecting terms varying by less than roughly $10^{-2}$ wavelengths over the ranges of variations of $t'$ and $\tau_n$, is

$$r(\tau_n + t') \approx h - (\eta + \tilde{z}) + \frac{\kappa \rho_0^2}{2h} + \left[ v_r - y_0 \frac{v_t}{h} - \tilde{w} \right] \tau_n + \left[ v_r - y_0 \frac{v_t}{h} \right] t' + \frac{v_s^2 \tau_n^2}{2h}.$$

The last term in this expression can be relevant in the Fully-Focused SAR altimetry context, for which coherent processing is used for $\tau_n$ as large as 1 s. On shorter time scales relevant for DDA, it remains negligibly small. The slant range from the radar to the facet is thus obtained as a function of the slow ($\tau$) and fast ($t'$) times as:

$$r(\tau_n + t') \approx \underbrace{h + \frac{\kappa \rho_0^2}{2h} - \eta - \tilde{z}}_{r_0} + \underbrace{\left[ v_r - y_0 \frac{v_t}{h} - \tilde{w} \right] \tau_n}_{r_0 \tau_n} + \underbrace{\left[ v_r - y_0 \frac{v_t}{h} \right] t'}_{r_0 t'} \quad (3)$$

The vertical component of the Lagrangian velocity of the facet manifests itself as a slow pulse-to-pulse phase drift. The other facet velocity components are much too small to be felt. The large and dominating contributions are due to the platform velocity and must be taken into account on both the slow (inter-pulse, $\tau$) and the fast (intra-pulse, $t'$) time scales.

### B. Radar signal generation and reception

Consider a burst of $N_p$ radar pulses, numbered from 0 to $N_p - 1$, transmitted every $1/f_p$ seconds from the satellite between $-T_b/2$ and $T_b/2$. Pulse number $n$ is transmitted at instant $\tau_n$. Time between pulse transmission and reception is measured with the “fast time” variable $t'$. Transmission is from $t' = -T_p/2$ to $t' = T_p/2$. The return signal is recorded during a time window centered around an estimate of the two-way pulse travel time. This is provided by an onboard tracker as $t'_{trk} = 2h_{trk}/c$. This time delay can be considered constant over the burst duration and the recording window is long enough to contain the entire useful length of the return signal<fnref n="1" />.

Between $t' = -T_p/2$ and $t' = T_p/2$, the instantaneous transmitted frequency changes linearly from $f_c - B/2$ to $f_c + B/2$. Introducing the chirp rate $s = B/T_p$, and denoting by $\mathbb{1}_{[a;b]}(t)$ the function equal to 1 for $a \le t \le b$ and 0 otherwise, the transmitted signal thus reads:

$$S_{TX}(t') = \mathbb{1}_{[-T_p/2;T_p/2]}(t')e^{j2\pi f_c t'}e^{j2\pi \frac{s t'^2}{2}}.$$

<footnote n="1">¹A number of delicate intermediate steps are required to make this technologically feasible, in order to accommodate digitizer bandwidth or memory constraints, such as introduction, and subsequent compensation of, a digitization window delay with respect to the pulse start, pulse interleaving control, PRI modulation to accommodate altitude variations, etc. These steps are very well controlled, and we will in the following assume that the artefacts they introduce have been perfectly compensated.</footnote>
At first order in $\dot{r}_0/c$, the signal received at time instant $\tau_n + t'$ left the instrument at time $\tau_n + t'(1 - 2\dot{r}_0/c) - 2\frac{\tau_n}{c}(1 - \dot{r}_0/c)$. Introducing a slowly varying complex amplitude $a(\tau_n)$ to keep track of the effect on the field of the reflection on the facet, and up to the geometric factors required to take into account the antenna radiation pattern, spherical divergence and propagation losses (which can all be reinstated at the end of the calculation), the received signal is:

$$S_{RX}(\tau_n + t') = a(\tau_n)S_{TX}\left(t'\left(1 - 2\frac{\dot{r}_0}{c}\right) - 2\frac{r_n}{c}\left(1 - \frac{\dot{r}_0}{c}\right)\right). \quad (4)$$

Upon reception, the signal is amplified and analogically multiplied by $e^{-j2\pi f_c t'}$, prior to recording. During the recording window, the return signal is “de-ramped” (for a description of the “full de-ramp” technique see e.g. [38]) by a replica of the transmitted chirp, either with hardware prior to digitization or, as in the Poseidon-4 instrument, numerically after digitization. This replica is centered on $t' = t'_{trk}$, and is better described as a function of yet another time scale, $t = t' - t'_{trk}$. The result is a digitized recording of

$$S(\tau_n, t) = e^{-j2\pi(f_c t' + \frac{st^2}{2})} S_{RX}(\tau_n + t'_{trk} + t).$$

### C. Range compression

Using the expression of $S_{RX}$ given in Eq. (4) and neglecting terms amounting to phase shifts smaller than $10^{-2}$ cycles, the expression of the digitized signal segments is

$$S(\tau_n, t) = a(\tau_n) \mathbb{1}_{[-T_p/2; T_p/2]} \left( t + 2\frac{h_{trk} - r_n}{c} \right) \\ \times e^{j4\pi \left[ s \frac{(r_n - h_{trk})^2}{c^2} - \frac{r_n}{\lambda} - \frac{r_0}{c} \frac{h_{trk}}{\lambda} \right]} \\ \times e^{-j4\pi \frac{st}{c} \left[ r_n + \frac{f_c r_0}{s} - h_{trk} \left(1 - \frac{r_0}{c}\right) \right]}.$$

The term on the first line represents the envelope of the pulse. The term on the second line is a phase term, constant over one pulse, but slowly changing from pulse to pulse. The final term describes the time dependence of the deramped echo on fast time. This time-dependence is used to obtain the response of the instrument to the facet as a function of observation range<fnref n="2" />, through a Fourier transform on $t$. Using the fact that the independent variable $r$ is linked to the analysis frequency $\omega$ by the relationship $r = h_{trk} - \frac{c\omega}{4\pi s}$, one obtains the radar

echo amplitude of the facet as a function of $r$ at time $\tau_n$ as:

$$S(\tau_n, r) = a(\tau_n) e^{j4\pi \left[ s \frac{(r_n - h_{trk})^2}{c^2} - \frac{r_n}{\lambda} - \frac{r_0}{c} \frac{h_{trk}}{\lambda} \right]} \\ \times \frac{1}{T_p} \int \mathbb{1}_{[-T_p/2; T_p/2]} \left( t + 2\frac{h_{trk} - r_n}{c} \right) \\ \times e^{j4\pi \frac{st}{c} \left[ r - r_n - \dot{r}_0 \left( \frac{h_{trk}}{c} + \frac{f_c}{s} \right) \right]} dt \\ = a(\tau_n) e^{j4\pi \left[ s \frac{(r_n - h_{trk})^2}{c^2} - \frac{r_n}{\lambda} - \frac{r_0}{c} \frac{h_{trk}}{\lambda} \right]} \\ \times \Upsilon_r \left( r - r_n - \dot{r}_0 \left( \frac{h_{trk}}{c} + \frac{f_c}{s} \right) \right).$$

In this expression $\Upsilon_r(\delta r) = \text{sinc}(\frac{2\pi s T_p}{c} \delta r)$, with $\text{sinc}(x) = \frac{\sin(x)}{x}$, denotes the range-compression amplitude impulse-response function, which is symmetric and peaking around 0. At each instant, the range-resolved waveform reaches its peak close to the actual facet range $r_n$, but includes cm-scale corrections for the motion of the instrument during the travel time of the radar pulse ($\dot{r}_0 h_{trk}/c$ term) and for delay-Doppler ambiguity ($\dot{r}_0 f_c/s$ term).

The first of these two corrections can be absorbed in a change of the definition of $\tau_n$, from the instant of pulse transmission to the instant the pulse strikes the facet (this convention has been used by [22]). The latter correction is of a different nature, and has its origin in the range-resolution method itself: the only differences allowing discrimination between time-shifted and Doppler-shifted versions of a chirp occur at its ends. This is a extremely short part of its total duration. A target observed from a moving instrument thus appears shifted from its actual location, by an amount that depends on the relative speed and on the chirp bandwidth and direction. The facet thus appears shifted by $\dot{r}_0 \delta t_{rr}$ in the range-resolved waveform, with $\delta t_{rr} = \frac{h_{trk}}{c} + \frac{f_c}{s}$ an effective time shift, equal to 3.13 ms for the S6-MF flight parameters. The sensitivity of this term on $h_{trk}$ is small, and it is adequately evaluated using the nominal flight altitude $h$. From Eq. (3), $\dot{r}_0 = v_r - y_0 v_t / h$. The portion of this range correction due to the radial velocity of the platform $v_r$ is uniform in $y$ and identically affects the whole scene. More importantly, the tangential velocity $v_t$ component depends on the along-track location of the facet. This can only be removed after it is estimated in the Doppler-resolution domain.

The waveform peak itself has a sinc shape, which stems from the amplitude envelope of the chirp. The Fourier-space representation of its square, which will be needed in the following, is readily derived to be the “triangle” function

$$\widehat{\Upsilon}_r^2(\mathcal{K}) = \sigma_r \mathbb{1}_{[-1;1]} \left( \frac{\mathcal{K} \sigma_r}{2\pi} \right) \left[ 1 - \frac{|\mathcal{K}| \sigma_r}{2\pi} \right], \quad (5)$$

with $\sigma_r = \frac{c}{2B}$.

For the Sentinel 6-MF parameters, the half-power width of this function is equal to $0.886 \times \frac{c}{2B} \approx 0.415$ m. Classically [2], [3], it is approximated as a Gaussian, though some confusion exists regarding the appropriate width. One choice is to select the value for which the -3 dB widths of the sinc shape and the

<footnote n="2">²In the original article of [10], the range compression is performed after the along-track Fourier transform. The order of these two stages is unconsequential, and we have based our description on the order used in the S6-MF onboard processor.</footnote>
Gaussian approximation coincide. For the power point target response (PTR) in the range direction, this yields

$$ \Upsilon_r^2(\delta r) \approx e^{-\frac{\delta r^2}{2\sigma_r^2}}, $$

using this time $\sigma_r = \frac{1}{\sqrt{2\log(2)}} \frac{0.886c}{4B} = 0.176$ m. The equivalent expression in Fourier space reads

$$ \widehat{\Upsilon}_r^2(\mathcal{K}) \approx \sqrt{2\pi}\sigma_r e^{-\frac{\kappa^2\sigma_r^2}{2}}. \quad (6) $$

### D. Doppler resolution

The next step in the evaluation of the contribution of the individual facet to the DDM is to study the effect of the Doppler frequency resolution. This step is performed through a weighted Fourier transform in slow time. The discrete sampling in $\tau$ tends to complicate notations. It can obscure those characteristics of the DDM that are due to geophysical processes with others that are due to technicalities (typically, aliasing in Doppler frequency). We will thus conduct the analysis as if the sampling was continuous in slow time, and discuss the changes introduced by discrete sampling when relevant. In this framework, the echo amplitude waveform due to the facet at slow time $\tau$ can be modeled as

$$ S(\tau, r) = a(\tau) e^{j4\pi \left[ s \frac{(r_0 - h_{trk})^2 + \tau^2 r_0^2}{c^2} - \frac{r_0}{\lambda} - \frac{h_{trk}}{\lambda} \frac{r_0}{c} \right]} \\ \times e^{-j4\pi \frac{\tau r_0}{c}} \left[ f_c - 2s \frac{r_0 - h_{trk}}{c} \right] \\ \times \Upsilon_r \left( r - [r_0 + \tau r_0' + \dot{r}_0 \delta t_{rr}] \right). $$

The phase term on the first line of this expression is essentially constant over the duration of one burst, and can in fact be subsumed in the complex return signal amplitude $a(\tau)$. The term on the last line encodes the fact that the return waveform is essentially zero, except in the close neighborhood of the actual facet range. Of the phase terms in the second line, it can be seen that only the carrier frequency term can induce a noticeable contribution. It is the Doppler shift term the DDA method is based on. The others can be neglected.

The amplitude signature of the facet in the Delay-Doppler plane<fnref n="3" /> is obtained by performing a weighted Fourier transform in $\tau$:

$$ S(f, r) = \frac{1}{T_b} \int w(\tau) a(\tau) e^{-j2\pi\tau \left(f + 2\frac{r_0}{\lambda}\right)} \\ \times \Upsilon_r \left( r - \left[r_0 + \tau r_0' + \dot{r}_0 \delta t_{rr}\right] \right) d\tau. $$

In this expression, $w(\tau)$ is a “window” weighting function, used to control the shape of the impulse response function along the $f$ dimension.

A small $\tau$-varying term appears in the argument of the range-compression impulse response function. For a fixed

value of $r$, this term amounts to a slowly-varying perturbation of the window function used in the Fourier transform. This gives more weight to the beginning, center, or end of the measurement interval depending on whether the facet is entering, dwelling in, or leaving the neighbourhood of $r$. The effect of this term is small, broadening slightly the result in the range direction for the large $r_0'$ bins. Neglecting it (which amounts to using for $\Upsilon_r$ its value at the middle of the observation interval), we obtain

$$ S(f, r) = \Upsilon_r \left( r - [r_0 + \dot{r}_0 \delta t_{rr}] \right) \\ \times \frac{1}{T_b} \int w(\tau) a(\tau) e^{-j2\pi\tau \left(f + 2\frac{r_0'}{\lambda}\right)} d\tau. $$

### E. Ensemble averaging over the sea surface realizations

Taking the squared modulus of $S(f, r)$, the power response for this realization of the facet is:

$$ |S|^2(f, r) = \Upsilon_r^2 \left( r - [r_0 + \dot{r}_0 \delta t_{rr}] \right) \\ \frac{1}{T_b^2} \iint w(\tau) w(\tau') a(\tau) a^*(\tau') e^{j2\pi(\tau'-\tau)} \left( f + 2\frac{r_0'}{\lambda} \right) d\tau d\tau'. $$

The mathematical expectation of the DDM contribution per unit projected area of the facets located around $(x_0, y_0, \Delta+\eta)$ must be obtained by averaging this individual contribution over the possible realizations of the sea surface.

Sources of randomness that must be considered are $N$, the actual number of facets, with probability $p_N$, dependent on the location and area of the patch, and for each facet, labelled $n$,

* its corresponding vertical displacement due to waves, $\tilde{z}_n$,

* its corresponding vertical Lagrangian velocity, $\tilde{w}_n$,

* its corresponding complex reflecting amplitude at times $\tau$ and $\tau'$, $a_n(\tau)$ and $a_n(\tau')$.

Hence, in the most general setting, the overall DDM contribution of a unit area patch located in the neighbourhood of $(x_0, y_0, \Delta+\eta)$ should be expressed as:

$$ \langle |S| \rangle^2 (f, r; x_0, y_0, \Delta + \eta) = \frac{1}{T_b^2} \sum_{N=0}^{\infty} p_N(x_0, y_0, \Delta + \eta) \sum_{n=0}^{N} \int \dots \int P_{sp}(\tilde{z}_n, \tilde{w}_n, a_n(\tau), a_n^*(\tau'); x_0, y_0, \Delta + \eta) \\ \Upsilon_r^2 \left( r + \tilde{z}_n - h + \eta - \frac{\kappa}{2h} \rho_0^2 - \dot{r}_0 \delta t_{rr} \right) \\ w(\tau) w(\tau') a_n(\tau) a_n^*(\tau') e^{j2\pi(\tau'-\tau)} \left( f - \frac{2}{\lambda} [\tilde{w}_n + y_0 \frac{v_t}{h} - v_r] \right) d\tilde{z}_n d\tilde{w}_n d a_n(\tau) d a_n^*(\tau') d\tau d\tau'. $$

The probability distribution function $\sum_N p_N \sum_n^N P_{sp}(\dots)$ is a very high-dimensional function, of which only specific dependencies are currently known:

* The number of facets $N$ and their complex reflection coefficient $a$ are not known separately, but [36] $\sum_N p_N \sum_n^N \int |a|^2_n P_{sp}(|a_n|^2) d|a_n|^2$ is proportional to the normalized sea surface backscatter cross section, $\sigma_0$, whose statistics and dependencies with respect to

<footnote n="3">³This is an instance where the discrete sampling in $\tau$ introduces a difference with the time-continuous analysis, by making the impulse-response periodic in $f$ with ambiguity period $f_p$. While this ambiguity is unconsequential for the high-$f_p$ instruments of CryoSat-2 and Sentinel-3, for which the sidelobes correspond to very large along-track distances and are rejected by the antenna diagram, the factor of two reduction in $f_p$ of the Poseidon-4 instrument lets very conspicuous sidelobes enter the DDM. The power content in those sidelobes significantly contributes to the stacked DDA power waveform. They can be accounted for analytically, as will be shown in section IV-D.</footnote>
environmental conditions and observation geometry have been extensively studied.

* The statistics of $\tilde{z}_n$, independently of the other variables, have also been studied extensively. A Gaussian prescription is a reasonable assumption, though higher-order contributions likely exist, e.g. altimetric Sea State Bias (SSB).

* Statistical correlations of $\tilde{z}_n$ and $\sigma_0$ have been studied, and are known to result in altimetric electromagnetic bias (EMB).

* The joint statistics of $\tilde{z}_n$ and $\tilde{w}_n$ have been studied theoretically [39]. In the Gaussian context<fnref n="4" />, these variables are independent.

* On the contrary, even in the Gaussian context, $\tilde{w}_n$ is not independent of the local sea surface slope $\nabla \tilde{z}$, which in turn is correlated with the backscatter cross-section. A correlation between surface backscatter and along-line-of-sight surface velocity projection is thus expected, giving rise to the “Geophysical Doppler” bias.

* The slow-time dynamics of the complex backscatter amplitude $a_n$ are difficult to examine experimentally, and thus remain elusive, despite their fundamental importance in determining the optimal time window for coherent-radar observations of the sea surface. Choosing a long time window for the coherent processing is not beneficial if the lifetime of the individual facets is short, as the expected resolution improvement is not achieved, but the opportunity to obtain independent looks is wasted. For instance, the effective Doppler (hence along-track) point target response function can be computed under the assumptions of a Gaussian dependence on $\tau - \tau'$ of the correlation function $\langle a_n(\tau)a_n^*(\tau') \rangle$ and of a Gaussian window function. This exercise shows that it is the shortest time scale (facet coherence time or window function width) that sets the effective Doppler resolution.

Given all these knowledge gaps, we pursue the analysis according to the following practical assumptions:

* We consider all facets to be independent. The sums over facet number $N$ and index $n$ reduce to multiplication by $\langle N \rangle$, the mathematical expectation of $N$. $\langle N \rangle$ depends on the observation geometry.

* We consider the elevations of the facets $\tilde{z}$ to be Gaussian and independent of all other variables (i.e. the study of SSB and EMB is left for future investigation).

* We consider the correlation time of the complex reflecting amplitude $a_n$ to be long enough with respect to the observation window duration, i.e. its dynamics can be neglected. This assumption is probably at least marginally correct over the duration of one burst (6.97 ms). With this assumption, the product of Fourier transforms over slow time simplifies into a Doppler resolution PTR,

$\Upsilon_f^2(\delta f) = \frac{1}{T_b^2} \left| \int_{-T_b/2}^{T_b/2} w(\tau) e^{-j2\pi\tau\delta f} \right|^2$. As discussed in [22], if the so-called “Hamming” window is used for $w(\tau)$, this PTR can be extremely well approximated by the Gaussian shape

$$ \Upsilon_f^2(\delta f) \approx \left(\frac{25}{46}\right)^2 e^{-\frac{\delta f^2}{2\sigma_f^2}}, $$

$$ \text{with } \sigma_f = \frac{1}{\sqrt{2 \log(2)}} \frac{1.293}{2T_b} \approx 78.74 \text{ Hz.} $$

* After [35], [36], the cross section of an isolated specular point is expressed as $|a_n|^2 = |R|^2\pi/|\tilde{\Omega}_n|$, with $|R|^2$ the normal-incidence Fresnel power reflection coefficient and $\tilde{\Omega}_n$ the Gaussian curvature at the specular point.

## F. Averaging over $\tilde{z}$: recovering the “triple-convolution model”

With these assumptions, the expression of the DDM contribution per unit area of the neighbourhood of $(x_0, y_0, \Delta + \eta)$ reads

$$ \begin{aligned} \langle |S|^2 (f, r; x_0, y_0, \Delta + \eta) \rangle &= \pi |R|^2 \langle N \rangle \\ & \times \int P_{sp}(\tilde{z}) \Upsilon_r^2 (r + \tilde{z} - h + \eta - \frac{\kappa}{2h} \rho_0^2 - r_0 \delta t_{rr}) \, d\tilde{z} \\ & \times \iint P_{sp}(\tilde{w}, \tilde{\Omega}; x_0, y_0) \Upsilon_f^2 \left( f - 2 \frac{\tilde{w}}{\lambda} + 2 \frac{r_0}{\lambda} \right) \frac{d\tilde{w} d\tilde{\Omega}}{|\tilde{\Omega}|}, \end{aligned} \quad (8) $$

where $\Upsilon_r^2$ and $\Upsilon_f^2$ are approximately Gaussian, with specified widths, and $P_{sp}(\tilde{z})$ and $P_{sp}(\tilde{w}, \tilde{\Omega})$ respectively provide the (Gaussian) pdf of the elevation and joint pdf of the vertical velocity and the total curvature of the surface at specular facets. Since we have chosen to neglect $\tilde{z}$'s correlations with the other variables, $P_{sp}(\tilde{z})$ is identical to the full-surface pdf of $\tilde{z}$, $P(\tilde{z})$.

The first integral is thus the convolution of a Gaussian with the instrument range PTR:

$$ \tilde{\Upsilon}_r^2(\delta r) = \int P(\tilde{z}) \Upsilon_r^2(\delta r + \tilde{z}) \, d\tilde{z}. \quad (9) $$

This integral defines an effective range PTR, accounting for both the effect of surface waves and the instrument PTR (such effective quantities will be denoted by a $\tilde{\cdot}$ symbol in the following). It is the “frozen-sea” contribution to the detected interface thickness. This is accounted for by the “triple convolution model” treatment of CA and DDA.

In the Gaussian sea state approximation,

$$ P(\tilde{z}) = \frac{1}{\sqrt{2\pi}\sigma_h} e^{-\frac{\tilde{z}^2}{2\sigma_h^2}}, $$

where the surface vertical displacement standard deviation $\sigma_h$ is related to the significant wave height $H_s$ by

$$ \sigma_h = \frac{H_s}{4}. $$

Using the Gaussian approximation of the instrument PTR given by Eq. (6), the effective range PTR can thus be readily approximated in physical space as

$$ \tilde{\Upsilon}_r^2(\delta r) \approx \frac{\sigma_r}{\tilde{\sigma}_r} e^{-\frac{\delta r^2}{2\tilde{\sigma}_r^2}}, $$

<footnote n="4">⁴By “Gaussian context” we mean the situation addressed in references [37], [39]–[41], in which the Fourier components of the sea surface elevation field are independently distributed. In this situation, the sea surface elevation and its derivatives with respect to time and space are jointly Gaussian, and their joint pdfs can be expressed in terms of their second-order moments only. An obvious limitation of this approach is that it can only represent statistically homogeneous surfaces. It can for instance not account for the modulation of the small-scale statistics of the surface by its long wavelength components.</footnote>
with

$$ \sigma_r = \frac{1}{\sqrt{2 \log(2)}} \frac{0.886c}{4B}, \quad \tilde{\sigma}_r = \sqrt{\sigma_r^2 + \sigma_h^2}. $$

A representation using the exact range PTR expression Eq. (5) can be obtained in Fourier space as

$$ \tilde{\Upsilon}_r^2(\mathcal{K}) = \sigma_r \mathbf{1}_{[-1;1]} \left( \frac{\mathcal{K}\sigma_r}{2\pi} \right) \left[ 1 - \frac{|\mathcal{K}|\sigma_r}{2\pi} \right] e^{-\frac{\kappa^2 \sigma_h^2}{2}}, $$

where $\sigma_r = \frac{c}{2B}$. Except where explicitly mentioned, all graphs in this article are produced using this formulation.

We now turn to the second integral in Eq. (8), which involves the probability density function $P_{sp}(\tilde{w}, \tilde{\Omega})$.

### G. Averaging over $\tilde{w}$: the effect of wave motion on the Doppler PTR

As shown in [40], [41], a specular-point pdf follows from the full-surface pdf as

$$ P_{sp}(\tilde{w}, \tilde{\Omega}) = \frac{|\tilde{\Omega}| \left[ 1 + \partial_x \tilde{z} \Big|_{sp}^2 + \partial_y \tilde{z} \Big|_{sp}^2 \right]^2}{\langle N \rangle} \times P(\partial_x \tilde{z} = \partial_x \tilde{z} \Big|_{sp}, \partial_y \tilde{z} = \partial_y \tilde{z} \Big|_{sp}, \tilde{w}, \tilde{\Omega}). $$

Under a Gaussian assumption, this equation shows that the curvature dependent terms in Eq. (8) in fact compensate: facets with small curvature have large cross-section, but occur less frequently by an exactly compensating amount. At lowest order in $\tilde{z}$, the curvature $\tilde{\Omega}$ is a product of second spatial derivatives of $\tilde{z}$. In a Gaussian framework, discussed in section 2.1 of [39], these derivatives are independent of both the wave slope and the vertical velocity. $\tilde{\Omega}$ can then be marginalized out and integrated. The surface slope at specular facets is fixed by the geometry, and expressed in terms of the observation angles $\psi$ and $\varphi$. There remains:

$$ \langle |S| \rangle^2 (f, r; x_0, y_0, \Delta + \eta) = \frac{\pi |R|^2}{\cos^4(\psi)} \quad (10) \\ \times \tilde{\Upsilon}_r^2 \left( r - \left[ h - \eta + \frac{\kappa}{2h} \rho_0^2 + \dot{r}_0 \delta t_{rr} \right] \right) \\ \times \int P(\partial_x \tilde{z} = \partial_x \tilde{z} \Big|_{sp}, \partial_y \tilde{z} = \partial_y \tilde{z} \Big|_{sp}, \tilde{w}) \\ \times \tilde{\Upsilon}_f^2 \left( f - 2 \frac{\tilde{w}}{\lambda} + 2 \frac{\dot{r}_0}{\lambda} \right) d\tilde{w}. $$

The pdf of $\tilde{w}$ is discussed at length in appendix A, where its expression is shown to read

$$ P(\tilde{w}, \psi, \varphi) = \frac{1}{\pi \sqrt{2\pi} \sigma_{\tilde{w}} \text{mss}_T} \exp\left(-\frac{\tan^2 \psi}{\text{mss}_T}\right) \\ \times \exp\left(-\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + U_{GD} \tan(\psi) \cos(\varphi - \varphi_{GD})]^2\right), $$

with mssT the total mean squared surface slope, and

$$ \mathbf{U}_{GD} = -\mathbf{M}_{ss}^{-1} \mathbf{msv} $$

the $\mathbf{U}_{GD}$ vector, with $U_{GD}$ and $\varphi_{GD}$ its magnitude and direction. The distribution for $\tilde{w}$ is Gaussian, but not centered, with variance

$$ \sigma_{\tilde{w}}^2 = m_{tt} - \mathbf{msv}^T \mathbf{M}_{ss}^{-1} \mathbf{msv}, $$

where $m_{tt}$, $\mathbf{M}_{ss}$ and $\mathbf{msv}$ are defined in terms of spectral moments of $\tilde{z}$ and its time and space derivatives. This variance is uniform on the sea surface, and is somewhat smaller than the full-surface variance of $\tilde{w}$, because part of the dispersion of $\tilde{w}$ is absorbed in the variations of the projection of $\mathbf{U}_{GD}$ on the radar line-of-sight. Introducing this expression in Eq. (10), one obtains:

$$ \langle |S| \rangle^2 (f, r; x_0, y_0, \Delta + \eta) = \frac{|R|^2}{\cos^4(\psi)} \frac{\exp\left(-\frac{\tan^2 \psi}{\text{mss}_T}\right)}{\text{mss}_T} \quad (11) \\ \times \tilde{\Upsilon}_r^2 \left( r - \left[ h - \eta + \frac{\kappa}{2h} \rho_0^2 + \dot{r}_0 \delta t_{rr} \right] \right) \\ \times \int \exp\left(-\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + U_{GD} \tan(\psi) \cos(\varphi - \varphi_{GD})]^2 \right) \\ \times \frac{1}{\sqrt{2\pi}\sigma_{\tilde{w}}} \tilde{\Upsilon}_f^2 \left( f - 2\frac{\tilde{w}}{\lambda} + 2\frac{\dot{r}_0}{\lambda} \right) d\tilde{w}. $$

With $\tilde{\Upsilon}_f^2$ known to be Gaussian, the last integral is the convolution of two Gaussian functions. It yields the Doppler Point Target Response function, accounting for the spreading effect of surface wave motion on the relative velocity of the satellite and the facet:

$$ \tilde{\Upsilon}_f^2(\delta f) = \left(\frac{25}{46}\right)^2 \frac{\sigma_f}{\tilde{\sigma}_f} e^{-\frac{\delta f}{2\tilde{\sigma}_f^2}}, \quad (12) $$

with

$$ \tilde{\sigma}_f^2 = \sigma_f^2 + \frac{4\sigma_{\tilde{w}}^2}{\lambda^2}. $$

Inserting this expression in Eq. (11), and introducing the instrument power radiation diagram G and the various factors required to account for the spherical divergence of the transmitted signal, the DDM contribution per unit transmitted power of the neighborhood of $(x_0, y_0, \Delta + \eta)$ finally reads

$$ \begin{aligned} \text{dDDM}(f, r; x_0, y_0, \Delta + \eta) = & \frac{\lambda^2 G^2(\theta, \varphi)}{(4\pi)^3 r_0^4} \sigma_0(\psi, \varphi) \text{dS} \quad (13) \\ & \times \tilde{\Upsilon}_r^2 \left( r - \left[ h - \eta + \frac{\kappa}{2h} \rho_0^2 + v_r \delta t_{rr} - \frac{y_0 v_t}{h} \delta t_{rr} \right] \right) \\ & \times \tilde{\Upsilon}_f^2 \left( f + \frac{2}{\lambda} \left[ v_r - \frac{y_0 v_t}{h} + \tan(\psi) U_{GD} \cos(\varphi - \varphi_{GD}) \right] \right). \end{aligned} $$

### H. Facet-level Model Discussion

Figure 3 presents the DDM obtained for 7 individual patches of the sea surface located every 3000 m along the satellite track, not accounting for (top graphs) and accounting for (bottom graphs) the impact of surface waves on the signal. $v_r$ has been set equal to 0 in this example.

A first striking feature is that the footprint of each patch is highly localized in the DDM: each sea surface patch produces only a very small footprint over a burst duration, and even if averaging over 7 consecutive bursts was performed to produce the DDM ($\frac{1}{20}$ s incoherent averaging time such as classicaly used), the signature of each patch would only describe a very short segment of the parabola. Note, the full parabolic signature, often shown in publications, results from the combination of many footprints of individual patches, and
<figure><img src="image_4.png" /><figcaption>Fig. 3. a) DDM contributions of 1 m² sea surface patches regularly spaced every 3000 m in the along-track direction between $y = -9000$ m and 9000 m, taking into account only the instrument point-target response. c) DDM contributions of 1 m² sea surface patches regularly spaced every 3000 m in the along-track direction between $y = -9000$ m and 9000 m, taking into account the smearing in the range and Doppler directions due to fully developed waves described by the spectral form [42], with a 12 m.s⁻¹ headwind (3.75 m significant wave height). b) and d): close-up views on the neighborhood of the $y = 9000$ m patch (red rectangles in subplots a) and c). On all the graphs, the color shades represent the DDM normalized by its maximal value. The dotted, dash-dotted, dashed and thick lines represent respectively the locus of the sea surface in the flat-Earth approximation, in the spherical-Earth approximation, in the spherical-Earth approximation, account taken of the $\delta t_{rr}$ term, and taking also into account the Geophysical Doppler term. The green line segments in subplots b) and d) represent the shift along the $f$ direction of the image of the patch in the interval separating one burst and the next and one burst and the sixth next (classical 20 Hz product). The thin white crosses in subplots b) and d) mark the tip of the Delay-Doppler PTR.</figcaption></figure>

the full DDM must be estimated as an integral over the sea surface.

A second feature, mostly visible in Fig. 3a and b, is that even in the “no-waves” case the frequency resolution achieved in a single burst is quite coarse in comparison with the total Doppler bandwidth (the ratio is equal to the number of inter-pulse intervals per burst, 63 in this case). The “Range Migration Correction” on which DDA processing is based, by applying a frequency-dependent range shift, leaks this Doppler spread into a large contribution to range spread. This can then lead to overestimation of the significant wave height during the geophysical parameter estimation step within waveform retracking. The effective widths of the range and Doppler point-target responses $\tilde{\Upsilon}_r^2$ and $\tilde{\Upsilon}_f^2$, estimated as a function of wind speed using the elevation spectrum model of [42] in its fully-developed seas limit, are represented in figure 4. The wave contribution to the effective range PTR width rapidly increases with wind speed, and becomes dominant for wind speed larger than 6 m.s⁻¹. This is consistent with the con- spicuous broadening of the PTR in the range direction, visible between figures 3b and 3d. By contrast, the instrumental contribution to the effective Doppler PTR width dominates up to wind speeds in excess of 13 m.s⁻¹. The overall Doppler PTR width only weakly depends on wind speed, remaining between 80 Hz, its low-wind value, and 120–130 Hz at high wind speed. Again, this is consistent with the only marginal width increase of the PTR in the Doppler direction between figures 3b and 3d.

Returning to Eq. (13), the effect of a non-zero $v_r$ on the DDM is minor. The two terms involving $v_r$ amount to shifts of the DDM; one in the range direction ($v_r \delta t_{rr}$ term in the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Effective range ($\tilde{\sigma}_r$, blue, left axis) and Doppler ($\tilde{\sigma}_f$, red, right axis) widths of the wave-accounting Point Target Responses $\tilde{\Upsilon}_r^2$ and $\tilde{\Upsilon}_f^2$ as a function of wind speed (bottom axis) and corresponding $H_s$ (top axis), for the elevation spectrum model [42] at infinite fetch. Dashed lines mark the instrumental lower bounds, thin lines mark the waves-only contribution, and thick lines represent the total effective widths.</figcaption></figure>

range PTR), and one in the Doppler direction ($2v_r/\lambda$ term in the Doppler PTR, leading to a “Doppler Centroid” correction). While not negligible, these two effects are readily corrected. In the following, we assume these corrections are applied and omit both terms.

More importantly, the various white lines in Fig. 3 mark
the locus of a flat ($\eta = 0$) sea surface on the DDM under a number of approximations. The dotted lines represent the flat-Earth approximation, which is clearly too coarse. The dash-dotted lines represent the spherical-Earth approximation, showing improved agreement with the observed behaviour, but still differing with respect to the center of the PTR in Fig. 3d. Adding the range-Doppler ambiguity correction $\delta t_{rr}r_0$ results in the dashed lines. While further improved, reaching near-perfect agreement also requires the $U_{GD}$ correction factor (solid line). This last correction should be applied as a stretch of the $f$ coordinate, but can be more conveniently expressed as a shift in the $r$ coordinate for given $f$: taking $x_0 = 0$ for simplicity, the contribution to the DDM observed at $f$ is the one that should appear at $f/(1 - \kappa \cos \varphi_{GD}U_{GD}/v_t)$ or, using the notations of Eq. (2), $f/(1 + \varepsilon)$. Its spherical Earth deflection correction $\Delta$ should thus not be evaluated as

$$ \Delta(f)|_{U_{GD}=0} = \frac{\kappa (fh\lambda)^2}{2h (2v_t)^2} $$

but as

$$ \Delta(f)|_{U_{GD}} = \frac{\kappa (fh\lambda)^2}{2h (2v_t)^2} \frac{1}{(1 + \varepsilon)^2} $$

The effect of $U_{GD}$ must thus be represented as a perturbation to the spherical Earth range deflection correction,

$$ \begin{aligned} \delta\Delta(f) &= \Delta(f)|_{U_{GD}} - \Delta(f)|_{U_{GD}=0} \\ &= \Delta(f)|_{U_{GD}=0} \left[ \frac{1}{(1 + \varepsilon)^2} - 1 \right] \\ &\approx -2\varepsilon\Delta(f)|_{U_{GD}=0} \end{aligned} $$

These corrections are represented as a function of the Doppler

<figure><img src="image_6.png" /><figcaption>Fig. 5. Range corrections as a function of Doppler frequency for a sea surface patch located on the satellite ground-track. The blue line (left scale) represents the spherical-Earth range deflection. The solid red line (right scale, note the very different scale) represents the range-Doppler ambiguity correction $r_0\delta t_{rr}$. The dashed red line (right scale) represents the $\delta\Delta$ correction for waves described by the spectral form [42], with a $12 \text{ m.s}^{-1}$ headwind and infinite fetch (3.75 m significant wave height).</figcaption></figure>

frequency in Fig. 5. The spherical Earth deflection correction (in blue) clearly dominates. Yet, the other corrections are far from being negligible in the very demanding context of precise nadir altimetry. The range-Doppler ambiguity correction ranges up to 15 cm at the ends of the unambiguous

along-track zone, and the $\delta\Delta$ correction accounting for the effect of the Geophysical Doppler contribution can reach several centimeters of magnitude in the $12 \text{ m.s}^{-1}$ situation (the discussion in Appendix A shows that this case is in fact representative of the frequently encountered ocean situations with $|U_{10}| \ge 7 \text{ m.s}^{-1}$).

## IV. NEAR-NADIR DELAY/DOPPLER RADAR: FROM DELAY/DOPPLER MAP TO ECHO WAVEFORM

Section III provided an in-depth development for the instrument response to an infinitesimal ocean surface patch. The upscaling operation from a single facet response to complete sea surface DDM - more specifically, the ways in which this response can be enfolded into the “triple convolution model” and production of the stacked echo waveform, are now addressed. All processing algorithms yield identical stacked waveforms in the along-track homogeneous situation. Selecting one or the other to construct our analysis leads to no restriction to its applicability other than the along-track homogeneity assumption. For computational convenience, this section essentially follows the rationale of the LR-RMC algorithm of [15].

### A. The Flat Surface Impulse Response

Eq. (13) expresses the DDM contribution of an infinitesimal ocean surface patch located at $(x_0, y_0, \Delta + \eta)$. Considering the sea surface to be flat, particularizing to the $v_r = 0$ case, assuming the antenna radiation diagram to be axisymmetric, the backscattering cross-section to be independent of $\theta$ and $\varphi$ in the narrow solid angle illuminated by the radar, approximating the $1/r_0^4$ term introduced by the spherical divergence by its value at nadir, integrating over the sea surface, denoting by primes the integration point coordinates, and finally introducing the expression (2) for the Doppler shift at the integration point, a simplified expression for the full DDM is:

$$ \begin{aligned} \text{DDM}(f, r) = & \frac{\lambda^2 \sigma_0}{(4\pi)^3 h^4} \iint_{Sea} G^2(\theta') \rho' d\rho' d\varphi' \\ & \times \widetilde{\Upsilon}_r^2 \left( r - h + \eta - \frac{\kappa}{2h} \rho'^2 + \rho' \cos(\varphi') \frac{v_t}{h} \delta t_{rr} \right) \\ & \times \widetilde{\Upsilon}_f^2 \left( f - \frac{2v_t}{\lambda h} \rho' (1 + \varepsilon) \cos(\varphi' - \varphi_\varepsilon) \right), \end{aligned} $$

It is a 2-dimensional convolution of a combined range and Doppler PTR that accounts for the effect of waves, both in terms of vertical displacement and Doppler broadening. This is performed with a Flat sea Surface Impulse Response (FSIR) that accounts for the Geophysical Doppler bias as:

$$ \begin{aligned} \text{FSIR}(f, r) = & \frac{\lambda^2 \sigma_0}{(4\pi)^3 h^4} \iint_{Sea} G^2(\rho'/h) \\ & \times \delta_r \left( r - h + \eta - \frac{\kappa}{2h} \rho'^2 + \rho' \cos(\varphi') \frac{v_t}{h} \delta t_{rr} \right) \\ & \times \delta_f \left( f - \frac{2v_t}{\lambda h} \rho' (1 + \varepsilon) \cos(\varphi' - \varphi_\varepsilon) \right) \rho' d\rho' d\varphi'. \end{aligned} $$

In this expression, $\delta_r$ and $\delta_f$ denote Dirac delta functions acting along the $r$ and $f$ directions.

At this point, a useful simplification is to replace the second integration coordinate, i.e. the azimuth $\varphi'$, by the Doppler
frequency at the integration point $\frac{2v_t}{\lambda h} \rho'(1 + \varepsilon) \cos(\varphi' - \varphi_\varepsilon)$. Next, we neglect the small misalignment angle $\varphi_\varepsilon$ and an order $\varepsilon$ correction in the range-Doppler ambiguity term. Introducing Eq. (1) for the antenna gain, one obtains a Doppler-distributed breakup of the FSIR as

$$ \text{FSIR}(f, r) = \frac{2\lambda^2 G_0^2 \sigma_0}{(4\pi)^3 h^4} \int_0^\infty \int_{-f_{\text{max}}}^{f_{\text{max}}} \delta_f (f - f') \frac{e^{-\frac{4}{h^2} \frac{\rho'}{\rho}}}{\sqrt{f_{\text{max}}^2 - f'^2}} \times \delta_r \left( r - h + \eta - \frac{\kappa}{2h} \rho'^2 + \frac{\lambda}{2} f' \delta t_{rr} \right) \rho' d\rho' d f', $$

with $f_{\text{max}}(\rho') = \frac{2v_t}{\lambda h}(1 + \varepsilon)\rho'$. Defining

$$ \mu_{\varepsilon} = \frac{\kappa h \lambda^2}{8v_t^2(1+\varepsilon)^2}, \quad \nu = \frac{8}{\gamma\kappa h}, \qquad (14) \\ A_{\varepsilon} = \frac{2\lambda^2 G_0^2 \sigma_0 \sqrt{\mu_{\varepsilon}}}{(4\pi)^3 \kappa h^3}, $$

and the frequency and range coordinates of the apex of the FSIR in the range/Doppler plane

$$ f_A = \frac{\lambda \delta t_{rr}}{4\mu_\varepsilon}, \quad r_A = h - \eta - \mu_\varepsilon f_A^2 $$

the FSIR becomes

$$ \text{FSIR}(f, r) = A_{\varepsilon} \frac{H(r - r_A - \mu_{\varepsilon}(f - f_A)^2)}{\sqrt{r - r_A - \mu_{\varepsilon}(f - f_A)^2}} e^{-\nu[r - r_A]} \quad (15) $$

with $H(\cdot)$ the Heaviside step function. Note that an additional small correction in the antenna diagram term has been neglected. The effect of a non-zero $\varepsilon$ in the expressions of $f_A$ and $r_A$ is assumed negligible. This FSIR expression is deemed to be suitable for following aspects of this study including the moving surface case.

Taking the S6-MF parameters, one finds $f_A = 8.447$ Hz, and $\mu_\varepsilon f_A^2 = 0.146$ mm. Though not fully negligible, these shifts can still be considered small and, being of technological origin, are deterministic and correctable. In the following, we consider them to be perfectly corrected ($r_A = h - \eta$, $f_A = 0$). This amounts to the assumption that the correction has been applied, as represented as a thick red line in Fig. 5.

## B. From the Flat Surface Impulse Response to the echo waveform

Going from the expressions of the Range/Doppler PTR and the FSIR to an expression for the DDA echo waveform requires the following steps:

* Computation of the Fourier transforms of the PTR and the FSIR along the range and Doppler dimensions.

* Multiplication of the Fourier transforms.

* Inverse Fourier transformation in the Doppler dimension.

* Multiplication by a Doppler-dependent mask to implement the Range Migration Correction.

* Integration in the Doppler dimension to perform the stacking operation.

* Inverse Fourier transformation in the range direction.

All but the last step can be analytically performed.

The first step is detailed in Appendices B-A and B-B. The Fourier transforms of the FSIR and the PTR are obtained as functions of two new independent variables $T$ (conjugate to $f$) and $K$ (conjugate to $r$) as

$$ \widehat{\text{FSIR}}(T, K) = A_0 \frac{\pi e^{-iK r_A} e^{-\frac{T^2}{4\mu_\varepsilon(\nu+iK)}}}{\sqrt{\mu_0}(\nu+iK)} $$

and

$$ \widehat{\widetilde{\Upsilon}_f^2 \Upsilon_r^2}(T, K) = \sqrt{2\pi}\sigma_f \left(\frac{25}{46}\right)^2 e^{-\frac{T^2 \sigma_f^2 + \kappa^2 \sigma_h^2}{2}} \widehat{\Upsilon_r^2}(K), $$

with $\widehat{\Upsilon_r^2}(K)$ either of the two Fourier-space representations Eq. (5) or Eq. (6) of the range PTR. Their product is:

$$ \widehat{\overline{\text{DDM}}}(T, K) = A_0 \frac{\pi \sqrt{2\pi} \sigma_f}{\sqrt{\mu_0}} \frac{25^2}{46^2} \frac{e^{-iK r_A - \frac{\kappa^2 \sigma_h^2}{2}}}{\nu + iK} \widehat{\Upsilon_r^2}(K) \times e^{-\frac{T^2}{2} \left[ \widehat{\sigma_f^2} + \frac{1}{2\mu_\varepsilon(\nu+iK)} \right]} \quad (16) $$

As a sanity check, this expression can be used to derive the CA echo waveform. This calculation, presented in Appendix B-C, indeed yields the expected results.

Noting this expression to be Gaussian in $T$, its inverse Fourier transform in this variable is

$$ \widehat{\text{DDM}}(f, K) = A_{\varepsilon} \sqrt{2\pi} \sigma_f \frac{25^2}{46^2} \frac{e^{-iK r_A - \frac{\kappa^2 \sigma_h^2}{2}}}{\sqrt{\nu + iK}} \widehat{\Upsilon_r^2}(K) \times \frac{e^{-\frac{\mu_\varepsilon(\nu+iK)f^2}{1+2\mu_\varepsilon(\nu+iK)\widehat{\sigma_f^2}}} {\sqrt{1+2\mu_\varepsilon(\nu+iK)\widehat{\sigma_f^2}}}. \quad (17) $$

The next processing stage is to multiply this function by a Doppler-dependent mask to apply the Range Migration Correction. The aim of this operation is to phase-shift the contributions at different $f$ in order to ensure a constructive summation in the stacking stage. Based on the original treatment by [10], the phase correction is $\exp(i\mu_0 K f^2)$. This step is discussed in Appendix B-D, and yields the DDA echo waveform in Fourier space as:

$$ \widehat{W_{DDA}}(K) = A_0 \frac{\sqrt{2\pi}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{46^2} \times \frac{e^{-iK r_A - \frac{\kappa^2 \sigma_h^2}{2}} \widehat{\Upsilon_r^2}(K)}{\sqrt{\nu + iK} \sqrt{\nu - 2iK(\varepsilon + \mu_0 \nu \widehat{\sigma_f^2}) + 2\mu_0 K^2 \widehat{\sigma_f^2}}} \quad (18) $$

Analytically inverting the Fourier transform to obtain the echo waveform as a function of $r$ is only straightforward for the case where $\varepsilon = 0$, $\widehat{\sigma_f} = 0$, and when using the Gaussian approximation of the range PTR. This is done in Appendix B-E, where the result by [21] is recovered. We have not been able to obtain an analytic form in a more general case.

## C. Full DDM Model Discussion

Figure 6a represents the FSIR in the $(f, r)$ plane taking account of the along-track sidelobes in the Doppler dimension (see section IV-D below). Restricting attention to the unambiguous region of the plane, the locus of the FSIR displays
<figure><img src="image_7.png" /><figcaption>Fig. 6. a) Graph of the Flat Surface Impulse Response Eq. (15) as a function of Doppler frequency shift $f$ and range $r - r_A$. FSIR replicas shifted by $+f_p$ and $-f_p$ have been added to represent the along-track sidelobes. b) Graph of the physical space DDM obtained by inverse-transforming the Fourier space expression Eq. (17) with $\bar{\sigma}_f^2 = 0$. This DDM takes into account the wave-induced broadening of the PTR in the range direction only, for 3.75 m significant wave height. Fourier-space DDM replicas shifted by $+f_p$ and $-f_p$ have been added before the inverse transform to represent the along-track sidelobes. c) Graph of the physical space DDM obtained by inverse-transforming the Fourier space expression Eq. (17). This DDM takes into account the wave-induced broadening of the PTR in the range and Doppler directions due to fully-developed seas described by the spectral form [42], with a $12 \text{ m.s}^{-1}$ headwind (3.75 m significant wave height). Fourier-space DDM replicas shifted by $+f_p$ and $-f_p$ have been added before the inverse transform to represent the along-track sidelobes. d) Graph of the Range Migration Corrected physical space DDM obtained by inverse-transforming the Fourier space expression Eq. (22), in the same geophysical conditions. Fourier-space DDM replicas shifted by $+f_p$ and $-f_p$ have been added before the RMC and inverse transform to represent the along-track sidelobes.</figcaption></figure>
the expected parabolic shape. It is slightly displaced due to the $U_{GD}$ effect, but this effect is too faint to be measurable in this representation. Taking slices for fixed $f$ (not shown), it varies from zero for $r-r_A \le \mu_\epsilon(f-f_A)^2$, and follows the expected $1/\sqrt{r-r_A-\mu_\epsilon(f-f_A)^2}$ singularity on crossing the locus. At each end of the interval, the $f > f_p/2$ (resp. $f < -f_p/2$) branch of the parabola is aliased into the $f \ge -f_p/2$ (resp. $f \le f_p/2$) part of the domain. Close to the edges, this induces a twofold increase of the recorded power density.

Fig. 6b shows the DDM in the $(f, r)$ plane, and under the $\tilde{\sigma}_f = 0$ Hz approximation for a $12 \text{ m.s}^{-1}$ headwind (3.75 m significant wave height). This case thus includes the effect of the instrument PTR and waves in the range direction only. It corresponds to conditions shown in [21]. As expected, the convolution by the wave-broadened instrument PTR introduces a strong smearing of the FSIR in the range direction.

Fig. 6c shows the DDM in the $(f, r)$ plane for $12 \text{ m.s}^{-1}$ headwind (3.75 m significant wave height), now taking account of both the instrument and waves PTR in the range and Doppler directions. The Doppler component of the PTR smears the DDM in the Doppler direction and strongly increases the thickness of the DDM. This effect is particularly noticeable for large values of $f$, where the FSIR locus slants more with respect to the iso-range lines. This leads to increased smearing in the $f$ direction.

Finally, Fig. 6d represents the range-migration corrected DDM in the $(f, r)$ plane and for the same environmental conditions. The unambiguous part of the DDM has successfully been flattened to the $r \approx r_A$ portion of the domain. In this portion of the DDM, the smearing in $f$ has been converted by the RMC into a strong excess of smearing in the $r$ direction, which can easily be confused with a strong excess of significant wave height.

The sidelobes, on the other hand, have been straightened (the parabolic part of the RMC applies for every $f$), but not folded to $r \approx r_A$. Their power density is thus missing from the $r \approx r_A$ region, and is contaminating the $r > r_A$ part of the domain. Thus the overall shape of the echo waveform is affected. This issue can analytically be modeled and is the subject of the next section.

## D. Accounting for Poseidon-4 Doppler Sidelobes

The Poseidon-4 open-burst acquisition mode imposes stringent constraints on the radar $f_p$: it must vary to accommodate changes in flight altitude along the orbit, and has to remain near the low value of 9 kHz. Consequently, signals scattered from distant up/down track areas can experience phase shifts of more than a half period between consecutive pulses. This introduces sidelobes in the DDM: a bin located at $(r, f)$ collects power that should appear at $(r, f)$, but also the power that should appear at $(r, f+f_p)$ and $(r, f-f_p)$. Higher-order ambiguities are rejected by the instrument radiation diagram and play a negligible role.

These sidelobes are inconsequential in the CA context, as the applied processing does not depend on $f$. In the DDA context, however, all echoes detected in the Doppler bin at $f$ are range-shifted by $\mu_0 f^2$. This includes those

which should have been shifted by $\mu_0(f + f_p)^2$ or by $\mu_0(f - f_p)^2$. The total stacked echo waveform is thus better reproduced by stacking $\overline{\text{DDM}}(f, \mathcal{K})e^{i\mu_0\mathcal{K}f^2}$ for $f \in [-\frac{f_p}{2}; \frac{f_p}{2}]$, $\overline{\text{DDM}}(f, \mathcal{K})e^{i\mu_0\mathcal{K}(f-f_p)^2}$ for $f \in [\frac{f_p}{2}; \frac{3f_p}{2}]$ and $\overline{\text{DDM}}(f, \mathcal{K})e^{i\mu_0\mathcal{K}(f+f_p)^2}$ for $f \in [-\frac{3f_p}{2}; -\frac{f_p}{2}]$. The expression for the full DDM including the sidelobes is given in Appendix B-F Eq. (23).

After DCA correction, the DDM is effectively even in $f$. The contributions from the two sidelobes are thus identical. The sidelobes-accounting Fourier space DDA waveform can thus be expressed as

$$ \begin{aligned} \widehat{W}_{DDA}^{SL}(\mathcal{K}) = & 2 \int_0^{\frac{f_p}{2}} e^{i\mu_0 \mathcal{K} f^2} \overline{\text{DDM}}(f, \mathcal{K}) \mathrm{d}f & (19) \\ & + 2 \int_{\frac{f_p}{2}}^{\frac{3f_p}{2}} e^{i\mu_0 \mathcal{K}(f-f_p)^2} \overline{\text{DDM}}(f, \mathcal{K}) \mathrm{d}f. & \end{aligned} $$

The calculation is detailed in Appendix B-G. Finally, the Fourier-space and sidelobe-accounting DDA waveform is obtained as

$$ \begin{aligned} \widehat{W}_{DDA}^{SL}(\mathcal{K}) = & A_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{462} & (20) \\ & \times \frac{e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}} \hat{\Upsilon}_r^2(\mathcal{K})}{\sqrt{\nu + i\mathcal{K}} \sqrt{\nu - 2i\mathcal{K}(\varepsilon + \mu_0\nu\tilde{\sigma}_f^2) + 2\mu_0\mathcal{K}^2\tilde{\sigma}_f^2}} \\ & \times \left[ \text{erf}\left(\frac{f_p \Xi_\varepsilon}{2}\right) \right. \\ & \qquad \left. + e^{-f_p^2 \frac{\Xi_\varepsilon}{4}} \text{erfc}\left(f_p \left[\frac{\Xi_\varepsilon}{2} + \frac{i\mu_0\mathcal{K}}{\Xi_\varepsilon}\right]\right) \right. \\ & \qquad \left. - e^{-f_p^2 \left[i2\mathcal{K}\mu_0 + \frac{9\Xi_\varepsilon}{4}\right]} \text{erfc}\left(f_p \left[\frac{3\Xi_\varepsilon}{2} + \frac{i\mu_0\mathcal{K}}{\Xi_\varepsilon}\right]\right) \right]. \end{aligned} $$

with $\overline{\text{erfc}}(x) = e^{x^2}\text{erfc}(x)$ the exponentially scaled complementary error function,

$$ \Xi_{\varepsilon} \approx \sqrt{\mu_{\varepsilon} \frac{\nu - 2i\mathcal{K}(\varepsilon + \mu_0\nu\tilde{\sigma}_f^2) + 2\mu_0\mathcal{K}^2\tilde{\sigma}_f^2}{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\tilde{\sigma}_f^2}}, $$

and we recall for the sake of compactness that

$$ \begin{aligned} A_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} &= \frac{\lambda^2 G_0^2 \sigma_0}{(8\pi)^{3/2} \kappa h^3}, \\ \varepsilon &= -\kappa \frac{U_{GD}}{v_t} \cos(\varphi_{GD}), \\ \mu_{\varepsilon} &= \frac{\kappa h \lambda^2}{8v_t^2(1+\varepsilon)^2}, \quad \nu = \frac{8}{\gamma \kappa h}, \\ \delta t_{rr} &= \frac{h}{c} + \frac{f_c}{s}, \quad f_A = \frac{\lambda \delta t_{rr}}{4\mu_{\varepsilon}}, \\ r_A &= h - \eta - \mu_{\varepsilon} f_A^2, \\ \sigma_r &= \frac{c}{2B}, \quad \sigma_f = \frac{1}{\sqrt{2\log(2)}} \frac{1.293}{2T_b}, \\ \tilde{\sigma}_f^2 &= \sigma_f^2 + \frac{4\sigma_{\overline{w}}^2}{\lambda^2}, \\ \hat{\Upsilon}_r^2(\mathcal{K}) &= \sigma_r \mathbf{1}_{[-1;1]} \left(\frac{\mathcal{K}\sigma_r}{2\pi}\right) \left[1 - \frac{|\mathcal{K}|\sigma_r}{2\pi}\right]. \end{aligned} $$
This analytically-derived final expression is the IASCO DDA
waveform model.

V. DISCUSSION

A. *Comparison between different waveform models*

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

Fig. 7. Graph of altimeter waveforms computed for a 12m.s⁻¹ crosswind, corresponding to Hs = 3.75 m, ε = 0, for the S6-MF flight parameters and for a number of different altimeter processing configurations. The cyan curve represents the CA waveform Eq. (21). The red curve represents the σ̃ = 0 DDA waveform discussed in Appendix B-E, under the assumption of very large fₚ. The green curve represents the DDA waveform Eq. (18), which takes into account the Doppler PTR width, but assumes large fₚ. The dashed blue curve represents the unambiguous Doppler frequency range [−fₚ/2; fₚ/2] contribution to the IASCO waveform model Eq. (20), and the continuous blue curve represents the full IASCO waveform model Eq. (20), including the effect of Doppler sidelobes. Finally, the magenta curves represent the multilook waveforms produced by the “PySAMOSA” [31] implementation of the established SAMOSA model [22] over the main lobe only (dashed magenta curve), and after modification to account for the sidelobe aliasing (continuous magenta curve). All analytical curves correspond to the same A₀ factor, chosen such that the σ̃ = 0 waveform is one at the epoch. The SAMOSA waveforms have been normalized to correspond to the same integrated energy as the blue waveforms.

Fig. 7 presents waveforms obtained using the different expressions developed in the text, for the case of a Hs = 3.75 m sea surface generated using the spectrum of [42] for 12 m.s⁻¹ cross wind and infinite fetch. All the waveforms are generated using the same value of A₀ for normalization. They contain identical integrated energy, except for the dashed blue and dashed magenta curves, which only account for the returns from the unambiguous Doppler frequency band. Aside from this particular case, the curves only differ by the expected distribution of the received energy as a function of time.

The cyan curve, which represents the expected waveform
from CA processing, has a very different shape from the
others: after the sharp initial rise at the encounter of the radar
pulse with the surface, the power return in all range bins
is expected to be essentially identical, except for the slow
decrease due to the antenna gain pattern.

On the contrary, DDA processing brings a much larger
portion of the received radar energy into the leading edge.
All DDA model waveforms reflect this general behaviour. The

significant DDA differences shown reflect the different sets of
hypotheses under which these waveforms are derived.

Typically, the waveform derived by [21] displays a much
sharper initial rise. That model accounts for the spreading
effect in range of the vertical distribution of backscattering
facets, but not for the relatively coarse frequency resolution
of the along-track Fourier transform and the broadening due
to the vertical orbital velocities (red curve). Its peak level is
also much higher than the other DDA waveforms. Accounting
for the additional smearing of the energy in the Doppler
dimension leads to a markedly different waveform depicted
by the green line. The main peak is largely reduced, and there
is a significantly earlier rise in energy (the “toe” mentioned
in [21]).

This effect of the finite Doppler resolution of the along-track Fourier transform is taken into account by the SAMOSA model [22] for multi-look SAR altimeter returns, as well as by the earlier semi-analytical models [19], [20], [23]. For these modeled waveforms, a correction serendipitously emerges from numerical integrations performed over azimuth. Yet the broadening effect of scattering facet motion is not explicitly taken into account by any of these models.

Figure 6 clearly shows that the halved $f_p$ value of the S6-MF instrument with respect to the CryoSAT and Sentinel-3A instruments introduces aliasing in the along-track Fourier transform. Consequently, only the energy returning from the strip of the sea surface with a Doppler frequency shift comprised between $-f_p/2$ and $+f_p/2$ is correctly redistributed by the RMC step. The waveform obtained by stacking over this reduced frequency band is represented as a blue dashed line in Figure 7. The height of the main peak is reduced by almost 20% with respect to the green curve, obtained by extending the integration to $\pm\infty$. Taking into account the energy contained in the sidelobes (continuous blue curve) only marginally increases the height of the main peak, but slows the decrease of the signal tail. We stress that the total energy content of this waveform is identical to that of the very different-looking CA (cyan) or $\tilde{\sigma}_f = 0$ (red) waveforms.

Waveforms produced using the well-established SAMOSA model [22] are depicted by magenta lines, either multi-looking over the unambiguous Doppler frequency range (dashed line) or, after *ad hoc* modification to take into account the aliasing effect and over the full frequency range (continuous line). In both cases, the SAMOSA waveforms feature a higher main peak and a somewhat reduced “toe” compared to their counterparts produced using Eq. (20) (blue lines). The waveforms rapidly match thereafter, showing similar behavior over the decay region. Overall, the IASCO waveform model Eq. (20) is in good agreement with SAMOSA over the late peak and decay region. Differences in the leading edge are however significant and are discussed next.

B. Effect of the Gaussian approximation of the range PTR

The red curve in Figure 8 represents the waveform obtained
using the modified multilook SAMOSA model in the same en-
vironmental conditions, together with the waveforms obtained
with Eq. (20) using the sinc² (blue curve) and Gaussian (black
<figure><img src="image_9.png" /><figcaption>Fig. 8. Graph of altimeter waveforms computed for a 12m.s⁻¹ cross wind, corresponding to Hs = 3.75 m, ε = 0, for the S6-MF flight parameters. The continuous blue curve represents the waveform obtained using the full DDA waveform Eq. (20), assuming the instrument range PTR has a squared sinc shape Eq. (5). The continuous black curve represents the waveform obtained in the same conditions but now using for a range PTR having a Gaussian approximation Eq. (6). The dashed blue curve represents the waveform obtained with the squared sinc shape, this time neglecting the broadening induced by the surface wave motions (ie in the $\bar{\sigma}_f = \sigma_f$ approximation). The dashed black curve represents the waveform obtained with the Gaussian approximation of the PTR and in the $\bar{\sigma}_f = \sigma_f$ approximation. Finally, the red curve represent the multilook waveforms produced using the “PySAMOSA” [31] implementation of the established SAMOSA model [22] after modification to account for sidelobe aliasing. The normalization is consistent with Figure 7. Zoomed-in plots of the waveforms at their peak and “toe” are shown in the two insets on the left of the figure.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Graph of normalized altimeter waveforms for different 10-m neutral wind speeds, with significant wave height and vertical velocity variance computed from the spectrum [42] in its infinite-fetch limit. All waveforms are computed with the full DDA waveform Eq. (20), considering for the instrument range PTR the exact squared sinc Eq. (5). The wind direction is orthogonal to the satellite track ($ε = 0$). The influence of wind on signal level through variations of the backscattering cross-section $σ_0$ is not accounted for.</figcaption></figure>

curve) approximations of the instrument range PTR. The effect of the Gaussian approximation is visible with a slight (~4%) increase in peak height, and a very slightly reduced “toe” with respect to the more exact squared sinc representation. Though this effect is clearly visible, it is not sufficient to explain the discrepancy with the SAMOSA waveform.

The dashed black curve represents the waveform obtained in the Gaussian range PTR approximation (consistent with SAMOSA), but neglecting the Doppler broadening due to surface motion (also consistent with SAMOSA). The models are now in close agreement.

This result is one validation of the IASCO waveform model Eq. (20), but it also hints at a strong impact of the scattering facet motion on the SAR altimeter waveform. As the difference between the continuous and dashed blue curves shows, neglecting this effect likely leads to a strong overestimation of the waveform peak height, as well as a noticeable underestimation of the waveform “toe” extension.

As a side remark, we notice that the Gaussian approximation of the instrument range PTR, which is not practically useful since the final Fourier transform is performed numerically anyway, actually shows a slight detrimental effect. Henceforth we do not use it.

### C. Sensitivity to environmental conditions

1. **Dependence on $U_{10}$ at infinite fetch:** The dependence of the IASCO waveform model Eq. (20) is represented as a function of 10-m neutral wind speed in Figure 9. The wind direction is orthogonal to the satellite track ($ε = 0$). The wind only acts through its influence on $H_s$ and $\sigma_{\tilde{w}}$, which are computed using the spectrum prescription of [42] in its infinite-fetch limit. The dependence of the peak height on wind through the backscatter cross-section $σ_0$ is not accounted for.

With increasing winds, the initial rise of the waveform progressively broadens, and the peak height progressively decreases. The leading edge half-power point also progressively shifts leftward. The decay part of the waveform is essentially impervious to parameter changes.

2. **Dependence on $H_s$ for fixed $\sigma_{\tilde{w}}$:** Figure 9 mixes the influences of $H_s$ and $\sigma_{\tilde{w}}$, the two parameters affecting the waveform. Figure 10 shows the effect of changing the significant wave height for a fixed value of $\sigma_{\tilde{w}}$, kept equal to its value in an infinite-fetch 12 m.s⁻¹ $U_{10}$ sea surface, 0.77 m.s⁻¹. $ε$ is kept equal to 0.

With increasing $H_s$, the peak height decreases, and the initial rise progressively slackens. The half-height point gradually shifts to earlier instants.

3. **Dependence on $\sigma_{\tilde{w}}$ for fixed $H_s$:** Figure 11 now conversely presents the effect of varying $\sigma_{\tilde{w}}$ for a fixed $H_s$, kept equal to its 12 m.s⁻¹ value of 3.75 m (refer to Figure 4 for intuition on the range of variations explored here). Clearly, the influence of $\sigma_{\tilde{w}}$ on the waveform is quite different from the influence of $H_s$: with increasing $\sigma_{\tilde{w}}$, the peak height decreases, the waveform “toe” lengthens, but this is achieved without shifting the half-height point of the peak.

Results imply a mechanism that is likely to bias the joint estimation of $H_s$ and $r_A$ from the waveform: in cases of particularly large $\sigma_{\tilde{w}}$, one might imagine the peak of the
<figure><img src="image_11.png" /><figcaption>Fig. 10. Graph of normalized altimeter waveforms for fixed $\sigma_{\tilde{w}} = 0.77 \text{ m.s}^{-1}$ and varying $H_s$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 12. Difference between SAMOSA retracked $H_s$ and “true” $H_s$ as a function of “true” $H_s$, obtained by varying $U_{10}$ in the infinite-fetch limit of spectrum [42].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Graph of normalized altimeter waveforms for fixed $H_s = 3.75 \text{ m}$ and varying $\sigma_{\tilde{w}}$.</figcaption></figure>

1) Significant height: Figure 12 presents the difference between “true” $H_s$ and $H_s$ estimates obtained running the SAMOSA retracker on IASCO waveforms generated with varying wind speeds. Again, the infinite-fetch form of the spectra is used to estimate $H_s$, $\sigma_{\tilde{w}}$ and $U_{GD}$ as functions of $U_{10}$. Clearly, the produced waveforms are consistently associated by SAMOSA to overestimated values of $H_s$. This bias varies rapidly for small $H_s$, and increases weakly afterwards, rising from 0.22 m for $H_s \approx 2 \text{ m}$ to 0.29 m for $H_s \approx 10.5 \text{ m}$. For a given wind speed, this bias is identical in the headwind, tailwind and crosswind situations, showing it is essentially independent of $\varepsilon$, and dominated by the effect of $\sigma_{\tilde{w}}$.

waveform to be quite low and the midpoint to be unaffected, which a tracker might mistake for a large $H_s$ situation with a depressed sea surface.

4) Effect of the “Geophysical Doppler” shift ($\varepsilon \neq 0$): The effect of a non-zero $\varepsilon$ is extremely small, and would be hard to capture in large-scale graphs such as Figures 9, 10 and 11. To provide insight into its influence on the waveform, the most practical approach is to study the biases its presence induces in the results of a waveform retracker such as SAMOSA. This is the subject of the next section.

<figure><img src="image_14.png" /><figcaption>Fig. 13. SAMOSA retracked Sea Level for a nominally flat surface located at $\eta = 0$, obtained by varying $U_{10}$ in the infinite-fetch limit of spectrum [42].</figcaption></figure>

## D. Tracker biases

The influences of $U_{GD}$ and $\sigma_{\tilde{w}}$ on the results of the SAMOSA retracker are studied in this section by generating IASCO waveforms using known parameter values and then examining biases in the retrieved values.

2) Sea Level: Fig. 13 presents the Sea Level (SL) estimates obtained using the SAMOSA retracker on IASCO waveforms generated for varying wind speeds blowing over a flat surface located at $\eta = 0$. Clearly, the biasing effect mentioned above
is present: a bias varying between −1 cm and 2 cm is present even for zero $\varepsilon$ (black curve). The mere existence of a non-zero $\sigma_{\tilde{w}}$ induces a bias in the retrieved SL. This time, however, the effect of a non-zero “Geophysical Doppler” shift, though very small, is clearly visible. The case of a headwind drives the tracker to overestimate SL, while the existence of a tail wind has the opposite effect. The azimuthal dependence of the bias (not shown) is essentially sinusoidal, but is not centered around 0, as the $\sigma_{\tilde{w}}$ effect introduces a (actually larger) isotropic bias.

### E. Effect of instrument parameters

The IASCO waveform can be applied with no modifications to assess CryoSAT-2 and Sentinel-3A cases. But because the altitude, flight velocity, and $f_p$ of these instruments all differ, it is non-trivial to disentangle parameter change impacts from simple waveform plot assessments. This section thus focuses on waveform sensitivity to the most important parameters, the pulse repetition frequency and the platform flight altitude. Two other key parameters are the antenna beamwidth and the flight velocity. The antenna beamwidth affects the waveform in a straightforward way, by changing the slope of the decay region of the waveform (the finer the beamwidth, the faster the decay). The flight velocity is not a very sensitive parameter over its quite limited range of variation, and mainly affects the balance between the main lobe and sidelobes of the along-track SAR processing - this by changing the portion of the instrument FoV that is associated to Doppler shifts smaller than the Nyquist frequency.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Graph of altimeter waveforms obtained for a 12m.s⁻¹ side wind, corresponding to Hs = 3.75 m, $\sigma_{\tilde{w}}$ = 0.77 m.s⁻¹, $\varepsilon$ = 0, for the S6-MF flight parameters and for instrument $f_p$ equal to varying fractions of the actual S6-MF $f_p$.</figcaption></figure>

1) Effect of instrument $f_p$: Fig. 14 shows the dependence of the IASCO waveform as a function of the instrument $f_p$, expressed in fraction of the actual S6-MF $f_p$, for constant burst duration (i.e. the number of pulses per burst is increased when $f_p$ increases). Reducing the instrument $f_p$ by a factor of two (cyan curve) reduces the effectiveness of the RMC and reduces the main peak height by a large factor with

respect to the nominal (orange) curve. Interestingly, increasing the $f_p$ by a factor of two (green curve), effectively makes S6-MF similar to CryoSAT-2 and Sentinel-3A $f_p$. This also seems to almost saturate this effect. Further increasing $f_p$ only brings minimal changes (the black curve represents the $f_p \to \infty$ limit). This is expected, as for large values of $f_p$ the ambiguous frequency ranges correspond to far-off distance in the along-track. Signal that is effectively rejected by the off-nadir antenna gain pattern.

<figure><img src="image_16.png" /><figcaption>Fig. 15. Graph of altimeter waveforms obtained for a 12m.s⁻¹ side wind, corresponding to Hs = 3.75 m, $\sigma_{\tilde{w}}$ = 0.77 m.s⁻¹, $\varepsilon$ = 0, for the S6-MF flight parameters and $f_p$, and for varying instrument $N_b$ (burst duration).</figcaption></figure>

2) Effect of instrument $N_b$: Conversely, Fig. 15 shows the dependence, for fixed $f_p$, of the IASCO waveform on the burst duration (number of pulses per burst). The burst duration influences the waveform through the along-track SAR processing resolution $\sigma_f$. Changing $N_b$ changes the measurement floor on $\sigma_{\tilde{w}}$, shown as the red dashed line in Figure 4. Higher $N_b$ values reduce the instrumental contribution to the overall $\tilde{\sigma}_f$. As can be seen in the figure, even though $N_b$ = 64, the value implemented in the S6-MF onboard processor, may be a bit marginal, the result obtained for a value of $N_b$ = 128 (green curve) is very close to the ideal curve obtained in the limit of vanishing instrumental contribution (black curve), at least in these conditions. Increasing this much further could increase the influence of facet range walk during the burst time window. Also, the finite lifetime of the individual scatters, which is currently not well known, sets limits to gains that can be achieved by increasing the coherent processing time window. Reprocessing S6-MF data with $N_b$ = 128 is feasible on ground, but can not be done using the current onboard processor.

3) Effect of platform altitude: Fig. 16 shows the dependence of the IASCO waveform on the platform altitude $h$, compensated for the expected $h^{-5/2}$ power decrease [10]. The remaining dependence after amplitude compensation is small, and is essentially due to changes in energy balance between the main lobe and the sidelobes of the along-track SAR processing. It almost disappears in the $\tilde{\sigma}_f \to 0$ limit. The
<figure><img src="image_17.png" /><figcaption>Fig. 16. Graph of altimeter waveforms obtained for a 12m.s⁻¹ side wind, corresponding to Hs = 3.75 m, ε = 0, for varying platform altitude, the S6-MF fₚ and constant 7200 m.s⁻¹ platform speed. A compensation factor of (h/h₀)⁵/² has been applied.</figcaption></figure>

only difference remains in the decay region. It is consistent with the change in the aliasing lobe slope in the DDM caused by the change of μ. In this case, the pulse repetition frequency is fixed and equal to its S6-MF value.

## VI. CONCLUSION

With the ambition to more precisely estimate the impact of ocean surface motions on Delay-Doppler altimetry, a new analytical (Fourier-space in range) Delay Doppler Altimeter waveform has been derived.

Developments start with the DDM signature of an isolated backscattering facet of the sea surface, the basic building block from which the aggregated DDM response of the sea surface is composed. Through a careful analysis of the joint statistics of its equivalent radar cross-section, its elevation with respect to the mean sea level, and its instantaneous velocity with respect to the radar instrument, the ensemble average of this signature is obtained under the Gaussian statistical approximation.

The convolution of this signature with the Flat Surface Impulse Response is computed, yielding an analytical expression of the full DDM, Eq. (23). Integrating this expression in Doppler, the IASCO waveform, an analytical expression of the DDA waveform, is finally obtained, Eq. (20). This waveform is validated against the well-established SAMOSA model waveform, modified to account for the SAR aliasing sidelobes.

Due to its analytical form, IASCO model waveforms and their relative sensitivities can be readily obtained with respect to change in varied technological and environmental parameters. In particular, the influence of surface facet orbital and mean motion on the retracking process has been studied, where the SAMOSA retracker is used as a benchmark. Surface phenomena not taken into account by the SAMOSA waveform model, but addressed by the IASCO waveform model, are found to potentially cause observable biases in geophysical

parameters of crucial scientific importance, such as significant wave height and altimeter range (and hence Sea Level).

Our study highlights the sensitivity of the DDA radar waveform to motion-related characteristics of the sea surface, which on the contrary have no influence on the CA radar waveform. In the context set by the open-burst acquisition mode of the S6-MF Poseidon-4 instrument, this opens the way to a coupled DDA and CA analysis of the radar observations, which will provide for a high-precision, unbiased continuation of the historical timeseries of σ₀, Hs and Sea Level, supplemented with new timeseries of kinematic quantities such as σᵢ and Uₐ:

* One option might be to combine the SAMOSA and IASCO waveforms in a retracker to retrieve and compare directly, through a Maximum-Likelihood approach, estimates of surface wave vertical motion standard deviation σᵢ and of the along-track projection of the “Geophysical Doppler” vector Uₐ.

* These new parameters could be derived and find scientific use in their own right, or be used as inputs to altimeter EMB- or SSB-mitigation routines.

* A hybrid approach might also be attempted, using CA estimates, which are impervious to the motion-induced biases DDA studied here, as first guesses or constraints for the DDA retracker. It will lead to the implementation of a Maximum A Posteriori optimization approach ingesting both the DDA and CA instrument waveforms.

* Another possible option would be to apply the fitting procedure to the DDM itself. Our Eq. (23) could serve as the basis for such an attempt.

Future investigations will be conducted to refine these theoretical and analytical formulations to provide a tool for present and future 2D off-nadir altimeter design and performance assessment. Efforts will also be conducted to optimize and extend DDA waveform processing strategies (e.g. [43]) to retrieve upper ocean sea surface height slopes and motions in an effort to better estimate upper ocean velocities.

## APPENDIX A

$$ \text{EXPRESSING } P(\partial_x \tilde{z} = \partial_x \tilde{z}|_{sp}, \partial_y \tilde{z} = \partial_y \tilde{z}|_{sp}, \tilde{w}). $$

From section 2.1 in [39], an expression can be simply derived for P(∂x̃z = ∂x̃z|sp, ∂ỹz = ∂ỹz|sp, w̃) in terms of spectral moments of the sea surface elevation. In the notations of [30], [44],

$$ m_{tt} = \langle \partial_t \tilde{z} \partial_t \tilde{z} \rangle, \quad \mathbf{msv} = \begin{vmatrix} \langle \partial_t \tilde{z} \partial_x \tilde{z} \rangle \\ \langle \partial_t \tilde{z} \partial_y \tilde{z} \rangle \end{vmatrix}, \\ \mathbf{M}_{ss} = \begin{bmatrix} \langle \partial_x \tilde{z} \partial_x \tilde{z} \rangle & \langle \partial_x \tilde{z} \partial_y \tilde{z} \rangle \\ \langle \partial_x \tilde{z} \partial_y \tilde{z} \rangle & \langle \partial_y \tilde{z} \partial_y \tilde{z} \rangle \end{bmatrix} $$

the correlation matrix of surface vertical velocity and surface slope components is:

$$ \mathbf{M}_{\text{tss}} = \left[ \begin{array}{c|c} \frac{m_{tt}}{\mathbf{msv}} & \mathbf{msv}^T \\ \hline \mathbf{msv} & \mathbf{M}_{\text{ss}} \end{array} \right] $$
and their joint pdf :

$$P(\partial_x \tilde{z}, \partial_y \tilde{z}, \tilde{w}) = \frac{1}{(2\pi)^{3/2} \sqrt{|\det(\mathbf{M}_{\text{tss}})|}} \times \exp \left( -\frac{1}{2} [\tilde{w}, \partial_x \tilde{z}, \partial_y \tilde{z}] \mathbf{M}_{\text{tss}}^{-1} \begin{vmatrix} \tilde{w} \\ \partial_x \tilde{z} \\ \partial_y \tilde{z} \end{vmatrix} \right).$$

The pdf of $\tilde{w}$ at $\partial_x \tilde{z} = \partial_x \tilde{z}|_{sp}$, $\partial_y \tilde{z} = \partial_y \tilde{z}|_{sp}$ is obtained by rearranging this expression. Using the Sherman-Morrison-Woodbury matrix inversion formula, the block inverse of $\mathbf{M}_{\text{tss}}$ reads:

$$\mathbf{M}_{\text{tss}}^{-1} = \frac{1}{\sigma_{\tilde{w}}^2} \left[ \frac{1}{-\mathbf{M}_{\text{ss}}^{-1} \mathbf{m}\mathbf{sv}} \begin{vmatrix} -\mathbf{m}\mathbf{sv}^T \mathbf{M}_{\text{ss}}^{-1} \\ \mathbf{C}_{\text{ss}} \end{vmatrix} \right],$$

with

$$\sigma_{\tilde{w}}^2 = m_{tt} - \mathbf{m}\mathbf{sv}^T \mathbf{M}_{\text{ss}}^{-1} \mathbf{m}\mathbf{sv}, \\ \mathbf{C}_{\text{ss}} = \sigma_{\tilde{w}}^2 \mathbf{M}_{\text{ss}}^{-1} + \mathbf{M}_{\text{ss}}^{-1} \mathbf{m}\mathbf{sv} \mathbf{m}\mathbf{sv}^T \mathbf{M}_{\text{ss}}^{-1}.$$

Also:

$$\det(\mathbf{M}_{\text{tss}}) = \sigma_{\tilde{w}}^2 \det(\mathbf{M}_{\text{ss}}).$$

One recognizes in $-\mathbf{M}_{\text{ss}}^{-1} \mathbf{m}\mathbf{sv}$ the “Geophysical Doppler” vector $\mathbf{U}_{\text{GD}}$ [30]. This vector is composed of the Total Surface Current Vector and a (usually dominant) “Wave Doppler” vector. The magnitude of $\mathbf{U}_{\text{GD}}$ and the azimuth towards which it points are respectively denoted as $\mathbf{U}_{\text{GD}}$ and $\varphi_{\text{GD}}$. Figure 17 represents $\mathbf{U}_{\text{GD}}$, $\sqrt{m_{tt}}$ and $\sigma_{\tilde{w}}$ as a function of

<figure><img src="image_18.png" /><figcaption>Fig. 17. Spectral moments of $\partial_t \tilde{z}$ as a function of wind speed, for the elevation spectrum model [42] at infinite fetch, expressed in m.s$^{-1}$ (left axis) and in Hz (right axis).</figcaption></figure>

wind speed, using the elevation spectrum model of [42] for infinite fetch, and under the assumption of infinitely short electromagnetic wavelength (see discussion at the end of this section). Already discussed in [30], [44], $\mathbf{U}_{\text{GD}}$ varies rapidly with wind speed up to $|U_{10}| \sim 7$ m.s$^{-1}$, after which it levels off on a slowly increasing plateau. The other quantities, $\sqrt{m_{tt}}$ and $\sigma_{\tilde{w}}$, increase almost monotonously, $\sqrt{m_{tt}}$ remaining slightly larger than $\sigma_{\tilde{w}}$.

The pdf of $\tilde{w}$ for a given surface slope $\nabla \tilde{z}|_{sp}$ is:

$$P(\tilde{w}, \nabla \tilde{z}_{sp}) = \frac{1}{(2\pi)^{3/2} \sigma_{\tilde{w}} \sqrt{|\det(\mathbf{M}_{\text{ss}})|}} \times \exp \left( -\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + \mathbf{U}_{\text{GD}} \cdot \nabla \tilde{z}|_{sp}]^2 \right) \times \exp \left( -\frac{1}{2} \nabla \tilde{z}|_{sp}^T \mathbf{M}_{\text{ss}}^{-1} \nabla \tilde{z}|_{sp} \right).$$

The term on the first line is a normalization factor, the term on the second line shows that the mode of the distribution of $\tilde{w}$ is displaced by the projection of $-\mathbf{U}_{\text{GD}}$ along the radar line-of-sight, and the final term accounts for the decreasing probability of finding specular facets away from nadir. The surface elevation gradient $\nabla \tilde{z}|_{sp}$ on specular facets at the current point is fixed by the observation geometry, and such that $\partial_x \tilde{z}|_{sp} = \tan(\psi) \sin(\varphi)$, and $\partial_y \tilde{z}|_{sp} = \tan(\psi) \cos(\varphi)$. The $\mathbf{U}_{\text{GD}} \cdot \nabla \tilde{z}|_{sp}$ term thus reduces to $\mathbf{U}_{\text{GD}} \tan(\psi) \cos(\varphi - \varphi_{\text{GD}})$. The surface slope covariance matrix $\mathbf{M}_{\text{ss}}$, being symmetric, can be diagonalized by a mere coordinate system rotation, and a maximal surface slope variance azimuth $\varphi_{\text{ss}}$ can be found such that:

$$\mathbf{M}_{\text{ss}} = \mathbf{R}_{\varphi_{\text{ss}}} \times \begin{bmatrix} \mathbf{m}_a & 0 \\ 0 & \mathbf{m}_x \end{bmatrix} \times \mathbf{R}_{-\varphi_{\text{ss}}},$$

with $\mathbf{R}_{\varphi_{\text{ss}}}$ the matrix describing rotation of angle $\varphi_{\text{ss}}$ in the $(x, y)$ plane, and $\mathbf{m}_a$ and $\mathbf{m}_x$ the mean squared surface slopes in the “along” and “across” directions. With these notations, one obtains

$$P(\tilde{w}, \nabla \tilde{z}_{sp}) = \frac{1}{(2\pi)^{3/2} \sigma_{\tilde{w}} \sqrt{\mathbf{m}_a \cdot \mathbf{m}_x}} \times \exp \left( -\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + \mathbf{U}_{\text{GD}} \tan(\psi) \cos(\varphi - \varphi_{\text{GD}})]^2 \right) \times \exp \left( -\frac{\tan^2 \psi}{2} \left[ \frac{\cos^2(\varphi - \varphi_{\text{ss}})}{\mathbf{m}_a} + \frac{\sin^2(\varphi - \varphi_{\text{ss}})}{\mathbf{m}_x} \right] \right).$$

As in [9], we denote the total mean squared surface slope by $\text{mss}_T = \mathbf{m}_a + \mathbf{m}_x$. The fractional difference between $\mathbf{m}_a$ and $\mathbf{m}_x$ is described by $\delta_{ss}$, such that $\mathbf{m}_a = \frac{\text{mss}_T}{2}(1 + \delta_{ss})$ and $\mathbf{m}_x = \frac{\text{mss}_T}{2}(1 - \delta_{ss})$. With these notations

$$P(\tilde{w}, \nabla \tilde{z}_{sp}) = \frac{1}{\pi \sqrt{2\pi} \sigma_{\tilde{w}} \text{mss}_T \sqrt{1 - \delta_{ss}^2}} \times \exp \left( -\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + \mathbf{U}_{\text{GD}} \tan(\psi) \cos(\varphi - \varphi_{\text{GD}})]^2 \right) \times \exp \left( -\frac{\tan^2 \psi}{\text{mss}_T(1 - \delta_{ss}^2)} [1 - \delta_{ss} \cos(2(\varphi - \varphi_{\text{ss}}))] \right).$$

The integrals yielding the mean squared slope statistics as moments of the surface elevation spectrum logarithmically diverge with wavevector $k$: they are thus dependent on the structure of the sea surface at all scales. This ultraviolet divergence can be regularized either by the roll-off of the elevation spectrum in the viscous regime, or by the filtering effect caused by the finiteness of the electromagnetic wavelength, depending on which limit is reached first. Numerous studies introduced $\text{mss}_{\text{shape}}$ to account for the electromagnetic filtering effect.
In the optical limit, [45]–[47] show $m_x/m_a$ to be weakly dependent on environmental conditions and close to 0.6. In the present notations, this translates to $\delta_{ss} \simeq 0.25$.

Admittedly, the relevance of these optical results to Ku-band physics remains to be firmly established, but point to a small impact of directional effects on the backscattering cross-section in the near-nadir geometry [48]. Taking the $\delta_{ss} \to 0$ limit in the above expression:

$$ P(\tilde{w}, \nabla \tilde{z}_{sp}) = \frac{1}{\pi \sqrt{2\pi} \sigma_{\tilde{w}} \text{mss}_T} \exp\left(-\frac{\tan^2 \psi}{\text{mss}_T}\right) \\ \times \exp\left(-\frac{1}{2\sigma_{\tilde{w}}^2} [\tilde{w} + U_{GD} \tan(\psi) \cos(\varphi - \varphi_{GD})]^2\right). $$

Note one should probably replace mssT by mssshape to empirically account for the electromagnetic cut-off effect. In our calculations, integrations are performed up to the viscous cut-off limit, using the spectral form proposed by [42], to keep mssT in the developments.

# APPENDIX B

## FROM FSIR TO ECHO WAVEFORM

### A. Fourier transformation of the FSIR

Starting from Eq. (15), the Fourier transform in $r$ and $f$ of the FSIR can be expressed as

$$ \widehat{\text{FSIR}}(\mathcal{T}, \mathcal{K}) = A_{\varepsilon} e^{-i[\mathcal{T}f_A + \mathcal{K}r_A]} \times \iint d\mathbf{r} d\mathbf{r}' e^{-i\mathcal{T}f - r(\nu + i\mathcal{K})} \frac{H(r - \mu_{\varepsilon}f^2)}{\sqrt{r - \mu_{\varepsilon}f^2}}. $$

Performing the change of variable $\sqrt{r - \mu_{\varepsilon}f^2} \to u$ turns this expression into

$$ \widehat{\text{FSIR}}(\mathcal{T}, \mathcal{K}) = A_{\varepsilon} e^{-i[\mathcal{T}f_A + \mathcal{K}r_A]} \times \int d\mathbf{r} e^{-i\mathcal{T}f - (\nu + i\mathcal{K})\mu_{\varepsilon}f^2} \int e^{-u^2(\nu + i\mathcal{K})} du. $$

The two Gaussian integrations can be performed independently. Introducing the approximations $\mu_{\varepsilon} \simeq \mu_0(1 - 2\varepsilon)$ and $A_{\varepsilon} \simeq A_0(1 - \varepsilon)$, one obtains the final result as

$$ \widehat{\text{FSIR}}(\mathcal{T}, \mathcal{K}) = A_0 \frac{\pi e^{-i[\mathcal{T}f_A + \mathcal{K}r_A]} e^{-\frac{\mathcal{T}^2}{4\mu_{\varepsilon}(\nu + i\mathcal{K})}}}{\sqrt{\mu_0}(\nu + i\mathcal{K})}. $$

### B. Fourier transformation of the wave-accounting Point Target Response

Starting from the expressions (9) and (12) of the wave-accounting PTRs, and using the fact that

$$ \int du e^{-iu\mathcal{K}} e^{-\frac{u^2}{2\sigma^2}} = \sqrt{2\pi\sigma} e^{-\frac{\kappa^2\sigma^2}{2}}, $$

the Fourier transform of the PTR is easily obtained as

$$ \widehat{\widehat{\Upsilon}_f^2 \Upsilon_r^2}(\mathcal{T}, \mathcal{K}) = \sqrt{2\pi} \sigma_f \frac{25^2}{46^2} e^{-\frac{\mathcal{T}^2 \sigma_f^2 + \kappa^2 \sigma_h^2}{2}} \widehat{\Upsilon_r^2}(\mathcal{K}). $$

Depending on the intended use of the expression, either the exact representation of the instrument PTR Eq. (5) or its Gaussian approximation Eq. (6) can be used for $\widehat{\Upsilon}_r^2(\mathcal{K})$.

### C. Computation of the CA echo waveform from Equation (16)

Going from expression (16) to the DDA waveform still requires significant work. The succession of steps is detailed in the main text: the Fourier-space DDM must still be inverse-transformed in the Doppler dimension, the RMC must be applied, the result must be stacked by integration in the Doppler dimension, and the result must be inverse-transformed in the range dimension to yield the waveform.

In the CA context, however, no RMC is applied before the integration. The range Fourier-transform of the CA waveform is equal to the Doppler-dimension integral of the DDM, which is easily obtained as the value of Eq. (16) along the $\mathcal{T} = 0$ line:

$$ \widehat{W_{CA}}(\mathcal{K}) = A_0 \frac{\sqrt{2\pi\pi\sigma_f} 25^2}{\sqrt{\mu_0} 46^2} e^{-i\mathcal{K}r_A} \frac{e^{-\frac{\kappa^2\sigma_h^2}{2}}}{\nu + i\mathcal{K}} \widehat{\Upsilon_r^2}(\mathcal{K}). $$

Assuming the Gaussian shape for $\widehat{\Upsilon}_r^2$, the inverse Fourier transform of this function can be obtained using Eq. 3.462.3 of [49] in terms of the parabolic cylinder function $D_{-1}$. This function can in turn be expressed in terms of the complementary error function, and the result is obtained as:

$$ W_{CA}(r) = A_0 \frac{\pi^2 \sigma_r \sigma_f 25^2}{\sqrt{\mu_0} 46^2} e^{-\frac{\nu^2 \sigma_r^2}{2}} \times \exp\left(-\nu \left[r - r_A - \nu \widehat{\sigma}_r^2\right]\right) \operatorname{erfc}\left(-\frac{r - r_A - \nu \widehat{\sigma}_r^2}{\sqrt{2} \widehat{\sigma}_r}\right). \quad (21) $$

This expression coincides with the classical result [3], [21]. The expression of $r_A$ contains a $\frac{\delta t_r^2 v_r^2}{2\kappa h}$ shift, which amounts to 0.146 mm: due to the range-Doppler ambiguity bias, the sea surface appears 0.146 mm closer to the satellite than it actually is. The effective interface thickness $\widehat{\sigma}_r$, accounting for both the instrument PTR and the sea surface height dispersion due to waves, has its usual expression.

### D. Range Migration Correction and stacking

In this step, the expression of $\widehat{\text{DDM}}(f, \mathcal{K})$ given by Eq. (17) has to be multiplied by $\exp(i\mu_0 \mathcal{K} f^2)$, and the result must be integrated ("stacked") in $f$. The product of the DDM Fourier transform and the RMC mask is:

$$ \widehat{\text{DDM}}_{\text{RMC}}(f, \mathcal{K}) = A_{\varepsilon} \sqrt{2\pi} \sigma_f \frac{25^2}{46^2} \frac{e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}}}{\sqrt{\nu + i\mathcal{K}}} \widehat{\Upsilon}_r^2(\mathcal{K}) \times \frac{e^{-f^2 \left[ \frac{\mu_{\varepsilon}(\nu + i\mathcal{K})}{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2} - i\mathcal{K}\mu_0 \right]}}{\sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}}. $$

Introducing $\Xi_{\varepsilon}$ such that

$$ \Xi_{\varepsilon}^2 = \frac{\mu_{\varepsilon}(\nu + i\mathcal{K})}{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2} - i\mathcal{K}\mu_0, $$

it can be expressed as

$$ \widehat{\text{DDM}}_{\text{RMC}}(f, \mathcal{K}) = A_{\varepsilon} \sqrt{2\pi} \sigma_f \frac{25^2}{46^2} \frac{e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}}}{\sqrt{\nu + i\mathcal{K}}} \widehat{\Upsilon}_r^2(\mathcal{K}) \times \frac{e^{-\Xi_{\varepsilon}^2 f^2}}{\sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}}. \quad (22) $$
Integrating this expression in $f$ yields the Fourier-transform in $r$ of the DDA waveform. The integral is Gaussian and is easily performed, and the result is obtained as

$$ \begin{aligned} \widehat{W}_{DDA}(\mathcal{K}) ={}& \mathcal{A}_{\varepsilon} \sqrt{2\pi^{3/2}} \sigma_f \frac{25^2 e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}}}{46^2 \sqrt{\nu + i\mathcal{K}}} \widehat{\Upsilon}_r^2(\mathcal{K}) \\ & \times \frac{1}{\Xi_{\varepsilon}\sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}}. \end{aligned} $$

This expression now deserves being simplified. Using again the approximations $\mu_{\varepsilon} \simeq \mu_0(1-2\varepsilon)$ and $\mathcal{A}_{\varepsilon} \simeq \mathcal{A}_0(1-\varepsilon)$, and neglecting terms of order $\varepsilon^2$, one obtains

$$ \Xi_{\varepsilon} \approx \sqrt{\mu_{\varepsilon} \frac{\nu - 2i\mathcal{K}(\varepsilon + \mu_0 \nu \widehat{\sigma}_f^2) + 2\mu_0 \mathcal{K}^2 \widehat{\sigma}_f^2}{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}} $$

and

$$ \begin{aligned} \widehat{W}_{DDA}(\mathcal{K}) ={}& \mathcal{A}_0 \sqrt{2\pi^{3/2}} \sigma_f \frac{25^2 e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}}}{46^2 \sqrt{\mu_0 \sqrt{\nu + i\mathcal{K}}}} \widehat{\Upsilon}_r^2(\mathcal{K}) \\ & \times \frac{1}{\sqrt{\nu - 2i\mathcal{K}[\varepsilon + \mu_0 \nu \widehat{\sigma}_f^2] + 2\mu_0 \mathcal{K}^2 \widehat{\sigma}_f^2}}. \end{aligned} $$

### E. Analytical inversion in the $\varepsilon = 0, \widehat{\sigma}_f = 0$ case

Inverting this Fourier transform in the general case is difficult. Assuming a Gaussian range PTR, the special case $\varepsilon = 0, \widehat{\sigma}_f = 0$ studied by [21] is however tractable, and provides a check of the expression. Starting from

$$ W_{DDA}(r) = \mathcal{A}_0 \frac{\pi \sigma_r \sigma_f 25^2}{\sqrt{\mu_0 \nu} 46^2} \int d\mathcal{K} \frac{e^{i\mathcal{K}(r-r_A)} e^{-\frac{\kappa^2 \widehat{\sigma}_r^2}{2}}}{\sqrt{\nu + i\mathcal{K}}}, $$

introducing the integration variable shift $\mathcal{K} \rightarrow \mathcal{K} + i\nu$ brings the integral to the form

$$ W_{DDA}(r) = \mathcal{A}_0 \frac{\pi \sigma_r \sigma_f 25^2}{\sqrt{\nu \mu_0} 46^2} e^{-\nu(r-r_A-\nu \frac{\widehat{\sigma}_r^2}{2})} \int (i\mathcal{K})^{-1/2} e^{i\mathcal{K}(r-r_A-\nu \widehat{\sigma}_r^2)} e^{-\frac{\kappa^2 \widehat{\sigma}_r^2}{2}} d\mathcal{K}. $$

Using Eq. 3.462.3 of [49] then yields

$$ W_{DDA}(r) = \mathcal{A}_0 \sigma_r \sigma_f \frac{25^2}{46^2} \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\widehat{\sigma}_r \nu \mu_0}} D_{-1/2} \left( -\frac{r-r_A-\nu \widehat{\sigma}_r^2}{\widehat{\sigma}_r} \right) \exp \left( -\frac{\nu^2 \widehat{\sigma}_r^2}{2} - \nu(r-r_A-\nu \widehat{\sigma}_r^2) - \frac{(r-r_A-\nu \widehat{\sigma}_r^2)^2}{4 \widehat{\sigma}_r^2} \right) $$

where $D_{-1/2}$ is the parabolic cylinder function of order $-1/2$.

Up to changes in the normalization constant and the transposition from $r$ to two-way travel time, this expression is identical to Eq. (A16) of [21].

As noted by these authors, numerical implementations of this parabolic cylinder function have a tendency to overflow for negative argument. This problem can be circumvented using Eqs. (19.3.7, 19.3.8, 19.27.4 and 19.27.5) from [50] to express it in terms of modified Bessel functions as:

$$ D_{-1/2}(x) = \sqrt{\frac{|x|}{2\pi}} K_{1/4}\left(\frac{x^2}{4}\right) + H(-x) \sqrt{\pi|x|} I_{1/4}\left(\frac{x^2}{4}\right). $$

A more uniform balancing of the different terms of the expression can then be achieved using the exponentially scaled implementations of the modified Bessel functions.

### F. Accounting for sidelobes in the DDM

For $f \in [-f_p/2; f_p/2]$, the full DDM, accounting for Doppler sidelobes, can be expressed as:

$$ \widehat{DDM}_{\text{RMC}}^{\text{SL}}(f, \mathcal{K}) = e^{i\mu_0 \mathcal{K} f^2} \times \left[ \widehat{DDM}(f, \mathcal{K}) + \widehat{DDM}(f+f_p, \mathcal{K}) + \widehat{DDM}(f-f_p, \mathcal{K}) \right] $$

Starting from Eq. (17), this is easily expressed as:

$$ \begin{aligned} \widehat{DDM}_{\text{RMC}}^{\text{SL}}(f, \mathcal{K}) ={}& \mathcal{A}_{\varepsilon} \sqrt{2\pi} \sigma_f \frac{25^2}{46^2} \widehat{\Upsilon}_r^2(\mathcal{K}) \\ & \times \frac{e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}} e^{-\Xi_{\varepsilon}^2 f^2}}{\sqrt{\nu + i\mathcal{K}} \sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}} \\ & \times \left[ 1 + e^{-(\Xi_{\varepsilon}^2 + i\mu_0 \mathcal{K})f_p(f_p+2f)} + e^{-(\Xi_{\varepsilon}^2 + i\mu_0 \mathcal{K})f_p(f_p-2f)} \right]. \end{aligned} \quad (23) $$

### G. Accounting for sidelobes in the DDA waveform

Starting back from Eqs. (19) and (17), the main lobe contribution to the sidelobes-accounting DDA waveform can be expressed as

$$ 2 \int_{0}^{\frac{f_p}{2}} e^{i\mu_0 \mathcal{K} f^2} \widehat{DDM}(f, \mathcal{K}) df = 2^{3/2} \mathcal{A}_{\varepsilon} \pi \sigma_f \frac{25^2}{46^2} \widehat{\Upsilon}_r^2(\mathcal{K}) \\ \times \frac{e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}}}{\sqrt{\nu + i\mathcal{K}} \sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}} \int_{0}^{f_p/2} e^{-\Xi_{\varepsilon}^2 f^2} df. $$

The integral can be expressed in terms of the error function. After simplification, this leads to

$$ 2 \int_{0}^{\frac{f_p}{2}} e^{i\mu_0 \mathcal{K} f^2} \widehat{DDM}(f, \mathcal{K}) df = \mathcal{A}_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{46^2} \\ \times \frac{\widehat{\Upsilon}_r^2(\mathcal{K}) e^{-i\mathcal{K}r_A - \frac{\kappa^2\sigma_h^2}{2}} \operatorname{erf}\left(\frac{f_p}{2}\right)}{\sqrt{\nu + i\mathcal{K}} \sqrt{\nu - 2i\mathcal{K}(\varepsilon + \mu_0 \nu \widehat{\sigma}_f^2) + 2\mu_0 \mathcal{K}^2 \widehat{\sigma}_f^2}}. $$

The contribution of the sidelobes, on the other hand, is

$$ 2 \int_{\frac{f_p}{2}}^{\frac{3f_p}{2}} e^{i\mu_0 \mathcal{K} (f-f_p)^2} \widehat{DDM}(f, \mathcal{K}) df = 2^{3/2} \mathcal{A}_{\varepsilon} \pi \sigma_f \frac{25^2}{46^2} \\ \times \frac{\widehat{\Upsilon}_r^2(\mathcal{K}) e^{-\frac{\kappa^2 \sigma_h^2}{2} + i\mathcal{K}[\mu_0 f_p^2 - r_A]}}{\sqrt{\nu + i\mathcal{K}} \sqrt{1 + 2\mu_{\varepsilon}(\nu + i\mathcal{K})\widehat{\sigma}_f^2}} \\ \times \int_{\frac{f_p}{2}}^{\frac{3f_p}{2}} e^{-\Xi_{\varepsilon}^2 f^2 - i2\mu_0 \mathcal{K} f f_p} df. $$
The integral can again be expressed in terms of the complementary error function erfc. The calculation is tedious, but straightforward, and the result is obtained as:

$$2 \int_{f_p/2}^{3f_p/2} e^{i\mu_0 K(f-f_p)^2} \overline{\overline{\text{DDM}}}(f, K) df = A_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{46^2} \times \frac{\hat{\Upsilon}_r^2(K)e^{-\frac{\kappa^2 \sigma_h^2}{2} + iK[\mu_0 f_p^2 - r_A]} - \frac{\kappa^2 \mu_0^2 f_p^2}{\Xi_\epsilon^2}}{\sqrt{\nu + iK}\sqrt{\nu - 2iK(\epsilon + \mu_0 \nu \hat{\sigma}_f^2) + 2\mu_0 K^2 \hat{\sigma}_f^2}} \times \left[ \text{erfc}\left(f_p\left[\frac{\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) - \text{erfc}\left(f_p\left[\frac{3\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) \right].$$

Finally, the complete Fourier space expression of the sidelobes-accounting DDA waveform is obtained as

$$\begin{aligned} \widehat{W}_{DDA}^{SL}(K) = & A_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{46^2} \hat{\Upsilon}_r^2(K) \\ & \times \frac{e^{-iKr_A - \frac{\kappa^2 \sigma_h^2}{2}}}{\sqrt{\nu + iK}\sqrt{\nu - 2iK(\epsilon + \mu_0 \nu \hat{\sigma}_f^2) + 2\mu_0 K^2 \hat{\sigma}_f^2}} \\ & \times \left[ \text{erf}\left(\frac{f_p \Xi_\epsilon}{2}\right) \right. \\ & \qquad \left. + e^{iK\mu_0 f_p^2 - \kappa^2 \frac{\mu_0^2 f_p^2}{\Xi_\epsilon^2}} \text{erfc}\left(f_p\left[\frac{\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) \right. \\ & \qquad \left. - e^{iK\mu_0 f_p^2 - \kappa^2 \frac{\mu_0^2 f_p^2}{\Xi_\epsilon^2}} \text{erfc}\left(f_p\left[\frac{3\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) \right]. \end{aligned}$$

Again, numerical implementations of the complementary error function have a tendency to overflow, and one has to resort to exponentially scaled implementations. Denoting $\overline{\text{erfc}}(x) = e^{x^2} \text{erfc}(x)$, a more stable expression is obtained as

$$\begin{aligned} \widehat{W}_{DDA}^{SL}(K) = & A_0 \frac{\sqrt{2\pi^{3/2}}}{\sqrt{\mu_0}} \sigma_f \frac{25^2}{46^2} \hat{\Upsilon}_r^2(K) \\ & \times \frac{e^{-iKr_A - \frac{\kappa^2 \sigma_h^2}{2}}}{\sqrt{\nu + iK}\sqrt{\nu - 2iK(\epsilon + \mu_0 \nu \hat{\sigma}_f^2) + 2\mu_0 K^2 \hat{\sigma}_f^2}} \\ & \times \left[ \text{erf}\left(\frac{f_p \Xi_\epsilon}{2}\right) \right. \\ & \qquad \left. + e^{-f_p^2 \frac{\Xi_\epsilon^2}{4}} \overline{\text{erfc}}\left(f_p\left[\frac{\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) \right. \\ & \qquad \left. - e^{-f_p^2 \left[i2K\mu_0 + \frac{9\Xi_\epsilon^2}{4}\right]} \overline{\text{erfc}}\left(f_p\left[\frac{3\Xi_\epsilon}{2} + \frac{i\mu_0 K}{\Xi_\epsilon}\right]\right) \right]. \end{aligned}$$

## ACKNOWLEDGMENT

This study was supported by the European Space Agency, through the Drift4SKIM and IASCO contracts (ESA contracts 4000126110/18/NL/FF/gp and 4000129945/19/NL/FF/gp, respectively).

It is a pleasure to acknowledge fruitful discussions with C. J. Donlon and A. Egido, from ESA, and A. Bohé, F. Boy and C. Maraldi, from CNES.

## REFERENCES

[1] G. D. Quartly, M. A. Srokosz, and A. C. McMillan, "Analyzing altimeter artifacts: Statistical properties of ocean waveforms," *Journal of Atmospheric and Oceanic Technology*, vol. 18, no. 12, pp. 2074 – 2091, 2001. [Online]. Available: https://doi.org/10.1175/1520-0426(2001)018(2074:AAASPO)2.0.CO;2

[2] D. E. Barrick, "Remote sensing of sea state by radar," in *Remote Sensing of the Troposphere*, V. Derr, Ed. Washington, D.C.: U.S. Govt. Printing Office, 1972, ch. 12, pp. 1–46.

[3] G. S. Brown, "The average impulse response of a rough surface and its applications," *IEEE Trans. Antennas Propag.*, vol. 25, no. 1, pp. 67–74, 1977. [Online]. Available: https://doi.org/10.1109/TAP.1977.114153

[4] D. B. Chelton, E. J. Walsh, and J. L. MacArthur, "Pulse compression and sea level tracking in satellite altimetry," *Journal of Atmospheric and Oceanic Technology*, vol. 6, no. 3, pp. 407 – 438, 1989. [Online]. Available: https://www.doi.org/10.1175/1520-0426(1989)006(0407:PCASLT)2.0.CO;2

[5] M. De Carlo, F. Ardhuin, A. Ollivier, and A. Nigou, "Wave groups and small scale variability of wave heights observed by altimeters," *Journal of Geophysical Research: Oceans*, vol. 128, no. 8, p. e2023JC019740, 2023. [Online]. Available: https://doi.org/10.1029/2023JC019740

[6] P. Thibaut, J. C. Poisson, E. Bronner, and N. Picot, "Relative performance of the MLE3 and MLE4 retracking algorithms on Jason-2 altimeter waveforms," *Marine Geodesy*, vol. 33, no. sup1, pp. 317–335, 2010. [Online]. Available: https://doi.org/10.1080/01490419.2010.491033

[7] B. Chapron, V. Kerbaol, D. Vandemark, and T. Elfouhaily, "Importance of peakedness in sea surface slope measurements," *Journal of Geophysical Research*, vol. 105, no. C7, pp. 17195–17202, 2000. [Online]. Available: https://doi.org/10.1029/2000JC000079

[8] D. Vandemark, B. Chapron, J. Sun, G. H. Crescenti, and H. C. Graber, "Ocean wave slope observations using radar backscatter and laser altimeters," *Journal of Physical Oceanography*, vol. 34, pp. 2825–2842, 2004. [Online]. Available: https://doi.org/10.1175/JPO2663.1

[9] 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 Geosci. Remote Sens. Lett.*, vol. 13, no. 9, pp. 1310–1314, 2016. [Online]. Available: https://doi.org/10.1109/LGRS.2016.2583198

[10] R. K. Raney, "The delay/Doppler radar altimeter," *IEEE Trans. Geosci. Remote Sens.*, vol. 36, no. 5, pp. 1578–1588, 1998. [Online]. Available: https://doi.org/10.1109/36.718861

[11] F. Boy, J.-D. Desjonquères, N. Picot, T. Moreau, and M. Raynal, "CryoSat-2 SAR-mode over oceans: Processing methods, global assessment, and benefits," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 1, pp. 148–158, 2017. [Online]. Available: https://doi.org/10.1109/TGRS.2016.2601958

[12] A. Egido and W. H. F. Smith, "Fully focused SAR altimetry: Theory and applications," *IEEE Geosci. Remote Sens. Lett.*, vol. 55, no. 1, pp. 392–406, 2017. [Online]. Available: https://doi.org/10.1109/TGRS.2016.2607122

[13] M. Kleinherenbrink, M. Naeije, C. Slobbe, A. Egido, and W. Smith, "Observations of polar ice fields," *Remote Sensing of Environment*, vol. 237, p. 111589, 2020. [Online]. Available: https://doi.org/10.1016/j.rse.2019.111589

[14] F. Collard, L. Marié, F. Nouguier, M. Kleinherenbrink, F. Ehlers, and F. Ardhuin, "Wind-wave attenuation in Arctic sea ice: A discussion of remote sensing capabilities," *Journal of Geophysical Research*, vol. 127, p. e2022JC018654, 2022. [Online]. Available: https://doi.org/10.1029/2022JC018654

[15] 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, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S027311772030911X

[16] C. Gommenginger, C. Martin-Puig, L. Amarouche, and R. K. Raney, "Review of state of knowledge for SAR altimetry over ocean, version 2.2," *EUMETSAT*, Tech. Rep., November 2013. [Online]. Available: https://nora.nerc.ac.uk/id/eprint/507775/1/SARAltimetry_Review_JasonCS_EUMETSAT.pdf

[17] R. K. Raney, "Maximizing the intrinsic precision of radar altimetric measurements," *IEEE Geoscience and Remote Sensing Letters*, vol. 10, no. 5, pp. 1171–1174, 2013. [Online]. Available: https://doi.org/10.1109/LGRS.2012.2235138
[18] 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, 2021. [Online]. Available: https://doi.org/10.1016/j.rse.2021.112395

[19] D. Wingham, L. Phalippou, C. Mavrocordatos, and D. Wallis, "The mean echo and echo cross product from a beamforming interferometric altimeter and their application to elevation measurement," *IEEE Trans. Geosci. Remote Sens.*, vol. 42, no. 10, pp. 2305-2323, 2004. [Online]. Available: https://doi.org/10.1109/TGRS.2004.834352

[20] A. Halimi, C. Mailhes, J.-Y. Tourneret, P. Thibaut, and F. Boy, "A semi-analytical model for Delay/Doppler Altimetry and its estimation algorithm," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 7, pp. 4248-4258, 2014. [Online]. Available: https://doi.org/10.1109/TGRS.2013.2280595

[21] E. S. Garcia, D. T. Sandwell, and W. H. Smith, "Retracking CryoSat-2, Envisat and Jason-1 radar altimetry waveforms for improved gravity field recovery," *Geophysical Journal International*, vol. 196, no. 3, pp. 1402-1422, 01 2014. [Online]. Available: https://doi.org/10.1093/gji/ggt469

[22] C. Ray, C. Martin-Puig, M. P. Clarizia, G. Ruffini, S. Dinardo, C. Gommenginger, and J. Benveniste, "SAR altimeter backscattered waveform model," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 2, pp. 911-919, 2015. [Online]. Available: https://doi.org/10.1109/TGRS.2014.2330423

[23] L. Recchia, M. Scagliola, D. Giudici, and M. Kuschnerus, "An accurate semianalytical waveform model for mispointed SAR interferometric altimeters," *IEEE Geosci. Remote Sens. Lett.*, vol. PP, pp. 1-5, 7 2017. [Online]. Available: https://doi.org/10.1109/LGRS.2017.2720847

[24] C. Buchhaupt, L. Fenoglio-Marc, S. Dinardo, R. Scharroo, and M. Becker, "A fast convolution based waveform model for conventional and unfocused SAR altimetry," *Advances in Space Research*, vol. 62, no. 6, pp. 1445-1463, 2018. [Online]. Available: https://doi.org/10.1016/j.asr.2017.11.039

[25] C. Buchhaupt, L. Fenoglio, M. Becker, and J. Kusche, "Impact of vertical water particle motions on focused SAR altimetry," *Advances in Space Research*, vol. 68, no. 2, pp. 853-874, 2021. [Online]. Available: https://doi.org/10.1016/j.asr.2020.07.015

[26] 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, 2023. [Online]. Available: https://doi.org/10.1016/j.asr.2022.12.034

[27] C. K. Buchhaupt, A. Egido, D. Vandemark, W. H. F. Smith, L. Fenoglio, and E. Leuliette, "Towards the mitigation of discrepancies in sea surface parameters estimated from low- and high-resolution satellite altimetry," *Remote Sensing*, vol. 15, no. 17, 2023. [Online]. Available: https://doi.org/10.3390/rs15174206

[28] A. Egido, C. Buchhaupt, F. Boy, C. Maraldi, E. Cadier, S. Dinardo, E. Leuliette, and T. Moreau, "A significant wave height correction to account for vertical wave motion effects in SAR altimeter measurements," in *Proceedings of the 2022 Ocean Surface Topography Science Team Meeting, Lido, Italy, 31 October - 4 November 2022*, 2022. [Online]. Available: https://doi.org/10.24400/527896/a03-2022.3460

[29] ESA, "Report for mission selection: SKIM," European Space Agency, Noordwijk, The Netherlands, Tech. Rep. ESA-EOPSM-SKIM-RP-3550, May 2019. [Online]. Available: https://esamultimedia.esa.int/docs/EarthObservation/EE9-SKIM-RfMS-ESA-v1.0-FINAL.pdf

[30] L. Marié, F. Collard, F. Nouguier, L. Pineau-Guillou, D. Hauser, F. Boy, S. Méric, P. Sutherland, C. Peureux, G. Monnier, B. Chapron, A. Martin, P. Dubois, C. J. Donlon, T. Casal, and F. Ardhuin, "Measuring ocean total surface current velocity with the KuROS and KaRADOC airborne near-nadir Doppler radars: a multi-scale analysis in preparation for the SKIM mission," *Ocean Science*, vol. 16, no. 6, pp. 1399-1429, 2020. [Online]. Available: https://os.copernicus.org/articles/16/1399/2020/

[31] F. Schlembach and M. Passaro, "PySAMOSA: An open source software framework for retracking SAMOSA-based, open ocean and coastal waveforms of SAR satellite altimetry," 2023. [Online]. Available: https://doi.org/10.5281/zenodo.8242782

[32] R. K. Raney, "CryoSat SAR-mode looks revisited," *IEEE Geoscience and Remote Sensing Letters*, vol. 9, no. 3, pp. 393-397, 2012. [Online]. Available: https://doi.org/10.1109/LGRS.2011.2170052

[33] ——, "Doppler properties of radars in circular orbits," *International Journal of Remote Sensing*, vol. 7, no. 9, pp. 1153-1162, 1986. [Online]. Available: https://doi.org/10.1080/01431168608948916

[34] B. Chapron, F. Collard, and F. Ardhuin, "Direct measurements of ocean surface velocity from space: interpretation and validation," *Journal of Geophysical Research*, vol. 110, no. C07008, 2005. [Online]. Available: https://doi.org/10.1029/2004JC002809

[35] R. Kodis, "A note on the theory of scattering from an irregular surface," *IEEE Trans. Antennas Propag.*, vol. 14, no. 1, pp. 77-82, 1966. [Online]. Available: https://doi.org/10.1109/TAP.1966.1138626

[36] D. Barrick, "Rough surface scattering based on the specular point theory," *IEEE Trans. Antennas Propag.*, vol. 16, no. 4, pp. 449-454, 1968. [Online]. Available: https://doi.org/10.1109/TAP.1968.1139220

[37] M. S. Longuet-Higgins, "Reflection and refraction at a random moving surface. I. pattern and paths of specular points," *J. Opt. Soc. Am.*, vol. 50, no. 9, pp. 838-844, Sep 1960. [Online]. Available: https://doi.org/10.1364/JOSA.50.000838

[38] J. MacArthur, "Design of the SEASAT-A radar altimeter," in *OCEANS '76, 1976, pp. 222-229*. [Online]. Available: https://doi.org/10.1109/OCEANS.1976.1154217

[39] M. S. Longuet-Higgins, "The statistical analysis of a random, moving surface," *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, vol. 249, no. 966, pp. 321-387, 1957. [Online]. Available: https://doi.org/10.1098/rsta.1957.0002

[40] ——, "The distribution of the sizes of images reflected in a random surface," *Mathematical Proceedings of the Cambridge Philosophical Society*, vol. 55, no. 1, p. 91-100, 1959. [Online]. Available: https://doi.org/10.1017/S0305004100033739

[41] R. G. Gardachov, "The probability density of the total curvature of a uniform random Gaussian sea surface in the specular points," *International Journal of Remote Sensing*, vol. 21, no. 15, pp. 2917-2926, 2000. [Online]. Available: https://doi.org/10.1080/01431160050121320

[42] T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, "A unified directional spectrum for long and short wind-driven waves," *Journal of Geophysical Research*, vol. 102, no. C7, pp. 15781-15796, 1997.

[43] Y. Quilfen and B. Chapron, "Ocean surface wave-current signatures from satellite altimeter measurements," *Geophysical Research Letters*, vol. 46, no. 1, pp. 253-261, 2019. [Online]. Available: https://doi.org/10.1029/2018GL081029

[44] F. Nouguier, B. Chapron, F. Collard, A. Mouche, N. Rascle, F. Ardhuin, and X. Wu, "Sea surface kinematics from near-nadir radar measurements," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 10, pp. 6169-6179, 2018. [Online]. Available: https://doi.org/10.1109/TGRS.2018.2833200

[45] C. Cox and W. Munk, "Measurements of the roughness of the sea surface from photographs of the Sun's glitter," *Journal of the optical Society of America*, vol. 44, no. 11, pp. 838-850, 1954. [Online]. Available: https://doi.org/10.1364/josa.44.000838

[46] W. Munk, "An inconvenient sea truth: spread, steepness and skewness of surface slopes," *Annual Review of Marine Science*, vol. 1, pp. 377-415, 2009. [Online]. Available: https://doi.org/10.1146/annurev.marine.010908.163940

[47] F. M. Bréon and N. Henriot, "Spaceborne observations of ocean glint reflectance and modeling of wave slope distributions," *Journal of Geophysical Research*, vol. 111, p. C0605, 2006. [Online]. Available: https://doi.org/10.1029/2005JC003343

[48] N. Tran and B. Chapron, "Combined wind vector and sea state impact on ocean nadir-viewing Ku- and C-band radar cross-sections," *Sensors*, vol. 6, no. 3, pp. 193-207, 2006. [Online]. Available: https://doi.org/10.3390/s6030193

[49] I. S. Gradshteyn and I. M. Ryzhik, *Table of integrals, series, and products*. Academic press, 2014.

[50] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions*. New York: Dover Publication, 1965.
<figure><img src="image_19.png" /></figure>

**Louis Marié** received the B.Eng. degree from Ecole Nationale Supérieure des Mines de Paris, France, in 1999, and the Ph.D. degree in physics (fluid dynamics) from the Université Denis Diderot, Paris, France, in 2003. He is 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 is a field oceanographer specialized in shelf seas dynamics. He has experience in non-linear physics, dynamo theory, turbulence theory and electromagnetism.

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

**Bertrand Chapron** received the B.Eng. degree from the Institut National Polytechnique de Grenoble, Grenoble, France, in 1984, and the Ph.D. degree in physics (fluid mechanics) from the University of Aix-Marseille II, Marseille, France, in 1988. He spent three years as a Post-Doctoral Research Associate with the NASA Goddard Space Flight Center, Greenbelt, MD, USA, and the Wallops Flight Facility, Wallops Island, VA, USA. 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.

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

**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.

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

**Doug Vandemark** received the bachelor's degree in physics from the Hope College, New York, NY, USA, in 1986, the master's degree in electrical engineering from the Microwave Remote Sensing Laboratory, University of Massachusetts at Amherst, Amherst, MA, USA, in 1988, and the Ph.D. degree in Earth science from the University of New Hampshire, Durham, NH, USA, in 2005. He was a Research Engineer and a Scientist with the Goddard Space Flight Center, National Aeronautics and Space Administration (NASA), Washington, DC, USA, from 1990 to 2005. He is active in studies tied to numerous NASA and European Space Agency Earth observing satellites and serves as a Principal Investigator on NASA's Ocean Surface Topography, Ocean Surface Salinity, Ocean Vector Wind, and Surface Water and Ocean Topography Science Teams. His research interests include the study of ocean and atmosphere boundary layer interactions and how to remotely measure near-surface ocean properties using radar and radiometer systems., Dr. Vandemark has served as a Strategic Advisory Team Member for the Northeast Regional Association of Coastal Ocean Observing Systems, the NASA-Jet Propulsion Laboratory (JPL) Physical Oceanography Data Archive Center, and the U.S. He has served as the Chair and a member of the NASA Earth Science Senior Reviews. He has served as an Associate Editor for the IEEE Geoscience and Remote Sensing Transactions.

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

**Fabrice Ardhuin** graduated from Ecole Polytechnique, Palaiseau, France, in 1997, and the Ph.D. degree in oceanography from the U.S. Naval Postgraduate School, Monterey, CA, USA, in 2001. He has been involved in ocean waves and related topics, from microseisms to remote sensing, with a strong focus on wave-current and wave-ice interactions, first at the French Navy Hydrographic and Oceanographic Service, Brest, France, then Institut Français de Recherche pour l'Exploitation de la Mer, Plouzané, France, and now CNRS, at the Laboratory for Ocean Physics and Satellite remote sensing. He led the SKIM proposal from 2016 to 2019 and is a member of the SWOT Science Team, the SEASTAR Mission Advisory Group, a co-investigator for the ODYSEA mission proposal, and the lead scientist for the ESA Sea State Climate Change Initiative project.