# Along-Track Resolution and Uncertainty of Altimeter-Derived Wave Height and Sea Level: Re-Defining the Significant Wave Height in Extreme Storms

Marine de Carlo, Fabrice Ardhuin

► To cite this version:

Marine de Carlo, Fabrice Ardhuin. Along-Track Resolution and Uncertainty of Altimeter-Derived Wave Height and Sea Level: Re-Defining the Significant Wave Height in Extreme Storms. Journal of Geophysical Research. Oceans, 2024, 129, 10.1029/2023jc020832. hal-04676630

HAL Id: hal-04676630

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

Submitted on 23 Aug 2024

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
**Key Points:**

* Retracking of altimeter waveforms yields fluctuations in wave height and sea level, correlated at the scale of the effective footprint

* A good approximation for the effective footprint diameter is the square root of the product of wave height and satellite altitude

* Estimating phenomenal wave heights precise to better than 2% requires filtering over wave groups, typically over a distance of 20–50 km

**Correspondence to:**

F. Ardhuin,
ardhuin@ifremer.fr

**Citation:**

De Carlo, M., & Ardhuin, F. (2024). Along-track resolution and uncertainty of altimeter-derived wave height and sea level: Re-defining the significant wave height in extreme storms. *Journal of Geophysical Research: Oceans, 129*, e2023JC020832. https://doi.org/10.1029/2023JC020832

Received 20 DEC 2023
Accepted 11 MAY 2024

**Author Contributions:**

**Conceptualization:** Fabrice Ardhuin

**Funding acquisition:** Fabrice Ardhuin

**Investigation:** Marine De Carlo

**Methodology:** Fabrice Ardhuin

**Software:** Marine De Carlo,
Fabrice Ardhuin

**Supervision:** Fabrice Ardhuin

**Validation:** Marine De Carlo

**Visualization:** Fabrice Ardhuin

**Writing – original draft:** Fabrice Ardhuin

**Writing – review & editing:**

Marine De Carlo, Fabrice Ardhuin

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

# Along-Track Resolution and Uncertainty of Altimeter-Derived Wave Height and Sea Level: Re-Defining the Significant Wave Height in Extreme Storms

Marine De Carlo¹ ID and Fabrice Ardhuin¹ ID

¹University Brest, CNRS, Ifremer, IRD, IUEM, Laboratoire d'Océanographie Physique et Spatiale (LOPS), Plouzané, France

**Abstract** Satellite altimeters are the most common source of wave measurement in phenomenal sea states, with significant wave heights exceeding 14 m. Unfortunately their data is still considered with skepticism, because there is usually no other data to verify the accuracy of the largest values. Here we investigate the self-consistency of the measurement, and their small scale variability, in order to define an estimate of satellite altimeter precision. Using numerical simulations of ocean surfaces and the processing involved in satellite retracking, we find that wave groups are responsible for a variance in estimated altimeter wave heights that is proportional to the square of the spectral peakedness parameter and the significant wave height. Additional variance induced by speckle noise is proportional to the wave height. The effect of wave groups generally dominates in the most severe storms. This variability requires a relatively large scale smoothing or filtering to yield accurate wave height estimates. For example, the largest ever reported 1 s average significant wave height from altimeters sampled by Jason-2 in the North Atlantic in 2011, at $\bar{H}_s = 20.1$ m, is now interpreted to correspond to a true wave height $H_s = 18.5 \pm 0.3$ m. The difference between 20.1 and 18.5 m is mostly due to wave group contributions to the raw measurement. We argue that wave group effects should not be included in the definition of the significant wave height, just like the maximum wave height differs from the significant wave height.

**Plain Language Summary** Over most of the past 30 years, satellite altimeters have been the only means to measure wave properties in the most severe ocean storms. How do we know that these data are trustworthy, and how can we define uncertainties? Here we show that as a satellite flies along its orbit, it reports wave height that fluctuate because of the random nature of the wavefield that can be organized in groups at the scale of a few kilometers. We are able to simulate the precision of the measurements, as a function of the wave height and the degree of organization of the wavefield, measured by a “spectral peakedness” parameter. This novel understanding can be used to define the precision of the measurements. For example, as far as we know, the largest reported value for a 1 s averaged satellite measurement of the significant wave height was $\bar{H}_s = 20.1$ m in a 2011 North Atlantic Storm, with no precision given. We can now re-interpret this data as evidence of a true significant wave height $H_s = 18.5 \pm 0.3$ m. The local fluctuations up to 20 m are caused by wave groups and should not be counted in the significant wave height.

## 1. Introduction

Satellite altimeters have been used over the past 30 years to measure sea level (Cazenave et al., 2018) and sea states (Ardhuin, Stopa, et al., 2019). These measurements are based on the estimated distances between a radar and the scattering elements at the sea surface, with a ‘local average’ related to the sea level and a ‘local standard deviation’ related to the significant wave height. This separation was understood well enough for most applications, but new instruments able to resolve shorter and shorter scales make it more important to clarify how the multi-scale ocean surface elevations and velocities contribute to the parameters estimated from altimeter data. In particular we explore the link between the underlying significant wave height, the ‘local standard deviation’ of the surface elevation and the altimeter measurements.

Our goal in this paper is to build a model for the small scale fluctuations in wave height estimates, given below by Equations 22 and 23. We apply this model to propose an uncertainty for altimeter measurements of large significant wave heights ($H_s > 8$ m) for which too few validation data exist (Dodet et al., 2020). Understanding these fluctuations is also relevant in the context of recent efforts to improve instruments and data processing techniques
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two groups of 50 consecutive CFOSAT/SWIM nadir waveforms spanning 11 s each (i.e., a distance of 75 km along-track), individual waveforms are color-coded with the estimated wave height $\hat{H}_s$. Both groups were acquired along the same descending orbit in the North Atlantic on 14 February 2020 around 9:10 UTC, on the edges of storm Dennis (a) around 61.5° N and (b) around 44.5° N. These are the echo_L1A variable, already corrected for the antenna pattern, and normalized by the estimated Level 2 Normalized Radar Cross Section. The white line represents the average waveform. The horizontal dashed lines represent the 98% confidence interval expected for random fluctuations due to speckle, assuming that 264 independent radar pulses are averaged for each range gate.</figcaption></figure>

to provide the highest possible resolution, in particular for coastal areas (Passaro et al., 2021; Vignudelli et al., 2018). The variability of sea level and $H_s$ estimates is generally well understood at scales larger than 30 km, where geostrophic currents and their effect on wave heights dominate (Ardhuin et al., 2017; Morrow & Le Traon, 2012). Extending this understanding toward high resolution requires a detailed analysis of the measurement system. The present paper extends the previous analysis by De Carlo et al. (2023), hereinafter DC23, with a particular emphasis on the correlation properties of the measured data.

The fundamental measurement of an altimeter is the power received from sea echoes as a function of delay time $t$, known as the waveform (Brown, 1977). In practice the time separation is not perfect and some blurring in time is caused by the finite frequency bandwidth of the radar, so that the measured waveform is a convolution of the true waveform and the instrument point target response (PTR). Example of waveforms from the China-France Ocean Satellite (CFOSAT) are shown in Figure 1. They correspond to the same average significant wave height, around 9.3 m, but very different sea states, as quantified by their spectral peakedness $Q_{kk}$, a parameter defined below, with a young wind sea on the left and a mature swell on the right. The method used to estimate sea level and wave height uses the fit of a theoretical waveform shape to the measured waveform. That theoretical shape is also, but not only, a function of sea level and wave height.

A well known source of deviations from the theoretical shape is the purely instrumental effect of 'speckle noise' which comes from Rayleigh fading: this noise is present when the propagation paths between the radar and individual scattering elements at the ocean surface have lengths that spread over a range much larger than the radar wavelength (Quartly et al., 2001). In both panels of Figure 1, speckle noise explains the fluctuations for range gates indices larger than 110: on average, for half of the range gates one waveform (out of 50) should exceed the upper dashed line and one should fall below the lower dashed line.

Another well known source of deviations from the theoretical waveform is the effect of ocean backscatter variability within the radar footprint, which is very important for wind speeds under 3 m/s (Dibarbourie et al., 2014), and in the presence of sea ice (Tourain et al., 2021). Both speckle noise and backscatter variability have been cited as the possible source of noise in the estimation of wave height and sea level, and Sandwell and Smith (2005) have explained the resulting correlation of these two retracked parameters. The observed correlation can be used empirically to reduce the noise in sea level estimates (Quartly et al., 2019; Zaron & DeCarvalho, 2016).

Here we find correlations of along-track variations of parameters, shown in Figure 2, that are consistent with previous studies, but with much larger relative fluctuations of both wave height and sea level estimates. We will
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Estimates of the (a) epoch and (b) wave heights for the 50 waveforms of Figure 1b, using two different cost functions: “LS” is a least squares 2-parameter fit to the theoretical waveform given by Equation A17 with $a = b = \xi = 0$, $\gamma = 1$ and “ML” is a maximum likelihood 2-parameter fit to the same theoretical waveform. The “native” data is shown for reference and is the operational method as described in Tourain et al. (2021). The good agreement of the “ML” retracking with the “native” data requires to ignore the first range gates using $k_{min} \approx 80$, or adapting $k_{min}$ such that $S(k_{min}) > r_{min} \max(S)$, with $r_{min} \approx 0.06$.</figcaption></figure>

show that this magnitude is specific to sea states with large wavelengths and narrow spectra, an effect that can uniquely explain why the waveforms in the two panels of Figure 1 differ so much for range gate indices 85 to 100. In fact, a much less researched source of deviations between measured waveforms and parametric models, is the non-uniform statistical distribution of the elevation of sea surface scatterers at the scale of the “instrument footprint” (to be precisely defined below). The first study of that effect was published by DC23, with the main results summarized at the beginning of Section 2. The waveforms in Figure 1a are consistent with the assumption of uniform wave heights across the footprint that was used to derive parametric theoretical waveforms, whereas in Figure 1b the waveforms are more different from the theoretical shape, which is typical of non-uniform conditions. Fitting these different waveforms with the theoretical shape gives wild variations of the wave height estimate $\hat{H}_s$, shown in Figure 2a, that may not be realistic, and even wilder variations of the sea level anomaly in Figure 2b, with differences up to 1.8 m for measurements only 19 km apart. This sharp gradient in that region of the ocean (middle of the North Atlantic Ocean) is clearly not realistic and can be a spurious effect of the violated elevation uniformity assumption.

In Section 2 we provide a basis for the understanding of perturbation of altimeter measurement associated to wave groups. This extends the approach of DC23 to more realistic representation of altimeter waveforms. Section 3 builds the uncertainty model, starting from the uncertainty of individual measurements and, after defining the relevant along-track correlation scales, defining an uncertainty for along-track averages. These effects are illustrated using simulated waveforms for the sea state conditions with strong wave groups corresponding to Figure 1b. The contribution of wave groups is a simple function of $H_s$, $Q_{kk}$, and instrument parameters. The same method can be applied to sea level fluctuations. This theoretical model is verified in Section 4 using simulated waveforms corresponding to a large range of combinations of $H_s$ and $Q_{kk}$ that may be found in the global ocean. This model is then applied to the estimation of $H_s$ uncertainty from a sequence of $\hat{H}_s$ estimates obtained from individual waveforms. Summary and conclusions follow in Section 5.

## 2. Waveforms and Their Retracking Over Wave Height Gradients

Here we extend the work of DC23 who neglected the Earth sphericity, assumed a broad antenna pattern and neglected effects of the PTR and of speckle noise. We relax these assumptions and investigate the influence of the choice of the cost function. We provide analytical derivations of the forward model (a generalized parametric waveform) with details in Appendix A. We emphasize that this parametric waveform was not designed for retracking, but rather to guide the interpretation of existing data sets that are the results of retracking with the usual Brown waveform. Indeed the inverse modeling (the retracking) could be done analytically in the case of DC23. For the generalized waveform the analytical retracking may not be feasible but for cost functions based on least
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Geometry of the measurement in the simplified case of a flat mean sea surface, and in the presence of wave groups giving a significant wave height $H_s = 2a$. (b) Corresponding waveforms with the x-axis showing the range distance from the satellite. In (a) the wave height is exaggerated and the satellite height is reduced by the same factor giving the correct Chelton radius $\rho_C \approx \sqrt{2hH_s}$ defined as the distance from nadir where the mean sea surface is at a distance $h + H_s$ from the satellite.</figcaption></figure>

squares the retracked wave height and sea level can still be expressed as functions of anomalies of the wave height field.

## 2.1. Footprints

The “radar footprint” is the region of the ocean that produces backscatter from a radar single pulse, and depends on the antenna aperture, satellite altitude $h$ and Earth radius $R_E$. In the case of the nadir beam on SWIM, this is typically a disc of radius 9.5 km centered on the nadir, where we have defined the footprint boundary as the location where the backscattered power drops to half the peak power at nadir. There are notable exceptions with strongly reflecting surfaces at high elevations above sea level (land or icebergs) that may corrupt the measurements even if they are further than 9.5 km from nadir. When we exclude these exceptional cases, the measurements are mostly sensitive to an area much smaller than the radar footprint. For those altimeters that only measure power as a function of time delay, Chelton et al. (1989) have argued that the estimates $\hat{H}_s$ for wave height and $\hat{z}_e$ for sea level are associated to the true physical values of the significant wave height $H_s$ and sea level $z_e$ within an “oceanographic footprint” which they defined to be a disc of radius

$$ \rho_C = \sqrt{2h(H_s + \delta_R)/(1 + h/R_E)}, \quad (1) $$

where the range resolution $\delta_R = c/(2 B)$ is defined by the radar bandwidth $B$ and the speed of light $c$. So far, all altimeters that use a Ku-band frequency have used $B = 320$ MHz giving $\delta_R = 0.47$ m so that the minimum radius $\rho_C$, corresponding to the lowest sea states, is of the order of 1 km. For a very large sea state with $H_s = 9$ m and the relatively low orbit height of 519 km of CFOSAT, one gets $\rho_C = 3.1$ km.

However, that estimate turns out to be very conservative. Data from the SWIM instrument on CFOSAT occasionally shows meaningful variations in $\hat{H}_s$ between consecutive measurements separated by only 1.7 km, for example, over coral reefs (Alice Dalphinet, personal communication). In Figure 2 these variations of $\hat{H}_s$ are as large as 3 m and may be caused by wave groups. So what is the effective diameter of the footprint of a satellite altimeter?

## 2.2. The Wave Group Effect and DC23 Results

Non-uniform wave heights occur even in homogeneous sea states represented by a single wave spectrum, due to the interference of wave trains with different frequencies and directions. This interference produces series of high waves known as wave groups (Arhan & Ezraty, 1978). We will therefore call this particular non-homogeneity the “wave group effect”. It is present for all sea states, albeit with different magnitudes. The most simple form of wave groups is shown in Figure 3 with the sum of two monochromatic wave trains, of wavenumbers $k_1$ and $k_2$ forming a beating pattern. It is obvious that waveforms obtained at times $t_1$ and $t_2$ are different: at $t_1$ the first echoes correspond to the distance $h - a$ whereas at $t_2$ the first echoes arrive later and correspond almost to $h$. As a result, the corresponding wave height estimates $\hat{H}_s$ differ by a factor 10, even though the sea state is “spatially
homogeneous”, in the sense that the corresponding wave spectrum and associated parameters, including the underlying significant wave height $H_s$, are constant.

It is common to study the properties of wave groups by introducing the wave envelope (Rice, 1944; Tayfun & Lo, 1989), which defines a local wave amplitude $\eta$ from the surface elevation $\zeta$, as represented in Figure 3. In the bi-chromatic case of Figure 3, the envelope varies at scales given by the wavenumber $K = |k_1 - k_2|$. A realistic sea state is the sum of many monochromatic components with a range of wavenumber vectors $k$. An important result is that the envelope contains all the spatial scales larger than the scale of dominant wave group, that is, with all the wavenumbers $K = k \pm k'$, including $K = 0$. Namely, whereas the elevation associated with a given sea state (outside of long swells or very severe storms) may not contain any wavelengths longer than say 400 m, the envelope of that same sea state does vary at all scales, including tens of kilometers.

From the envelope $\eta$, obtained from the analytical form of the surface elevation (see DC23 for details), we defined a local wave height

$$H_l = 4\sqrt{2/\pi} \times \eta, \qquad (2)$$

With a scaling a little different from Janssen (2014) so that the large-scale average of $H_l$ is the usual significant wave height defined from the average of the surface elevation variance $\langle H_l \rangle = H_s = 4\sqrt{\langle \zeta^2 \rangle}$. More specifically, the Power Spectral Density (PSD) of the envelope, and thus the PSD of the local wave height $H_l$, is proportional to the convolution of the double-sided wave spectrum $E(k_x, k_y)$ by itself (Tayfun & Lo, 1989). In particular, for a Gaussian wave spectrum, the envelope spectrum is also a Gaussian, but centered on $K = 0$, as detailed in DC23.

When concerned about fluctuations of the wave height $H_l$ filtered at scales much larger than the dominant wave groups, one can approximate the PSD of $H_l$ as a constant, and obtain the variance of $H_l$ as the value of the $H_l$ PSD at $K = 0$ times the spectral integral of the filter response $\Delta_k^2$. DC23 showed that this gives the following variance associated to wave groups

$$\text{var}_{\text{wg}}(H_l) = (4 - \pi)H_s^2 Q_{kk}^2 \Delta_k^2, \qquad (3)$$

Where they have defined the spectral peakedness as

$$Q_{kk}^2 = \frac{\iint_{\mathbb{R}^2} E^2(k_x, k_y) dk_x dk_y}{\left(\iint_{\mathbb{R}^2} E(k_x, k_y) dk_x dk_y\right)^2}. \qquad (4)$$

The link between these properties and altimeter data was made explicit by DC23 who showed that the output of the simplest 2-parameter (wave height and sea level) altimeter retracker can be expressed analytically as a spatial filter of the $H_l$ field.

More precisely they have shown that the magnitude of $\hat{H}_s$ fluctuations is consistent with smoothing the $H_l$ field with a two-dimensional Gaussian filter of parameter $\sigma \approx \rho_c/4.5$. Furthermore, in their Annex A, they introduced an ideally perturbed waveform, and analytically calculated the impact of the perturbation on a 2-parameter least-square fit with a non-perturbed waveform. In that work the authors neglected the Earth sphericity, assumed a broad antenna pattern and neglected effects of the PTR and of speckle noise. Moreover, they only considered a very simple altimeter retracker.

Here we extend the work of DC23 by relaxing these assumptions and investigating the influence of the choice of the cost function for retracking. To do so, in Section 2.3 we define the different retracking cost functions we will be using and in Section 2.4 we extend the DC23 perturbation theory to more realistic waveforms. In Section 2.5, this will lead to new functions representing the impact of these perturbation on retracked parameters. And the generalization of the spatial filtering of $H_l$ is shown in Section 2.6.

## 2.3. Definitions of Retracking Cost Functions

In the following, we will fit the measured or simulated waveform ($y_k$) with a parametric form ($s_k$), for all range gate indices $k$ between $k_{min}$ and $k_{max}$. One possibility is to use a least squares cost function
equivalent footprint:
radius ρ<sub>C</sub>/2, area A<sub>eq</sub>=πρ<sub>C</sub><sup>2</sup>/4, wave height: H<sub>s</sub>

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

perturbation: area A₀ distance ρ₀ from nadir,
wave height: (1+Δ)Hₛ

<caption>**Figure 4.** Schematic of idealized sea surface anomalies located at a given distance ρ₀. This distribution is obviously impossible to obtain with real waves: a real wavefield will have a smoothly varying distribution of Hᵢ as a function of ρ.</caption>

$$
C_{LS} = \sum_{k=k_{\min}}^{k_{\max}} w(k)(y_k - s_k)^2 \quad (5)
$$

With the default weights w set to 1 for all k. This is used in the very common so-called “3 or 4-parameters Maximum Likelihood Estimator” (MLE3/ MLE4). In the present context a better name for these would be LS3 and LS4, with the 3 fitted parameters being Ĥ<sub>s</sub>, ˆz<sub>e</sub>, and the normalized radar cross section ˆσ<sub>0</sub>, and the fourth parameter is generally the antenna mispointing angle (Schlembach et al., 2020).

We may also use a maximum likelihood (ML) fit, first introduced by Rodriguez (1988) and particularly developed for the ERS-1 altimeter by Challenor and Srokosz (1989), and later used by Gómez-Enri et al. (2007). ML is the optimal method for a uniform sea state with fluctuations in the waveforms dominated by speckle noise. In the limit of a large number of looks it takes the following form,

$$
C_{\text{ML}} = \sum_{k=k_{\text{min}}}^{k_{\text{max}}} \frac{y_k + \varepsilon}{s_k + \varepsilon} - \log\left(\frac{y_k + \varepsilon}{s_k + \varepsilon}\right), \quad (6)
$$

Where we have introduced ε = 10⁻⁵ to reduce the influence of numerical errors. We also define a parameter r_min
with values between 0 and 1, to exclude the first range gates that have a power lower r_min times the maximum
value of the waveform, up to the index k_min namely,

$$
S(k_{\text{min}}) > r_{\text{min}} \max(y_k). \tag{7}
$$

The ML-type cost function is used in the “adaptative” method used to produce the “native” CFOSAT data (Tourain et al., 2021), but it is not clear which are the actual range gates used in practice. Although we initially used r_min = 0, we found a generally good agreement with the native CFOSAT data when using r_min = 0.06. In the example on Figures 5b and 5d, using r_min = 0.06 corresponds to fitting only the part of the waveform that is above the horizontal dashed line. More details on the sensitivity of results to the value of r_min are given in Appendix C. An intermediate cost function is used in the WHALES retracker (Passaro & Algorithm Development Team, 2021; Schlembach et al., 2020). It is a weighted least squares with much larger weights w(k) for the early part of the waveform, defined by the inverse of the standard deviation of waveform residuals caused by speckle noise for the same wave height.

Taking the waveforms of Figure 1b as an example, LS-based retracking has less variability than the ML-based result, their mean values differ by 27 cm, and both retrackers give a strong correlation between epoch and Hs anomalies, shown in Figure 2b. Using simulated sea surfaces and altimeter waveforms (see details in Appendix B for the simulation method), we will show in Section 4 that this example is actually representative of large sea states with narrow wave spectra. In these cases speckle noise is a less important source of waveform deviations than the wave group effect, and ML-based retrackers are not optimal.

2.4. Wave Height Gradient Effect on Waveforms: Beyond DC23

Following DC23, we start by deriving an analytical perturbed waveform in the presence of an unrealistic localized anomaly in surface elevation. This is detailed in Appendix A. Our anomaly consists of a change in significant wave height Hs, defined as 4 times the standard deviation of the surface elevation, from a background value Hs to an anomalous value Hs(1 + Δ) over an area A0 centered at the distance from nadir ρ0. Both the normal and anomalous sea levels are taken to be Gaussian distributed. The distance from nadir ρ0 correspond to a distance h + R0 from the altimeter at the mean sea level. The wave height anomaly can be localized or distributed over a ring, as shown in Figure 4. More realistically, the anomaly is the ring-average of the true local wave heights Hl for the distance ρ0. After dealing with the kind of anomalies represented in Figure 4, we will consider the superposition of all the distributions at all the distances from nadir. Those anomalies can be caused by wave groups but
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Example of theoretical perturbed waveforms based on Equation A17 for $H_s = 10$ m, $\tau = 0$, $\xi = 0$, $N_t = 0.001$, $a = 0.3$ and $b = 0$ and its comparison to unperturbed waveforms with $a = 0$. The x axis is the delay time normalized by $\sigma_s = 2H_s/c$. Left panels (a, b) are obtained in the limit of very broad antenna pattern, and (c, d) correspond to the real SWIM antenna pattern with $\theta_{3dB} = 1.6^\circ$. The bottom panels (b, d) correspond to the same waveforms but plotted with a logarithmic y axis, and the shaded area is ignored in the ML fit using $r_{min} = 0.06$. All waveforms use the PTR given by Equation A16.</figcaption></figure>

also by many other processes (wave breaking over a coral reef, dissipation over a mud bank, wave-current interactions ...). We define the equivalent footprint area $A_e = \pi h H_s (1 + h/R_E)/2$, and find that the local wave height anomaly adds a perturbation to the usual waveforms, as given by Equations A10 and A17. The two adimensional parameters that define magnitude and location of the perturbation are.

$$a = \Delta \frac{2A_0}{\pi h H_s / (1 + h / R_E)} = \Delta \frac{A_0}{A_{eq}} \qquad (8)$$

$$b = \frac{R_0}{H_s} = \frac{\rho_0^2}{2hH_s(1 + h/R_E)} \approx (\rho_0 / \rho_C)^2. \qquad (9)$$

One example of this theoretical perturbed waveforms is shown in Figure 5 for $a = 0.3$ and $b = 0$, compared to the unperturbed Brown waveform plotted with a dashed line, and different attempts at fitting it with a Brown waveform. Taking into account the full complexity of the waveform (right panels) does not change the qualitative impact on the simplest possible waveforms (left panels) used in DC23. In our example perturbation, $a > 0$ means that the wave heights are locally higher, which tends to shift some of the echoes to shorter and larger ranges: the black curve is higher than the dashed curve for $|t - \tau| \approx 2\sigma_s$.

We have plotted the waveforms using both linear (top panels) and logarithmic (bottom panels) coordinates to illustrate the fact that the Maximum Likelihood cost function uses ratios instead of differences and gives a better fit in logarithmic coordinates.

When $b = 0$, our perturbed waveform is identical to the waveform for a uniform but non-Gaussian skewed surface elevation distribution, with skewness parameter $\lambda_{3,0,0} = 6a$ (Hayne, 1980; Srokosz, 1986). Although $a = 0.3$ is a fairly large but not impossible wave group effect, it would correspond to an impossibly large $\lambda = 1.8$. Hence the
<caption>**Table 1**</caption>

Fitted Parameters (in Meters), $\hat{H}_s$ and Epoch Converted to Distance, for the ($a = 0.3$, $b = 0$) Waveform Shown in Figures 5a and 5b, and a Few Other Examples, Using Three Different Cost Functions: The Least Square of Equation 5 and Two Versions of the Maximum Likelihood of Equation 6, With $R_{min} = 0$ or $R_{min} = 0.06$

<table><thead><tr><th rowspan="2">(a, b)</th><th colspan="2">LS fit</th><th colspan="2">Equation 10</th><th colspan="2">ML fit</th><th colspan="2">ML, r<sub>min</sub> = 0.06</th></tr><tr><th>\hat{H}<sub>s</sub></th><th>\hat{z}<sub>e</sub></th><th>\hat{H}<sub>s</sub></th><th>\hat{z}<sub>e</sub></th><th>\hat{H}<sub>s</sub></th><th>\hat{z}<sub>e</sub></th><th>\hat{H}<sub>s</sub></th><th>\hat{z}<sub>e</sub></th></tr></thead><tbody><tr><td>(0.00, 0.00)</td><td>10.0</td><td>0.00</td><td>10.0</td><td>0.0</td><td>10.0</td><td>0.00</td><td>10.0</td><td>0.00</td></tr><tr><td>(0.30, 0.00)</td><td>9.5</td><td>0.42</td><td>10.0</td><td>0.38</td><td>13.5</td><td>0.29</td><td>14.9</td><td>0.79</td></tr><tr><td>(-0.03, 0.00)</td><td>10.0</td><td>-0.04</td><td>10.0</td><td>-0.04</td><td>8.9</td><td>-0.25</td><td>9.6</td><td>-0.08</td></tr><tr><td>(0.30, 0.25)</td><td>12.8</td><td>0.22</td><td>12.9</td><td>0.14</td><td>10.1</td><td>-0.33</td><td>11.7</td><td>0.05</td></tr><tr><td>(-0.30, 0.25)</td><td>7.1</td><td>-0.12</td><td>7.1</td><td>-0.14</td><td>9.4</td><td>0.23</td><td>8.10</td><td>-0.05</td></tr></tbody></table>

wave group effect can be locally much larger than the skewness effect. We note that the spurious perturbations on the sea level estimate $\hat{z}_e$ is a 'tracker bias' (not a true physical effect) since the model waveforms correspond to a zero sea level and we do not take into account non-uniform scattering along the surface (the electro-magnetic bias). However, when averaging waveforms along the altimeter track, skewness persists but the wave group effect should vanish because $a$, the amplitude of wave group effects, is symmetrically distributed around zero. Consequences for retracking with generalized waveforms are discussed in Appendix D.

## 2.5. Influence of Idealized Wave Field Anomalies on Retracked Parameters

Our wavefield anomaly of amplitude $a$ and location $b$ gave us the perturbation to the waveform that in turn produces a perturbation of the retracked parameters: the estimates of wave height $\hat{H}_s$ and sea level $\hat{z}_e$. Using the analytical form of the perturbed waveform (obtained in the limit of a broad antenna pattern, i.e. using $v = 0$ in Equation A10, and ignoring the PTR), DC23 have computed the cost function $C_{LS}$ analytically, replacing the discrete sum by an integral over all ranges from minus infinity to plus infinity. Taking its derivative with respect to $\hat{H}_s$ and $\hat{z}_e$, they found that the cost function is minimum for these values of the retracked parameters.

$$ \hat{H}_s = H_s + \frac{aH_s}{2} J_H(b), \qquad (10) $$

$$ \hat{z}_e = -c\hat{\tau}/2 = -\frac{aH_s}{16} J_z(b), \qquad (11) $$

With $b = \rho_0^2/\rho_C^2 = \rho_0^2/(2hH_s)$ and.

$$ J_H(b) = 2b(6 - 16b^2)e^{-4b^2}, \qquad (12) $$

$$ J_z(b) = (2 - 16b^2)e^{-4b^2}. \qquad (13) $$

The perturbation caused by the wave anomaly on the retracked parameter is proportional to $aH_s$ and, a function of the off-nadir distance $\rho_0$ which we normalize as $b$. Further interpretation of $J_H$ is given in Section 2.3. In simple terms, the large values of $J_H(b)$ for $b$ up to 0.30 are the main driver of the along-track correlation scale, as will be explained below. The analytical perturbations in Equations 10 and 11 are typically accurate within 10% for LS retracking and $a < 0.2$, with some examples given in Table 1, showing that it is in fact fairly robust up to $a = 0.3$.

For more realistic waveforms or different cost functions such as ML or WHALES, there are no simple analytical solution. One can still analyze the perturbations of the retracked parameters and interpret the results by computing the following functions

$$ \hat{J}_H(b) = 2(\hat{H}_s - H_s)/a, \qquad (14) $$

$$ \hat{J}_z(b) = 16(c\hat{\tau}/2)/(aH_s), \qquad (15) $$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Influence of a local wave height anomalies as given by Equations 8 and 9 when retracking with a waveform that uses $a = b = 0$, for (a) wave height (b) epoch as a function of the perturbation distance from nadir defined by the parameter (b) All results were obtained for $a = 0.1$, except for the last curve with $a = -0.02$. The analytic expression are given by Equations 12 and 13 and are independent of (a) The numerical evaluations are given by Equations 14 and 15 and were obtained for different waveforms, with a broad ($v(t) = 0$ in Equation A10) or realistic beam, with or without PTR, with either LS or ML fitting. For realistic beams, the waveform power is also fitted, with very little influence on the adjustment of $\hat{H}_s$ and $\hat{t}$. For the WHALES retracker the weights in the least-square cost function are the weights used for Jason-2 retracking in the SeaState CCI-V3 dataset, for a wave height of 10 m.</figcaption></figure>

with results shown in Figure 6. These results demonstrate that the analytical derivation in DC23 does not exactly correspond to realistic waveforms, but it is qualitatively correct. The numerical estimates of the wave height perturbation $\hat{J}_H$ and sea level perturbation $\hat{J}_z$ were obtained for $a = 0.1$, progressively relaxing the different simplifying assumptions on the waveform: considering the instrument PTR, using a finite radar beam width. We also tested different cost function options: Least Squares (LS), and Maximum Likelihood (ML) with $r_{min} = 0$ and $r_{min} = 0.06$, and the MMSE of the WHALES retracker using the actual weights used for retracking Jason-2 in the SeaState CCI dataset. We note that relaxing the assumptions on the PTR has no visible effect when using LS fitting, and using the real radar beam width $\theta_{3dB} = 1.6^\circ$ instead of $v = 0$ also has a limited impact, especially for significant wave heights lower than the 10 m used here (not shown).

The $\hat{J}_H$ and $\hat{J}_z$ functions obtained with ML are very different from those obtained with LS: they are both maximum and larger for perturbations near nadir ($b = 0$), which explains the stronger correlation between epoch and wave height anomalies when using ML fitting, as shown in Figure 2c. As discussed above, the ML cost function introduces a very strong sensitivity to the early part of the waveform, and hence to near-nadir perturbations
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Surface elevation $ζ(x, y)$ and (b) local wave height $H_l$ around (x = 12, y = 43) for the sea state used in Figure 9. (c) Shows the values of $J_H(b)$ with $b = ρ^2/ρ_c^2$ and $ρ$ the distance from nadir. (d) $J_H$ multiplied by the local wave height anomaly, (e) sums of anomalies for each distance-to-nadir $b$, (f) waveform simulated from $ζ(x, y)$ and fits with different waveforms.</figcaption></figure>

($b < 0.15$): the estimated $\hat{H}_s$ can be corrupted by a very small area with very large waves. We have thus introduced the $r_{min}$ parameter as defined above in Section 2.3. Even with this adjustment, the ML-estimated $\hat{H}_s$ is a non-linear function of the perturbation amplitude $a$, as shown with the dotted line in Figure 6, obtained with a negative wave height anomaly. Finally the WHALES retracker gives results close to ML retracking but are linear: the $\hat{J}_H$ and $\hat{J}_z$ are independent of the amplitude $a$ (not shown).

## 2.6. Generalization to Any Wave Field and $H_l$ Pattern

The analysis of localized wave height anomalies generalizes to any combination of anomalies when using LS fitting, because the retracked values are linear combinations of the perturbations for each anomaly. DC23 demonstrated that a good estimate of the retracked values $\hat{H}_s$ can be obtained directly by filtering the map of local wave heights $H_l$ using the functions $J_H$, without performing any retracking.

In the rest of this section, we have taken the most simple waveforms, as done in DC23, generated from the same surface used in that paper and in the next section. Figure 7 shows details of the surface and corresponding waveform simulations with a nadir position at (x = 11.7, y = 43.2). This is the location where LS retracking gives the highest value of $\hat{H}_s$. From the surface elevation in Figure 7a to the waveforms in Figure 7f, the altimeter processing can be approximated with the following steps. First we may ignore the phases and only consider local wave heights $H_l$ shown in (b), then we filter using the $J_H(b)$ filter in (c) to produce amplified anomalies in (d) that can be averaged for each normalized radius $b$ into a value $H_b(b)$, as shown in (e), before summing the contributions for all radii to provide the local estimate 10.8 m. The local retracked value is $\hat{H}_s = 12.2$ m and both are significantly larger than the true wave height $H_s = 9.3$ m. This large local value is explained by the positive $H_l$ anomalies around $b = 0.25$ ($ρ = ρ_c/2$) where $J_H$ is positive and maximum, and the negative anomalies around $b = 0.9$ where $J_H$ is negative. The sum can also be done directly on all pixels of panel (d), in which case the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Scatter plots of retracked (a) wave height and (b) sea level for a 35 km all altimeter nadir position every 346 m in both x and y directions, giving 11,130 waveforms (without any noise added), compared to estimates using Equations 10 and 11. Red lines are best fit to the data, and the pixel circled in pink correspond to the case in Figure 7. (a) and (b) are obtained from the surface elevation $\zeta(x, y)$, (c) and (d) include a phase-average of two realizations $\zeta(x, y)$ and $-\zeta(x, y)$. Averaging over more realizations with different phases does not reduce further the scatter.</figcaption></figure>

equivalent perturbation amplitude is $a(x, y) = [H_l(x, y) - H_s]dx dy/(H_s A_{eq})$, and the contribution to $\hat{H}_s$ of each pixel is $J_H(b) \times a(x, y) \times H_s/2$, as given by Equation 10.

We note that if we multiply the surface elevation by -1, the crests become troughs and vice versa, leading to a slightly different waveform shown with the dashed line in Figure 7f, and a slightly different retracked value $\hat{H}_s^{-} = 11.8$ m, even though that surface has the exact same local wave heights $H_l$. It thus appears that the $J_H$ filter can give an interesting approximation of the altimeter result, but it cannot be exact, due to phase effects that it does not represent.

Testing further this $J_H$ filter idea, and the equivalent $J_z$ filter for the sea level, gives results shown in Figure 8, now looking at all 11,000 waveforms obtained from the same sea surface with nadir positions at all values of x and y spanning 35 km in each dimension. The right-most pixel of (a) corresponds to the case detailed in Figure 7. Our filter theory does not reproduce all the details of the variability in $\hat{H}_s$ and $\hat{z}_e$ estimates, but it explains 80%-90% of the variance. Here again, we have verified that changing the sign of the surface elevation gives a different
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a, b) Maps of sea level estimated from waveform retracking using ML or LS cost functions (c, d) maps of spatially filtered wave heights, using a two-dimensional Gaussian filter with a parameter $\sigma_l = 175$ m or $\sigma = \rho_c/4.5 = 708$ m (e, f) maps of retracked $\hat{H}_s$ corresponding to (a, b). (g) Is the retracking of the same waveforms but with $r_{min} = 0.06$ in ML cost function, and (h) is obtained with the LS cost function retracking of waveforms without speckle noise. The dashed circle around the filtered peak at $x = 19$ km, $y = 20$ km, has a radius $\rho_c/2 = 1.6$ km.</figcaption></figure>

estimation $\hat{H}_s^-$ and $\hat{\tau}_s^-$. Interestingly, the theoretical value is very close to the average of $\hat{H}_s$ and $\hat{H}_s^-$, as illustrated in Figures 8c and 8d.

# 3. A Model for Small Scale $\hat{H}_s$ Fluctuations

## 3.1. Retracking of Realistic Waveforms

In all the simulations discussed in this paper, the perturbation of the estimated epoch that is proportional to $H_s$ is completely spurious: it is a tracker error. In reality it will be combined with a true millimeter-scale sea level variation that is expected to scale like $H_s^2$ (Ardhuin et al., 2004; Longuet-Higgins & Stewart, 1962). In contrast, a large part of the variability of the wave height is real and may be of geophysical interest. In the case of wave groups, true perturbations of the local wave height $H_l$ travel at the group speed and do not persist for more than a few minutes. For other sources of gradients in wave heights, $H_l$ variability is related to a spatial gradient of $H_s$, and may persist longer and may be visible from one satellite pass to the next in the case of gradients caused by refraction over bottom topography or dissipation over mud banks. Whatever the source of the gradient in $H_l$, what is the smallest scale that can actually be resolved, and how well are the true patterns visible in altimeter data? To answer this question we use numerical simulations as described in Appendix B, starting from a directional wave spectra. In this section we use the same spectrum derived from CFOSAT L2S data that was already used in DC23.

Different retracking results are shown in Figure 9 and compared to smoothed local wave heights in panel (c) and (d), which are obtained directly from the surface elevation without any retracking. The smoothed local wave height $H_{l,\sigma_l}$, are obtained from $H_l$, given by Equation 2, with a two-dimensional Gaussian smoothing with a parameter $\sigma_l$.
The first striking result is that the retrackers give sea level estimates in panels (a) and (b) with significant variability (of the order of 40 cm for the ML fit), whereas the true sea level is actually flat in the simulation. That variability generally follow the large scales of the envelope in Figure 9d, and miss smaller details present in 9c.

For the wave heights, all retracking options used here give results that are visually clearly different from a simple Gaussian filter applied to the map of local wave heights $H_l$ in (c) or (d). The ML-retracked $\hat{H}_s$ in (e) is the most similar to the large-scale filtered local height $H_l$ in (d), with a maximum near $x = 19$ km and $y = 20$ km that has a similar shape, but this is not the case for other localized maxima in (d) that have ring shapes in (e). These ring shapes are much more present with LS retracking due to the shape of the $J_H$ function that is maximum for off-nadir perturbations, as further discussed in Appendix B. The radius of these rings is clearly related to the Chelton radius $\rho_c$, given by Equation 1, with a ring radius $\rho_c/2$ for LS retracking, corresponding to $b = 0.24$, and a smaller radius for ML retracking. Panel (g) was also obtained with ML but with $r_{min} = 0$, which gives different patterns. For example, the maximum at $x = 19$ km, $y = 20$ km gives a pattern reminiscent of a Mickey Mouse face with ears much more prominent than the nose. The fact that the “ears” are more prominent comes from the presence of two higher but much more narrow peaks in $H_l$, above 16 m. These higher values are not visible in Figure 9c because the color bar is saturated. ML fitting with $r_{min} = 0$ therefore gives a very particular distortion of the true map of wave heights that strongly emphasizes very high peaks even if they are very narrow. Finally we have also included in (h) one example with LS retracking of waveforms that do not include speckle, with the ring shapes now appearing more clearly than in (f). With ML the speckle has no visible impact for this sea state (not shown).

## 3.2. Linking Standard Deviation of $\hat{H}_s$ to Wave Spectral Shape

Because satellite altimetry is a technique more recent than in situ buoy measurement, the uncertainty of satellite data has generally been estimated based on buoy data (Abdalla et al., 2011; Dodet et al., 2020). These analyses have struggled to account for the fact that the two measurements cannot represent the same space-time coverage of the wavefield. Section 2.6 has now clarified that, for least squares cost functions, the spatial coverage of altimeters can be interpreted as a $\hat{J}_H$-filtered map of the local wave heights $H_l$, with the $\hat{J}_H$ functions shown in Figure 6 and $H_l$ defined by Equation 2. We can now generalize the analysis of the statistical uncertainty of integrals of buoy spectra given by Young (1986) to provide estimates of uncertainties for these spatially filtered wave heights.

Indeed, Young (1986) showed that any integral of the wave spectrum $E(f)$, obtained from a time series is $\chi^2$-distributed. For the particular case of the zeroth moment of the wave spectrum $m_0 = \int E(f)df$ the number of degrees of freedom is related to the record length $\tau$, $\nu_{m0}(\tau) = 2\tau/Q_f^2$ with the spectral frequency peakedness $Q_f$ defined by

$$Q_f^2 = \frac{\int_0^\infty E^2(f) df}{\left(\int_0^\infty E(f) df\right)^2} \qquad (16)$$

Because the buoy estimate of the significant wave height is $\hat{H}_\tau = 4\sqrt{m_0}$, it implies that $\hat{H}_\tau$ is $\chi$-distributed, and, assuming error-free measurements in the time series,

$$\frac{\text{std}(H_\tau)}{\text{mean}(H_\tau)} = \sqrt{\frac{\Gamma^2(\nu_{m0}(\tau)/2)\nu_{m0}(\tau)}{2\Gamma^2((\nu_{m0}(\tau)+1)/2)}} - 1 \approx 0.5Q_f/\sqrt{\tau}, \qquad (17)$$

Where $\Gamma$ is the Euler gamma function.

Following Young (1986), if we had a perfect spatial mapping of the surface elevation $\zeta(x, y)$ over a square of side length $L$, then $\hat{H}_L$ is a $\chi$-distributed random variable with $\nu_{m0}(L) = 2L^2/Q_{kk}^2(2\pi)^2$ degrees of freedom, giving the uncertainty

$$\frac{\text{std}(H_L)}{\text{mean}(H_L)} = \sqrt{\frac{\Gamma^2(\nu_{m0}(L)/2)\nu_{m0}(L)}{2\Gamma^2((\nu_{m0}(L)+1)/2)}} - 1 \approx \pi Q_{kk}/L. \qquad (18)$$
$Q_{kk}$ is given by Equation 4, and is analogous to $Q_f$ but defined from the double-sided wavenumber spectrum $E(k_x, k_y)$, instead of the single-sided frequency spectrum $E(f)$.

In our example, with $Q_{kk} = 43$ m and $Q_f = 4$ s$^{0.5}$, a standard 20 min buoy record gives a relative uncertainty std ($H_r/H_s = 0.058$, and it would take a square of side length $L = 2.4$ km to obtain the same uncertainty.

More generally, the same relative uncertainty is given by equating Equations 17 and 18 giving the spatio-temporal equivalence between observations of spatial scales $L$ and time scale $\tau$,

$$L = 2\pi \frac{Q_{kk}}{Q_f} \sqrt{\tau}. \qquad (19)$$

However, for an altimeter single measurement, our simulations in Figure 9e give a relative uncertainty of 0.085 that is equivalent to a square side $L = 1.6$ km. That scale is about $\rho_C/2$, and thus covers the same area as a disk of radius $\rho_C/(2\sqrt{\pi})$.

Alternatively, this result can be obtained by integrating the PSD of $H_l$. In the limit of statistics taken over scales $d_1$ much larger than $\rho_C$, eq. (step 4 bis) in DC23 gives,

$$\text{std}(\hat{H}_s) \approx \frac{\sqrt{2(4 - \pi)Q_{kk}H_s}}{4.5\rho_C} \approx \frac{4.2Q_{kk}\sqrt{H_s}}{\sqrt{h}} \qquad (20)$$

The combination of $L \approx \rho_C/2$ with Equation 18 gives the same result with the factor 4.2 replaced by 4.4.

Because altimeter data is generally averaged or filtered along-track in order to reduce the uncertainty of the measurements (Schlembach et al., 2020), we will now examine the uncertainty of the resulting along-track averages. For this we first need to investigate along-track correlations and define an effective resolution.

### 3.3. Along-Track Correlation and Effective Resolution of $\hat{H}_s$

The best retracker for sea level in Figure 9 is the one that will give the smallest values of $\hat{z}_e$, hence the LS retracker with results shown in (b). For the wave heights, it is unclear what are the retracking options that give the most accurate representation of the variability of the local wave $H_l$. Here we propose that wave heights $\hat{H}_s$ should be as similar as possible to the along-track sampling of the smoothed local wave height $H_{l,\sigma_l}$, obtained from $H_l$ with a two-dimensional Gaussian smoothing with a parameter $\sigma_l$. A reasonable expectation, consistent with the power spectrum of along-track $\hat{H}_s$ (Figure 10a), is that the “reasonable truth” is given by filtering with a scale $\sigma_l \approx \rho_C/4$, which is $\sigma_l = 800$ m in the example considered here. We may dream of being able to resolve smaller details, that can be found for example, in $H_{l,\rho_C/10}$, corresponding to $\sigma_l = 320$ m, but that “dream” is out of reach of altimeters, given the sensitivity kernels $\hat{J}_H$ for different retracking options, as shown in Figure 6. Indeed, the least squares cost function leads to a maximum sensitivity at $b = 0.25$ which corresponds to a distance from nadir $\rho = \rho_C/2$ where $\hat{J}_H$ is maximum.

In terms of sea level, instead of a zero value which is our input to the simulation, the retracked sea level exhibits a plateau at wavenumbers under 0.2 cpk, corresponding to the well-known 'hump' in along-track sea level spectra (Dibarboure et al., 2014). In terms of wave heights, considering the LS retracking, the black curve in Figure 10a, we find that its spectral level is similar but slightly higher than the “reasonable truth”, consistent with the analysis by DC23 who found the same variance as the “truth” when using $\sigma_l = \rho_C/4.5$. The ML retracking was adjusted, with $r_{min} = 0.06$ to give a similar PSD, while maximizing the coherence with the reasonable truth. This is detailed in Appendix C.

The coherence and phase are useful to quantify the distortion effect caused by the maximum of the $J_H$ function away from $b = 0$, and its change of sign. In along-track spectra this leads to scales for which the coherence goes to zero for LS retracking, here at $k \approx 0.22, 0.5$, and $0.71$ count per km (cpk), corresponding to wavelengths $L \approx 1.5, 0.7$ and $0.46$ times $\rho_C$. For these wavenumbers, the coherence phase jumps from near zero to near $180^\circ$ and back. As a result, LS retracking does not reproduce correctly any detail at wavelengths shorter than about $1/(0.22$ cpk) $\approx \rho_C/0.7$. $H_l$ perturbations at these scales have inverted signs when oriented across-track, and the along-track spectra also contain the projection of shorter scales on the along-track direction. With a two-dimensional map
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Along-track spectra of (a) $\hat{H}_s$ and (b) $\hat{z}_e = -\hat{t}_c/2$ estimates with different retracking options in the case of noiseless waveforms. (c) and (d) shown their coherence and phase shift relative to a “reasonable truth” $H_{l,\rho/4}$ obtained by filtering the local wave height $H_l$ with a 2-dimensional Gaussian filter of width $\sigma_{l,4} = \rho_c/4$ centered at the nadir point, while the “alternative truth” that contains much smaller detail is $H_{l,\rho/10}$.</figcaption></figure>

shown in Figure 9 one could invert the true map by a deconvolution with the $\hat{J}_H$ filter such as the example shown in Figure 7c. In one dimension, the phases of the different cross-track perturbations scales are mixed up. The shape of the $\hat{J}_H$ filter gives a stronger 2D response for LS compared to ML at short components, down to $L = \rho_c/2.5$ ($k = 0.7$ cpk) and those perturbations that are not in the along-track direction project on the x-axis at lower wavenumbers (larger scales): as a result LS retracking contains many along-track fluctuations associated to $H_l$ cross-track patterns, which are mixed up with “true alongtrack variability”, explaining the lower coherence of LS for $k < 0.2$ cpk. Assuming that we need at least two independent measurements per resolved wavelength, we may define an effective along-track resolution $\rho_{eff} = \rho_c/\alpha$. At this stage we expect that $\alpha \approx 1.4$, with different values for different retracking methods.

### 3.4. Uncertainty of Averaged Estimated Wave Heights $\bar{H}_s$

Now that we know how $\hat{H}_s$ estimates are correlated along-track we can estimate the uncertainty of $\bar{H}_s$, the average of $n (\langle \cdot \rangle_n)$ consecutive values of $\hat{H}_s$. For independent measurements this reduction would be a factor $1/\sqrt{n}$, but because the successive footprints overlap, there is only a $\sqrt{n_f/n}$ reduction where $n_f$ is the number of data points per effective footprint,

$$n_f = \frac{\rho_{eff}}{V_n/f_s} \approx \frac{\sqrt{2H_s h}}{\alpha V_n/f_s}, \qquad (21)$$

Where $V_n \approx 7$ km/s is the velocity of the satellite nadir on the ground, $\alpha$ is an along-track de-correlation parameter, $f_s$ is the sampling frequency of the measurement, typically $f_s = 20$ Hz for most altimeters. In practice we have found $\alpha \approx 1.5$ to be a good approximation for both ML and LS retracking, with possibly a weak dependence on $Q_{kk}$ and a range of possible values from 0.5 to 2.

This gives an expected variance of $(\bar{H}_s)$ caused by wave groups,

$$\text{var}_{\text{wg}}(\bar{H}_s) = \text{var}(\langle \hat{H}_s \rangle_n) \approx \frac{4.2^2 Q_{kk}^2 n_f H_s}{nh}, \qquad (22)$$
<figure><img src="image_13.png" /><figcaption>Figure 11. Map of the mean values of $Q_{kk}$, in meters, simulated for January–July 2023.</figcaption></figure>

All these calculations assumed noise-free measurements, but the interference of radar waves causes speckle noise, just like the interference of waves make groups. Speckle gives an extra term in the cost function that is a sum of $\chi^2$-distributed independent variables, and thus also $\chi^2$-distributed, as detailed in Appendix A.3 of DC23. The corresponding variance of fluctuations induced by speckle noise is given by,

$$ \mathrm{var}_s(\bar{H}_s) = s \times H_s/n \qquad (23) $$

With $s$ a function of the number of pulses $N_p$ per measurement. For least-square fitting with broad antenna patterns DC23 found

$$ s \approx s_0/N_p, \qquad (24) $$

With $s_0 = 5$ m. This expression gives $s = 0.019$ m for the LS fit of CFOSAT waveforms. There is a priori no such simple expression for ML retracking. It should be noted that the variance induced by speckle averages out faster than the wave group effect, like $1/n$, without the correction factor $n_f$. We note that for wave heights under 3 m, the speckle effect is further influenced by the discretization of the waveform and typically gives higher values of $s_0$.

We may assume that both effects are uncorrelated giving a total variance,

$$ \mathrm{var}(\bar{H}_s) = \mathrm{var}_{\mathrm{wg}}(\bar{H}_s) + \mathrm{var}_s(\bar{H}_s). \qquad (25) $$

## 4. Verification Over a Wide Range of Simulated Sea States

Although we looked in detail at a single and very particular sea state, we expect that our uncertainty model is applicable to any sea state, which is uniquely characterized by two parameters: the significant wave height $H_s$ and the wavenumber spectral peakedness $Q_{kk}$. The uncertainty model is also a function of the satellite instrument configuration through the altitude $h$ and number of pulses averaged $N_p$.

Given the prominent role of the peakedness, it is interesting to show the expected variability of $Q_{kk}$. We have used a 0.5° resolution global WAVEWATCH III model configuration with wave generation and dissipation source terms parameterized following the T702GQM option described in Alday and Ardhuin (2023), using a quasi-exact calculation of wave-wave interactions (Gagnaire-Renou et al., 2010; Lavrenov, 2001). As expected from its definition, $Q_{kk}$ is much larger for swells than for wind seas, and generally larger for longer dominant periods. Figure 11 gives average values of $Q_{kk}$ over a time period corresponding to the Austral summer and Fall.

We chose that time frame to minimize the effect of sea ice in the Southern Ocean: the presence of sea ice strongly damps the shorter wave components, leading to very large values of $Q_{kk}$ and very small wave heights. Outside of ice-covered regions, $Q_{kk}$ is typically under 10 m for enclosed seas and fetch-limited regions, and increases to 15–
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Distribution of modeled peakedness $Q_{kk}$ for ice-free conditions for January–July 2023, against usual parameters (a) $H_s$ and (b) Energy period $T_{m0,-1}$. (c) Compares $Q_{kk}$ to the peakedness for the frequency spectrum $Q_f$.</figcaption></figure>

40 m in the swell pool of the Eastern Pacific. Besides these mean values there is a significant variability, with a general increase with wave height shown in Figure 12a. Among the usual sea state parameters, $Q_{kk}$ is best correlated to the square of the so-called energy period $T_{m0,-1}$ (Figure 12b). When comparing the uncertainty of wave heights from buoy measurements and satellite data, it is also useful to know that there is not a simple correspondence between $Q_{kk}$ and $Q_f$ (Figure 12c).

We have simulated waveforms for 250 different sea states selected to fill a gridded histogram of $H_s$ and $Q_{kk}$. We insist that our sea state selection maximizes the ranges of $H_s$, from 0.5 to 12.5 m, and $Q_{kk}$ varying from 3 to 110. Most of these selected sea states are extremely unlikely, as shown in Figure 12a. A first display of the variability for the sea level and wave height is shown as a function of the wave height in Figure 13. In each panel, each dot corresponds to a different sea state with a given value of $H_s$ and $Q_{kk}$. For each dot, 11,000 waveforms were simulated from the same sea surface, shifting the nadir position (as we did for Figure 9) and retracked. Waveforms were simulated with and without speckle noise, and each was retracked with both LS and ML cost functions, using $r_{min} = 0.06$. The variability generally increases with wave height. For sea level, in panels a and b, it is of the order of 1%–3% of $H_s$, with some enhancement caused by speckle noise. We note that ML-based retracking is more noisy than LS for $\hat{z}_e$ retrieval, with occasional outliers. For wave heights, in panels c and d, the variability is generally higher with LS retracking once speckle noise is taken into account (panel d). We have found that even for ML fitting, we may use Equation 23 with a variance reduced by a factor 5 compared to LS fitting, as shown in Figure 13e.

At any given wave height, the variability can take values that differ by a factor 4 or more, as we expect from our analysis and the range of possible $Q_{kk}$ values. We verify our uncertainty model given by Equation 20 by plotting the uncertainties, now normalized by $H_s$ to a power $p$ against $Q_{kk}$ in Figure 14. We first note that the sea level, Figures 14a and 14b, scales with $p = 1$, probably due to the zero average of $J_z$, but we have not investigated this question further. For wave heights, Equation 20 gives a good representation of the data without speckle noise (see Figure 14c), and Equation 25 is a good model for the full simulation that includes noise (see Figure 14d), both uncertainties scaling with $p = 1/2$. In both cases there is an underestimation of the variability for high values of $Q_{kk}$.

We finally estimate along-track averages of 20 consecutive values to simulate 1 Hz averages. In Figure 15, we compare our error models given by Equations 22 and 23 to the variation of 1 Hz average simulation outputs. It shows that the error model given by Equation 22 is generally correct ($R^2 = 0.99$ for the selected sea states).

**5. Discussions and Applications**

The uncertainty model proposed in Equations 22 and 23 and verified with Figures 14d and 14h, appears robust, and is explained by the correlation structure that we understand well for the Least Square cost function. It also seems to hold for our adaptation of the Maximum Likelihood cost function. Some persistent biases may be refined. For example, the speckle contribution is underestimated for large wave heights. This is possibly due to the use a broad antenna pattern in DC23: for the largest wave heights and narrow radar beams the different shape of the waveform will give a different value of $s_0$, which can possibly be obtained analytically or numerically. Another bias is found for $Q_{kk} > 50$, with an overestimation of the wave group contribution given by Equation 20. In these cases the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Variability of (a, b) estimated sea level $\hat{z}_e$, and against $Q_{kk}$ (c, d) wave height $\hat{H}_s$ against $H_s$. (e) Speckle contribution to the variance of wave height.</figcaption></figure>

spectrum of the surface envelope must be very narrow, possibly narrower than the altimeter transfer function (the Fourier transform of $J_H$), and the approximation proposed in DC23, that the envelope PSD is constant, is likely to overestimate the variability of $\hat{H}_s$. This may be corrected by computing the spectral convolution (Step 3 in DC23), or using a better approximation for the envelope spectrum, not as a constant but for example, a two-dimensional Gaussian function. For our objective, very few conditions are concerned as the median value of $Q_{kk}$ is under 60 m, even for wave heights up to 18 m. This will be a real issue when extending the present work to Delay-Doppler altimetry, as the effective footprint becomes very small in the along-track direction.

## 5.1. Other Satellite Missions

We are now in a position to propose a clear trade-off between precision and resolution for storm conditions for CFOSAT data, and possibly extrapolate this to other satellite missions, as illustrated in Figure 16. And it would be interesting to compare the measured values of $\text{std}(\hat{H}_s)$ at cross-overs from different missions.

In order to arrive at the same uncertainty level as the buoy data, we find that we need to average around $n = 3$ points with CFOSAT's 4.5 Hz sampling, and $n = 12$ points at 20 Hz for the storm case considered here, for both cases this is an along-track length of 4.2 km. Due to the scaling of the effective footprint with $\rho_C = \sqrt{2hH_s}$, the lowest altitude of CFOSAT allows it to have a higher resolution with $h = 519$ km, compared to the $h = 1,336$ km
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Variability of (a, b) estimated sea level $\hat{z}_e$ against $Q_{kk}$ (c, d) wave height $\hat{H}_s$ against $H_s$. This is the same data as in Figure 13 but rescaled and plotted against different variables. Smaller dots correspond to cases with $H_s < 2$ m.</figcaption></figure>

of Jason 3. The sampling $f_s = 4.5$ Hz of CFOSAT is particularly efficient, with measurements that are more independent than with $f_s = 20$ Hz. For CFOSAT the available time between independent samples is well used by scanning the ocean with off-nadir beams to measure the wave spectrum (Hauser et al., 2021) that can be used to estimate $Q_{kk}$, as in DC23, and other properties useful to interpret nadir altimetry such as the skewness and the slope-sea level correlations (Janssen, 2014; Srokosz, 1986). Future missions can use the same type of nadir + off-nadir design to also measure ocean currents (Ardhuin, Brandt, et al., 2019).

Even for a low wave height of 1 m, at the Jason 3 altitude, the 20 Hz data is useful for reducing speckle but otherwise produces highly correlated errors. The value $n_f = 3$ means that, without speckle, a 6 Hz sampling would be enough to sample the variability induced by wave groups. ML retracking, with $s_0 \approx 1$ m (see Figure 16b) can also be used to reduce noise levels, in particular for wind seas with low wave heights (solid and dotted lines). However, when data are averaged over 1 Hz, the speckle contribution is less important, especially for swell-dominated conditions (dashed lines in Figure 16c). In that case a higher orbit provides averaging over a wider area, both along-track and across-track.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Variability of along-track 1 Hz averages $\bar{H}_s$ plotted against our predicted variance from Equations 22 and 23. The black line shows the $x = y$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Resolution and uncertainty for wave height measurements extrapolated to other satellite missions, taking into account their different altitudes and acquisition chronogram (number of pulses per burst $N_p$ and number of bursts per second $f_s$), but neglecting the effect of the antenna aperture: (a) spatial resolution estimated as $1.5\rho_c$ (b) normalized standard deviation of the measurement, (c) normalized standard deviation of 1 Hz along-track averaged measurements. We have used typical values of $Q_{kk} = 2H_s$ for wind-seas and a higher value for very long periods or narrow swells $Q_{kk} = 60$. The speckle parameter $s_0 = 5$ m is typical of least squares (LS) fitting, while $s_0 = 1$ m corresponds to maximum likelihood (ML).</figcaption></figure>

## 5.2. Re-Defining Significant Wave Heights

Looking back at Figure 3, there is a need for defining the underlying wave height from fluctuating measurements. The obvious solution is to average the data along-track and estimate the precision of the average using our uncertainty model. We give here two examples.

In the case shown in Figure 2, the maximum estimate $\hat{H}_s$ is 11.8 m, using LS fitting. When averaging over the 50 bursts, and considering our sampling error model we get the following estimate of the underlying true wave height (removing speckle and wave group effects), $H_s = 9.2 \pm 0.3$ m.

Hanafin et al. (2012) reported the highest-ever wave height measurement at $\bar{H}_s = 20.1$ m, using a Jason-2 data over storm Quirin on 14 February 2011 at 11:05 UTC, in the North Atlantic, with a relative precision $\langle\text{std}(\hat{H}_s)/\hat{H}_s\rangle = 8.9\%$ for the neighboring values. This is a 1 Hz-averaged data. Due to a different retracker, called WHALES, the maximum value for this event was revised at $\bar{H}_s = 19.7$ m in the version three of the Sea State CCI dataset (Schlembach et al., 2020), with a relative precision of 6.4%. Based on Figure 6, we expect the WHALES retracker to provide an effective resolution $\rho_{\text{eff}}$ in between the ML and LS retrackers, so that the uncertainty model, Equations 20 and 23 should apply. We thus expect the effective Jason resolution to be close to 4.5 km. Without a specific wave model hindcast of that storm we may expect $Q_{kk} \approx 60$ based on Figure 12a. With that value, our uncertainty model, Equations 20 and 23 using $N_p = 90$ and $h = 1,336$ km, gives $\text{std}(\hat{H}_s) = 1.43$ m for a single 20 Hz estimate. For reference, the value provided in the CCI dataset is $\text{std}(\hat{H}_s) = 0.58$ m. That value is anomalously low compared to the neighboring 1 Hz record with the following sequence of 9 values centered on the record with the maximum wave height $\text{std}(\hat{H}_s) = 1.43, 1.35, 1.31, 1.08, 0.58, 1.09, 1.3, 1.6, 1.1$ m, corresponding to 1 Hz averages $\bar{H}_s = 19.7, 17.6, 18.8, 19.3, 19.7, 17.6, 17.2, 18.3$ and $17.8$ m. The number of valid waveforms was also minimum (13 out of 20) for that record with the lowest variability.

Our model uncertainty for the 1 Hz average, Equations 22 and 23 gives 0.90 m, or about 5% of the measurement, with wave groups alone accounting for 0.87 m. We may average over a longer distance to get a mean value of $H_s$ and the corresponding uncertainty. Averaging over 54 km (9 points at 1 Hz) reduces the uncertainty to 0.29 m and gives an average of 18.5 m.

Hence, what should be reported as the maximum value of $H_s$? Is it $19.7 \pm 0.9$ m, from the 1 Hz record, or $18.5 \pm 0.3$ m from the 54 km average? From our analysis the first number is likely to be strongly impacted by wave groups: it may be correct for defining a local wave height that is physically correct, there is indeed a region with very high waves over a few kilometers of the satellite track, just like on Figure 3. However if we want to compare to numerical wave models that ignore wave groups, the longer along-track average is a better choice.
Alternatively, filtering small-scale variations in $\bar{H}_s$ can be done using Empirical Mode Decomposition (Dodet et al., 2020; Quilfen et al., 2018). That procedure gives $H_s = 18.7 \pm 0.3$ m, a value also reported in the CCI dataset, which is consistent with our estimate. Further analysis of other storm events will be useful for better understanding of the output of denoising using Empirical Mode Decomposition (Quilfen et al., 2018).

## 6. Conclusions and Perspectives

Following the demonstration by De Carlo et al. (2023) that the sampling uncertainty in the presence of wave groups is a significant source of along-track fluctuations in altimeter measurements, we have explored how we may interpret these fluctuations and define an uncertainty for the underlying true significant wave height. Our argument is that the contribution of wave groups to the local wave height should be removed when estimating a significant wave height due to their fast propagation: they are not relevant for most applications. That approach is consistent with phase-averaged wave modeling in which wave group fluctuations are absent. We have confirmed the analysis by DC23 for a wide range of realistic waveforms and retracking methods: the amplitude of small-scale fluctuations caused by wave groups is proportional to the peakedness parameter $Q_{kk}$ and the square root of the wave height. These fluctuations are spatially correlated through the effective footprint width that can be approximated as $\rho_c/1.5$, with some small differences depending on the details of the retracking method. This provides a useful scale to count the number of independent data in a satellite segment. The along-track distance $1.5\rho_c$ is also a good estimate of the shortest wavelength that can be resolved in the spatial pattern of the local wave height, including wave groups when they are present. This finest resolution is achieved when using some form of Maximum Likelihood cost function that is more sensitive than the least squares cost function to perturbations near nadir, and the weighted least squares used in WHALES provides an interesting intermediate method. These prediction could be tested with cross-overs between Jason-3 and SWOT which carries a nadir Poseidon-3 altimeter that is a copy of the Jason-3 instrument, but at a different altitude. It should also be possible to see that speckle noise is decorrelated between measurements from satellites flying in tandem with a 30 s time separation, whereas the effect of wave groups should be persistent (Rieu et al., 2021).

Our implementation of a Maximum Likelihood cost function may provide more realistic estimates of wave heights, but it generally led to larger errors in the sea level. The WHALES retracker is an interesting candidate for obtaining both accurate sea level and wave heights. An alternative approach was timidly explored in Appendix D: one may add more degrees of freedom to the waveform shape to properly handle their more complex shapes, including wave groups and skewness effects.

For very large wave heights, say $H_s > 15$ m, we find that the effective altimeter along-track resolution is of the order of 6 km or more, depending on the satellite altitude. Any estimate of wave heights with an accuracy of 3% or better typically requires along-track averaging or spatial filtering methods. With this kind of post-processing (averaging or filtering), the effect of speckle noise is less important, and we might even make a meaningful use of the C-band instruments that are also present on most satellite altimeters in addition to the Ku-band data that were discussed here.

The effects of wave groups on Delay-Doppler altimetry are not obvious a priori, and will require a dedicated investigation. As noted by Moreau et al. (2018), the anisotropic measurement geometry of Delay-Doppler altimetry introduces the difficulty that narrow directional swells are now part of the sea level fluctuations when propagating along-track, whereas they are still very much part of the sea state when propagating cross-track. That difficulty may be leveraged to provide some advantage, for example, for swell detection (Altiparmaki et al., 2022; Collard et al., 2022).

## Appendix A: Derivation of a Theoretical Waveform

In the following we shall use the same notations as in Tourain et al. (2021). We generalize the usual approach by Brown (1977), allowing the vertical distribution of scattering elements, denoted as PDF, to be a function of both the horizontal distance to nadir $\rho$ and the time $t$. In practice we start with Gaussian surface elevation PDF with a standard deviation $\sigma_H$, which translates to a standard deviation in the arrival time of the echo $\sigma_s = 2\sigma_H/c$, with $c$ the speed of light. At nadir, $\rho = 0$, the epoch $\tau$ defines the local mean sea level and we have,
$$PDF(\rho = 0, t) = G(\sigma_s, \tau, t) = \frac{e^{-(t-\tau)^2/2\sigma_s^2}}{\sigma_s\sqrt{2\pi}}. \qquad (A1)$$

Off-nadir this generalizes to (Chelton et al., 1989)

$$PDF(\rho, t) = G(\sigma_s, \tau + t_p, t) \qquad (A2)$$

Where $R_E$ is the Earth radius, $h$ is the satellite altitude, and the radius-dependent time shift is

$$t_p = \frac{\rho^2}{ch}(1 + h/R_E). \qquad (A3)$$

This gives the theoretical waveform as

$$S(t) = PDF(\rho = 0, t) * FSSR(t) * PTR(t). \qquad (A4)$$

Where FSSR stands for Flat Sea Surface Response and PTR stands for Point Target Response. The first convolution represented by the symbol * is in fact the sea surface response SSR obtained from the integral over the distance from nadir $\rho$ (Brown, 1977). We introduce a local perturbation of $\sigma_s$ which becomes $\sigma_s$ ($1 + \Delta$) and this perturbation affects an area $A_0$ of the ocean centered at the distance from nadir $\rho_0$, which correspond to a range $h + R_0$ in the absence of waves. DC23 showed that the PDF ($\rho, t$) could be assumed Gaussian for each value of $\rho$. The perturbed surface elevation pdf is now.

$$PDF'(\rho,t) = PDF(\rho,t) + \frac{A_0\delta(\rho - \rho_0)}{2\pi\rho_0} [G((1 + \Delta)\sigma_s, \tau, t - t_p) - G(\sigma_s, \tau, t - t_p)] \qquad (A5)$$

$$\approx PDF(\rho,t) + \frac{\Delta A_0 \sigma_s \delta(\rho - \rho_0)}{2\pi\rho_0} \frac{\partial G(\sigma_s, \tau, t - t_p)}{\partial \sigma_s},$$

$$= PDF(\rho,t) + \frac{\Delta A_0 \delta(\rho - \rho_0)}{2\pi\rho_0} G(\sigma_s, \tau, t - t_p) \frac{(t - \tau - t_p)^2 - \sigma_s^2}{\sigma_s^2}. \qquad (A6)$$

In the usual expressions, $\rho$ is transformed to a time $t_p$ with the following expression on the sphere

$$\frac{dt_p}{d\rho} = \frac{2\rho(1 + h/R_E)}{ch}. \qquad (A7)$$

Using this relation, we may now replace $\delta(\rho - \rho_0)$ by $\delta(t_p - t_{p0}) dt_p/d\rho$, to get the perturbed pdf as a function of two time scales,

$$PDF'(t_p, t) = G(\sigma_s, t) + a\delta(t_p - 4b\sigma_s) p(t - t_p), \qquad (A8)$$

With $a$ and $b$ defined by Equations 8 and 9.

$$p(t) = \frac{1}{\sqrt{2\pi}} e^{-0.5 \left(\frac{t-\tau}{\sigma_s}\right)^2} \left[ \left(\frac{t-\tau}{\sigma_s}\right)^2 - 1 \right]. \qquad (A9)$$

The dimensionless parameter $a < 1$ is the product of the relative wave height change and the ratio of the area $A_0$ affected by that change and an equivalent footprint area $A_{eq} = \pi c \sigma_s h / (1 + h/R_E)$ which is close to $\pi(\rho_c/2)^2$ or one quarter of the area of the oceanographic footprint defined by Chelton et al. (1989).
The convolution of FSSR(t) and PDF'(t) corresponds to an integration over the time $t_p$, which is Equation 2 in Brown (1977). It is thus the sum of two parts, the unperturbed part, and the perturbation given by the second term in Equation A8,

$$SSR(t) = A\sigma_0 \left\{ ap(t - 4b\sigma_s) + \frac{1}{2}[1 + \operatorname{erf}(u(t))] \right\} e^{-v(t)} + N_t, \quad (A10)$$

With.

$$u(t) = (t - \tau - c_{\xi}\sigma_s^2)/(\sqrt{2}\sigma_s), \quad (A11)$$

$$v(t) = c_{\xi}(t - \tau - c_{\xi}\sigma_s^2/2), \quad (A12)$$

$$A = \exp(-4\sin^2 \xi/\gamma), \quad (A13)$$

$$c_{\xi} = \frac{4c(4mss \cos(2\xi) + \gamma)}{h4\gamma mss} \quad (A14)$$

$$\gamma = \frac{2}{\ln(2)} \sin^2(\theta_{3dB}), \quad (A15)$$

Where $\xi$ is the antenna mispointing angle, $\theta_{3dB}$ is the antenna pattern parameter, $N_t$ is the mean thermal noise, and mss is the mean square slope (Tourain et al., 2021).

Equation A10 corresponds to a modification of the adaptive model in Tourain et al. (2021) with the perturbation function $p(t - 4b\sigma_s)$ that is the difference of two Gaussian PDFs centered at $t = \tau + 4b\sigma_s$, with standard deviation $\sigma_s$ and $(1 + \Delta)\sigma_s$.

In the case $b = 0$, we note that Equation A10 is equivalent to the effect of surface elevation skewness derived by Hayne (1980), with $\lambda = 6a$, re-derived by Srokosz (1986) and used by Gómez-Enri et al. (2007).

The full waveform is finally obtained by convolution with the instrument PTR. In the absence of more information we have used,

$$PTR(t) = \mathrm{sinc}^2(\pi B t), \quad (A16)$$

Giving the waveform

$$S(t) = SSR(t) * PTR(t). \quad (A17)$$

## Appendix B: Waveform Simulation, Retracking and Verification

The waveforms are obtained from a realization of the sea surface elevation map using random phases over 4,096 by 4,096 points with a horizontal resolution of 14 m, hence covering 56 by 56 km. Each point of the surface is given a radar power based on the 2-way antenna pattern, and the waveform is given by the power-weighted histogram of the distribution of travel times between each point of the surface and the satellite. When speckle noise is included it corresponds to $N_p \times 320/400$ looks, which is the number of pulses per cycle, $N_p = 264$ for CFOSAT (Tourain et al., 2021), corrected for the resampling factor of the waveform, from 320 to 400 MHz. A “local significant wave height” is defined at each point as $H_l(x,y) = \sqrt{32/\pi \times \eta(x,y)}$ where $\eta$ is the surface envelope (De Carlo et al., 2023), so that the spatial mean of $H_l$ is the true underlying wave height $H_s$. The retracked wave height $\hat{H}_s$ and epoch $\tau$ are computed for discrete satellite positions on a two-dimensional grid with a resolution of 350 m, as if the ocean were sampled by 106 satellites flying side by side and with a waveform computed every 0.05 s (a rate of 20 Hz) along each track. The result is a map of estimated parameters.
<caption>Table B1</caption>

Statistics for Wave Height and Epoch, for a Surface With Strong Wave Groups. Starting From the Idealized Waveform Simulation Which Does Not Include the PTR at the Top, We Progressively Add the PTR (No Noise), Then Thermal Noise, Then the Speckle

<table><thead><tr><th rowspan="2">PTR</th><th rowspan="2">Thermal noise</th><th rowspan="2">Speckle</th><th rowspan="2"></th><th colspan="2">H&#x0302;<sub>s</sub> (m)</th><th colspan="2">z&#x0302;<sub>e</sub> (m)</th></tr><tr><th>Mean</th><th>Std</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td rowspan="2">X</td><td rowspan="2">X</td><td rowspan="2">X</td><td>ML, r<sub>min</sub> = 0</td><td>8.88</td><td>1.00</td><td>-0.122</td><td>0.24</td></tr><tr><td>LS</td><td>9.23</td><td>0.70</td><td>0.000</td><td>0.07</td></tr><tr><td rowspan="2">✓</td><td rowspan="2">X</td><td rowspan="2">X</td><td>ML, r<sub>min</sub> = 0</td><td>9.14</td><td>0.96</td><td>0.328</td><td>0.22</td></tr><tr><td>LS</td><td>9.27</td><td>0.70</td><td>0.001</td><td>0.07</td></tr><tr><td rowspan="2">✓</td><td rowspan="2">✓</td><td rowspan="2">X</td><td>ML, r<sub>min</sub> = 0</td><td>9.14</td><td>0.96</td><td>-0.045</td><td>0.22</td></tr><tr><td>LS</td><td>9.27</td><td>0.70</td><td>0.001</td><td>0.07</td></tr><tr><td rowspan="2">✓</td><td rowspan="2">✓</td><td rowspan="2">✓</td><td>ML</td><td>9.14</td><td>0.96</td><td>-0.045</td><td>0.23</td></tr><tr><td>ML, r<sub>min</sub> = 0.06</td><td>9.23</td><td>0.78</td><td>-0.009</td><td>0.11</td></tr><tr><td></td><td></td><td></td><td>LS</td><td>9.28</td><td>0.77</td><td>0.002</td><td>0.08</td></tr></tbody></table>

Figure 9 shows some examples of such maps for different simulation settings, and fitting with ML or LS cost functions. Statistics for the retrieved parameters shown in Figure 9, are summarized in Table B1. We find that the root mean square (rms) wave height is underestimated with ML compared to LS, and the standard deviation of wave heights is larger with ML compared to LS, consistent with the retracking of the true waveforms in Figure 2. We also note that the ML retrieved epoch and wave height are strongly correlated with *r* = 0.85, which is comparable to *r* = 0.81 in Figure 2. For the purpose of reducing the epoch noise, for example, taking *z*′<sub>e</sub> = *z*<sub>e</sub> − α(*H*<sub>s</sub> − *H*<sub>s</sub>), the LS data give lower noise residuals than the ML data.

## Appendix C: Influence of R<sub>min</sub>

When using ML retracking, one may optimize the contribution of the lowest range gates used in ML fitting. Here we investigate the influence of the choice of *r*<sub>min</sub>, and try to maximize the coherence with our “reasonable truth” given by *H*<sub>I,σ=ρc/4</sub> over the widest possible range of scales, while keeping a near-zero phase shift, and getting perturbations on the epoch as low as possible. The value *r*<sub>min</sub> = 0.06 that gave results similar to the CFOSAT in Figure 2 appears to be a good compromise. Figure C1 shows that lower values of *r*<sub>min</sub> will all produce more noise

<figure><img src="image_19.png" /><figcaption>Figure C1. Same layout as Figure 10, with additional lines for different values of *r*<sub>min</sub>.</figcaption></figure>
in the epoch. Higher values reduce the range of wavenumbers with high coherence, giving results closer to the LS retracking, with a wider range of short scales for which the retracked values are out of phase of the true perturbations.

## Appendix D: Consequences for Retracking

This paper dealt with existing data sets, already retracked with existing methods, but our results may be used to refine retracking methods and better interpret alternatives. On a basic level, it is possible that some averaging before retracking may provide more robust results. Also, the results on along-track correlations and uncertainty model may differ for methods we did not cover: for example, the use of range-dependent weights in the WHALES retracker (Schlembach et al., 2020) modifies the $\hat{J}_H$ function as shown in Figure 6. Alternatively when a skewness parameter is added to the set of fitting parameters, following Hayne (1980), it will catch the waveform distortion caused by wave groups near nadir. One could imagine adding more degrees of freedom to the waveforms with a sum of wave group contributions $a(b)$ for each range $b$, and ideally one may want to estimate these values of $a(b)$ for each discrete range, and inverting the black curve $H_b(b)$ in Figure 7e from the waveform in Figure 7f.

Since the possible adjustment to retracking methods are endless and best choices probably depend on the chosen application (e.g., characterizing wave group properties, reducing noise on sea level estimates ...) we will not go down this path here. Instead we just illustrate how a modified retracker may better fit the waveform: we have chosen 2 variants on the LS and ML retrackers (here termed LS2 and ML2 for clarity) used in Figure 2. In LS3 and ML3 we add a the skewness parameter $\lambda_{3,0,0}$ as defined by Srokosz (1986), which is the skewness of the surface elevation points of zero slope, and in our model waveform corresponds to 6 times the amplitude of wave group perturbations at nadir $\lambda_{3,0,0} = 6a$ ($b=0$). This is the approach followed by Callahan and Rodriguez (2004) and Gómez-Enri et al. (2007), with the minor difference is that we use SWIM L1B data in which the antenna pattern and power have been corrected for, so that we do not have to deal with the usual other unknowns that are the mispointing and $\sigma_0$. In LS4 and ML4 we allow the Brown waveform to have one wave group perturbation of amplitude $a$ but that can be at any range $b$: because the possible waveforms fits with LS3 and ML3 are a subset of those for LS4 and ML4, the fits are at least as good with that other option, in the cases where the minimization method found the global minimum.

Results in Figure D1 show the values of wave heights, and $\lambda_{3,0,0}$ (or $6 \times a$ for LS4 and ML4), and two waveforms corresponding to the two highest values of $\hat{H}_s$ in the native CFOSAT data, here corresponding to the 20th and 49th (last) waveform in that sequence. We recall that the acquisition rate is 4.5 Hz, so that the nadir positions for consecutive waveforms are separated by about 1.5 km.

The first clear outcome shown in panel (a) is that LS3 and ML3 give wave heights very close to LS2 and ML2, with slightly lower extremes, and LS4 and ML4 give much lower value for the extremes (waveforms 20 and 49) but often fail to converge to reasonable values (waveforms 39 for LS4, 16, 22, 23 ... for ML4): this is not proposed as a practical retracking method but as a tool to understand some of the parameter variations. In panel (b), $\lambda_{3,0,0}$ fluctuations (-1 to 2) are much larger than its mean value of 0.17 with LS3 and 0.09 with ML3. We note that the mean value of $\lambda_{3,0,0}$ from ML3 is consistent with the skewness of 0.08 estimated from the CFOSAT directional spectrum using the method of Srokosz (1986) as implemented in the ECWAM model (Janssen, 2014). We suspect that the larger mean value for LS3 is caused by speckle noise, and estimating $\lambda_{3,0,0}$ is probably better done by first averaging several waveforms before retracking.

For waveforms number 20 and 49, we may expect that there is some significant wave group contribution, with a maximum $H_l$ that could be close to nadir for #20 because of the stronger value of $\hat{H}_s$ with ML2 compared to LS2, based on the different shapes of $\hat{J}_H$ in Figure 6. When the fitting waveform is allowed to have some skewness effect, the value of $\hat{H}_s$ is reduced from 13.9 with ML2 to 12.3 m with ML3 as we are effectively removing the effect of wave groups at nadir, and $\hat{J}_{H,ML2}$ is sensitive to these. This is not the case when changing from LS2 to LS3. However, if we allow the wave group perturbation to be away from nadir, then $\hat{H}_s$ drops to 10.2 m with LS4 compared to 11.5 with LS2, and the optimal perturbation position is placed at $b=0.13$. Things are a little different for waveform #49, presumably because the perturbations are located further from nadir. In that case ML3 is not very different from ML2, but both LS4 and ML4 give a much lower wave height, at 9.5 and 9.2 m respectively, with $b=0.14$.
<figure><img src="image_20.png" /></figure>

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

<figure><img src="image_22.png" /><figcaption>**Figure D1.** Retracked values of (a) $\hat{H}_s$ and (b) $\lambda_{3,0,0}$ or $6 \times a$ using 2, 3 or 4 parameter retracers applied to the same CFOSAT waveforms (c) and (d) waveform number 30 and 49 in that sequence, and associated fits.</figcaption></figure>

Data Availability Statement

The L2 and L2P SWIM data sets used here corresponds to the files reprocessed in version 5.1.2 and made available by CNES on the ftp server of AVISO+ (ftp-access.aviso.altimetry.fr, directories cfosat/swim_l2_op05 and cfosat/swim_l2p_box_nrt/), accessible to anyone after registration. Surface simulation and retracking software was developed in Python using elements from the WHALES retracker provided by M. Passaro, and is available at from the permanent link https://doi.org/10.17882/97944. A more up to date package can be found at https://github.com/ardhuin/wavesALTI.

The L2 SWIM dataset used here corresponds to the files reprocessed by CNES (2020) in version 5.1.2 and made available by CNES on the ftp server of AVISO+ (ftp-access.aviso.altimetry.fr, directory cfosat/swim_l2_op05), accessible to anyone after registration.
The L2P SWIM dataset used here corresponds to the files reprocessed by CNES/CLS (2021) in version 1.2 and made available by CNES on the ftp server of AVISO+ (ftp-access.aviso.altimetry.fr, directory cfosat/swim_l2p_box_nrt/), accessible to anyone after registration.

The L2S SWIM dataset used in this paper corresponds to the files reprocessed by Ifremer/CERSAT (2022) in version 1.0 and available at https://data-cersat.ifremer.fr/projects/iwwoc/swi_l2s.

## Acknowledgments

All CFOSAT data are provided by courtesy of CNSA and CNES. This research was made possible by support from ESA as part of the Sea State CCI project. MDC was supported by a postdoctoral grant from the Centre National d'Études Spatiales (CNES). MDC would like to thank Jean-François Piolle from Ifremer for the data management support. We are very grateful to Annabelle Ollivier, Alejandro Bohe, Beatriz Molero, Adrien Nigou, Bertrand Chapron and Marcello Passaro for in depth discussions that have led to a clearer exposition of the results.

## References

* Abdalla, S., Janssen, P., & Bidlot, J.-R. (2011). Altimeter near real time wind and wave products: Random error estimation. *Marine Geodesy*, 34(3-4), 396-406. https://doi.org/10.1080/01490419.2011.585113

* Alday, M., & Arduhuin, F. (2023). On consistent parameterizations for both dominant wind-waves and spectral tail directionality. *Journal of Geophysical Research*, 128(4), e2022JC019581. https://doi.org/10.1029/2022JC019581

* Altiparmaki, O., Kleinherenbrink, M., Naeije, M., Slobbe, C., & Visser, P. (2022). Sar altimetry data as a new source for swell monitoring. *Geophysical Research Letters*, 49(7). https://doi.org/10.1029/2021GL096224

* Arduhuin, F., Brandt, P., Gaultier, L., Donlon, C., Battaglia, A., Boy, F., et al. (2019). Skim, a candidate satellite mission exploring global ocean currents and waves. *Frontiers in Marine Science*, 6, 124. https://doi.org/10.3389/fmars.2019.00209

* Arduhuin, F., Chapron, B., & Elfouhaily, T. (2004). Waves and the air-sea momentum budget, implications for ocean circulation modelling. *Journal of Physical Oceanography*, 34(7), 1741-1755. https://doi.org/10.1175/1520-0485%282004%29034%3C1741%3AWATAMB%3E2.0.CO;3B2

* Arduhuin, F., Rascle, N., Chapron, B., Gula, J., Molemaker, J., Gille, S. T., & Rocha, C. (2017). Small scale currents have large effects on wind wave heights. *Journal of Geophysical Research*, 122(C6), 4500-4517. https://doi.org/10.1002/2016JC012413

* Arduhuin, F., Stopa, J. E., Chapron, B., Collard, F., Husson, R., Jensen, R. E., et al. (2019). Observing sea states. *Frontiers in Marine Science*, 6, 124. https://doi.org/10.3389/fmars.2019.00124

* Arhan, M., & Ezraty, R. (1978). Statistical relations between successive wave heights. *Oceanologica Acta*, 1, 151-158.

* Brown, G. S. (1977). The average impulse response of a rough surface and its applications. *IEEE Journal of Oceanic Engineering*, 2(1), 63-67. https://doi.org/10.1109/JOE.1977.1145328

* Callahan, P. S., & Rodriguez, E. (2004). Retracking of jason-1 data. *Marine Geodesy*, 27(3-4), 391-407. https://doi.org/10.1080/01490410490902098

* Cazenave, A., Palanisamy, H., & Ablain, M. (2018). Contemporary sea level changes from satellite altimetry: What have we learned? What are the new challenges? *Advances in Space Research*, 62(7), 1639-1653. https://doi.org/10.1016/j.asr.2018.07.017

* Challenor, P. G., & Srokosz, M. A. (1989). The extraction of geophysical parameters from radar altimeter returns from a nonlinear sea surface. In S. R. Brooks (Ed.), *Mathematics in remote sensing* (pp. 256-268). Clarendon Press.

* Chelton, D. B., Walsh, E. J., & MacArthur, J. L. (1989). Pulse compression and sea level tracking in satellite altimetry. *Journal of Atmospheric and Oceanic Technology*, 6(3), 407-438. https://doi.org/10.1175/1520-0426(1989)006(0407:pcaslt)2.0.co;2

* CNES. (2020). Wave products from SWIM Level 2 (L2), version 5.1.2. [Dataset]. Retrieved from ftp://ftp-access.aviso.altimetry.fr/(after-registration,directory:cfosat/swim_l2_op05). AVISO+ with support from CNES

* CNES/CLS. (2021). Wave products from SWIM level 2+ box off nadir NRT (L2P), version 1.2. [Dataset]. Retrieved from ftp://ftp-access.aviso.altimetry.fr/. (after registration, directory: cfosat/swim_l2p_box_nrt). AVISO+ with support from CNES

* Collard, F., Marié, L., Nouguier, F., Kleinherenbrink, M., Ehlers, F., & Arduhuin, F. (2022). Wind-wave attenuation in arctic sea ice: A discussion of remote sensing capabilities. *Journal of Geophysical Research*, 127(7), e2022JC018654. https://doi.org/10.1029/2022JC018654

* De Carlo, M., Arduhuin, F., Ollivier, A., & Nigou, A. (2023). Wave groups and small scale variability of wave heights observed by altimeters. *Journal of Geophysical Research*, 133(8), e2023JC019740. https://doi.org/10.1029/2023JC019740

* Dibarbouré, G., Boy, F., Desjonqueres, J. D., Labroue, S., Lasne, Y., Picot, N., et al. (2014). Investigating short-wavelength correlated errors on low-resolution mode altimetry. *Journal of Atmospheric and Oceanic Technology*, 31(6), 1337-1362. https://doi.org/10.1175/JTECH-D-13-00081.1

* Dodet, G., Piolle, J.-F., Quilfen, Y., Abdalla, S., Accensi, M., Arduhuin, F., et al. (2020). The sea state cci dataset v1: Towards a sea state climate data record based on satellite observations. *Earth System Science Data*, 12(3), 1929-1951. https://doi.org/10.5194/essd-12-1929-2020

* Gagnaire-Renou, E., Benoit, M., & Forget, P. (2010). Ocean wave spectrum properties as derived from quasi-exact computations of nonlinear wave-wave interactions. *Journal of Geophysical Research*, 115(C12), C12058. https://doi.org/10.1029/2009JC005665

* Gómez-Enri, J., Gommenginger, C. P., Srokosz, M. A., Challenor, P. G., & Benveniste, J. (2007). Measuring global ocean wave skewness by retracking RA-2 Envisat waveforms. *Journal of Atmospheric and Oceanic Technology*, 24(6), 1102-1107. https://doi.org/10.1175/JTECH2014.1

* Hanafin, J., Quilfen, Y., Arduhuin, F., Sienkiewicz, J., Queffeulou, P., Obrebski, M., et al. (2012). Phenomenal sea states and swell radiation: A comprehensive analysis of the 12-16 february 2011 North Atlantic storms. *Bull. Amer. Meterol. Soc.*, 93(12), 1825-1832. https://doi.org/10.1175/BAMS-D-11-00128.1

* Hauser, D., Tourain, C., Hermozo, L., Alraddawi, D., Aouf, L., Chapron, B., et al. (2021). New observations from the SWIM radar on-board CFOSAT: Instrument validation and ocean wave measurement assessment. *IEEE Transactions on Geoscience and Remote Sensing*, 59(1), 5-26. https://doi.org/10.1109/tgrs.2020.2994372

* Hayne, G. (1980). Radar altimeter mean return waveforms from near-normal-incidence ocean surface scattering. *IEEE Transactions on Antennas and Propagation*, 28(5), 687-692. https://doi.org/10.1109/tap.1980.1142398

* Ifremer / CERSAT. (2022). Global Ocean directional wave parameters level 2S from SWIM onboard CFOSAT for IWWOC project, version 1.0. [Dataset]. https://doi.org/10.12770/12cfed8d-7645-442b-b8ef-a8d08decbaed. Ifremer, Plouzane, France

* Janssen, P. A. E. M. (2014). On a random time series analysis valid for arbitrary spectral shape. *Journal of Fluid Mechanics*, 759, 236-256. https://doi.org/10.1017/jfm.2014.565

* Lavrenov, I. V. (2001). Effect of wind wave parameter fluctuation on the nonlinear spectrum evolution. *Journal of Physical Oceanography*, 31(4), 861-873. https://doi.org/10.1175/1520-0485(2001)031(0861:EOWWPF)2.0.CO;2

* Longuet-Higgins, M. S., & Stewart, R. W. (1962). Radiation stresses and mass transport in surface gravity waves with application to 'surf beats. *Journal of Fluid Mechanics*, 13(4), 481-504. https://doi.org/10.1017/s0022112062000877
Moreau, T., Tran, N., Aublanc, J., Tison, C., Gac, S. L., & Boy, F. (2018). Impact of long ocean waves on wave height retrieval from sar altimetry data. *Advances in Space Research*, 62(6), 1434–1444. https://doi.org/10.1016/j.asr.2018.06.004

Morrow, R., & Le Traon, P.-Y. (2012). Recent advances in observing mesoscale ocean dynamics with satellite altimetry. *Advances in Space Research*, 50(8), 1062–1076. https://doi.org/10.1016/j.asr.2011.09.033

Passaro, M., & Algorithm Development Team. (2021). Algorithm theoretical basis document (atbd), sea state climate change initiative (Tech. Rep.). European Space Agency. Retrieved from https://climate.esa.int/media/documents/Sea_StateCCI_ATBD_v3.0-signed.pdf

Passaro, M., Hemer, M. A., Quartly, G. D., Schwatke, C., Dettmering, D., & Seitz, F. (2021). Global coastal attenuation of wind-waves observed with radar altimetry. *Nature Communications*, 12(1), 3812. https://doi.org/10.1038/s41467-021-23982-4

Quartly, G. D., Smith, W. H. F., & Passaro, M. (2019). Removing intra-1-Hz covariant error to improve altimetric profiles of σ₀ and sea surface height. *IEEE Transactions on Geoscience and Remote Sensing*, 57(6), 3741–3752. https://doi.org/10.1109/TGRS.2018.2886998

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

Quilfen, Y., Yurovskaya, M., Chapron, B., & Ardhuin, F. (2018). Storm waves sharpening in the agulhas current: Satellite observations and modeling. *Remote Sensing of Environment*, 216, 561–571. https://doi.org/10.1016/j.rse.2018.07.020

Rice, S. (1944). Mathematical analysis of random noise. In N. Wax (Ed.), *Noise and stochastic processes* (pp. 133–294). Dover Publications Inc. (published 1954).

Rieu, P., Moreau, T., Cadier, E., Raynal, M., Clerc, S., Donlon, C., et al. (2021). Exploiting the sentinel-3 tandem phase dataset and azimuth oversampling to better characterize the sensitivity of sar altimeter sea surface height to long ocean waves. *Advances in Space Research*, 67(1), 253–265. https://doi.org/10.1016/j.asr.2020.09.037

Rodriguez, E. (1988). Altimetry for non-Gaussian oceans: Height biases and estimation of parameters. *Journal of Geophysical Research*, 93(C11), 14107–14120. https://doi.org/10.1029/JC093iC11p14107

Sandwell, D. T., & Smith, W. H. F. (2005). Retracking ers-1 altimeter waveforms for optimal gravity field recovery. *Geophysical Journal International*, 163(1), 79–89. https://doi.org/10.1111/j.1365-246X.2005.02724.x

Schlembach, F., Passaro, M., Quartly, G. D., Kurekin, A., Nencioli, F., Dodet, G., et al. (2020). Round robin assessment of radar altimeter low resolution mode and delay-Doppler retracking algorithms for significant wave height. *Remote Sensing*, 12(8), 1254. https://doi.org/10.3390/rs12081254

Srokosz, M. A. (1986). On the joint distribution of surface elevation and slopes for a non linear random sea, with an application to radar altimetry. *Journal of Geophysical Research*, 91(C1), 995–1006. https://doi.org/10.1029/jc091ic01p00995

Tayfun, A., & Lo, J.-M. (1989). Wave envelope and related spectra. *J. Of Waterway, Port Coast. Ocean Engineering*, 115(4), 515–533. https://doi.org/10.1061/(ASCE)0733-950X(1989)115:4(515)

Tourain, C., Piras, F., Ollivier, A., Hauser, D., Poisson, J. C., Boy, F., et al. (2021). Benefits of the adaptive algorithm for retracking altimeter nadir echoes: Results from simulations and CFOSAT/SWIM observations. *IEEE Transactions on Geoscience and Remote Sensing*, 59(12), 9927–9940. https://doi.org/10.1109/TGRS.2021.3064236

Vignudelli, S., Birol, F., Benveniste, J., Fu, L., Picot, N., Raynal, M., & Roinard, H. (2018). Satellite altimetry measurements of sea level in the coastal zone. *Surveys in Geophysics*, 40(6), 1319–1349. https://doi.org/10.1007/s10712-019-09569-1

Young, I. R. (1986). Probability distribution of spectral integrals. *J. Of Waterway, Port Coast Ocean Engineering*, 112(2), 338–341. https://doi.org/10.1061/(asce)0733-950x(1986)112:2(338)

Zaron, E. D., & DeCarvalho, R. (2016). Identification and reduction of retracker-related noise in altimeter-derived sea surface height measurements. *Journal of Atmospheric and Oceanic Technology*, 33(6), 201–210. https://doi.org/10.1175/JTECH-D-15-0164.1