the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A method for constructing directional surface wave spectra from ICESat2 altimetry
Christopher Horvat
Sea ice is important for Earth's energy budget as it influences surface albedo and air–sea fluxes in polar regions. On its margins, waves heavily impact sea ice. Routine and repeat observations of waves in sea ice are currently lacking, and therefore a comprehensive understanding of how waves interact with sea ice and are attenuated by it is elusive.
In this paper, we develop methods to separate the twodimensional (2D) surface wave spectra from seaice height observations made by the ICESat2 (IS2) laser altimeter, a polarorbiting satellite. A combination of a linear inverse method, called generalized Fourier transform (GFT), to estimate the wave spectra along each beam and a Metropolis–Hastings (MH) algorithm to estimate the dominant wave's incident angle was developed. It allows us to estimate the 2D wave signal and its uncertainty from the highdensity, unstructured ATL03 ICESat2 photon retrievals. The GFT is applied to rebinned photon retrievals on 25 km segments for all six beams and outperforms a discrete Fourier transform (DFT) in accuracy while having fewer constraints on the data structure.
The MH algorithm infers wave direction from beam pairs every 25 km using coherent crests of the most energetic waves. Assuming a dominant incident angle, both methods together allow a decomposition into 2D surface wave spectra with the advantage that the residual surface heights can potentially be attributed to other seaice properties. The combined GFT–MH method shows promise in routinely isolating waves propagating through sea ice in ICESat2 data. We demonstrate its ability on a set of example ICESat2 tracks, suggesting a detailed comparison against in situ data is necessary to understand the quality of retrieved spectra.
 Article
(4619 KB)  Fulltext XML

Supplement
(1759 KB)  BibTeX
 EndNote
Sea ice covers up to 9 % of the world's oceans and plays an important role in the energy balance of Earth's climate. Even though sea ice damps ocean surface waves (Squire, 2007), broad regions along the periphery of the seaicecovered ocean are continually under the influence of surface waves (Rapley, 1984; Horvat et al., 2020; Thomson, 2022; Horvat, 2022). These regions are collectively referred to as the Marginal Ice Zone (MIZ). In the MIZ, waves influence sea ice's thermodynamic and dynamic properties and impact the coupled exchange between atmosphere and ocean. Currently, we do not have reliable global observations of waves in sea ice and hence are unable to sufficiently understand air–sea exchange and wave propagation in the MIZ. This paper describes how ICESat2 (IS2) altimeter observations can be used to record wave spectra in the MIZ and how to infer additional seaice properties for building parametrizations of wave attenuation in sea ice.
Models of wave propagation in sea ice typically evolve the ocean surface wave spectrum, ${\stackrel{\mathrm{\u0303}}{S}}_{h}\left(k\right)$ (m^{2} k^{−1}, k is the wavenumber), which is attenuated when it comes into contact with sea ice. There has been significant debate over the functional form and dependencies of this attenuation (Squire, 2018; Thomson et al., 2021). Yet as it controls how deep waves reach into sea ice, it is vital for modeling MIZ variability and coupled feedbacks in the polar seas.
Constraining iceinduced wave attenuation is challenging because wave observations in ice are difficult to make at scale. A majority of observations of waves in ice are carried out using ships or arrays of floating buoys deployed by ships (or by helicopters from ships; see, for example, Thomson, 2022, and references therein). While such observations provide hightemporalfrequency observations of wave spectra, they only cover a limited geographic domain and are limited by the seaice types and conditions at the original buoy locations. Recently satellite remote sensing technologies have shown promise for describing wave spectra in seaice regions. Synthetic aperture radar (SAR) imagery is capable of observing wave crests as they move into the MIZ, and the twodimensional wave spectrum can be constructed in good agreement with spectra observed in situ if the sea ice is not rough (Stopa et al., 2018; Ardhuin et al., 2017). However, SAR alone cannot observe continuous spectra as they propagate into the sea ice.
The IS2 altimeter has the potential to greatly increase the quantity of available observations of wave–ice interactions, either alone or in combination with other remote sensing instruments (Collard et al., 2022). IS2 carries a single measurement tool, ATLAS, a sixbeam laser oriented in three weak–strong beam pairs (Fig. 1a, colored lines) offset at nearuniform 3 km on the ground, with a weak–strong beampair lateral offset of about 90 m. ATLAS measures the return time of individual photons to infer the height of the ice–ocean surface. Typical alongtrack photon spacings can be centimeters or smaller, and so IS2 is capable of directly sampling ocean surface waves, particularly over reflective sea ice.
Recent studies have examined waves in sea ice using IS2, basing their results on a higherorder seaice height product derived from photon retrievals (known as ATL07). Horvat et al. (2020) identified the capability of IS2 to retrieve ocean waves by examining a storm in the Barents Sea in 2019 and used a simple threshold to establish where and when waves were observed in the sea ice to produce global maps of the MIZ. In Collard et al. (2022), IS2 retrievals during this Barents Sea storm were shown to compare well with model and SARbased observation data. Brouwer et al. (2022) selected a series of Southern Hemisphere IS2 retrievals, analyzing wave attenuation using direct spectral transform methods. Both found that areas affected by waves were common in both hemispheres, with repeated measurements of waves hundreds of kilometers into the seaice zone, particularly in the Southern Hemisphere.
Three challenges limit the direct comparison of IS2derived wave spectra to observations and models. First, waves propagate at an angle θ relative to the alongtrack direction of the satellite (Fig. 1b), and observable wavelengths λ are aliased by an unknown factor cos (θ)^{−1} (Rapley, 1984; Horvat et al., 2019; Yu et al., 2021). Second, observed surface height variability is a convolution of the dynamic ocean topography, seaice topography, surface waves, and noise. The surface wave signal can then be successfully reconstructed if these other signals are on a different scale, like dynamic ocean topography, or not periodic, like the seaice topography and noise. Third, the fractured nature of sea ice, the influence of clouds, and changing surface albedo cause gaps or irregularities in IS2 photon retrieval rates, creating a highdensity but irregular observation. The method must be applicable to irregular data without generating spurious sources of variance, i.e., artificial wave energy. The above factors complicate direct assessments of spectra and their attenuation in sea ice.
Here we demonstrate a method for producing anglecorrected, twodimensional (2D) wave spectra in sea ice using photon height data from IS2. We partition surface height variability into waves and seaice or noiserelated components. It permits a direct assessment of the most significant wave energy along each track to record wave attenuation and evaluate numerical attenuation schemes. We show this partitioning allows for significantly improved seaice height estimates in the MIZ and may also allow for expanding existing higherlevel IS2 products to broader icecovered regions.
In this paper, we demonstrate this method on four example cases, Tracks 1 to 4 (details in Table S1 in the Supplement. Their granule, i.e., their identification number, is also given in each figure). We describe the preprocessing of IS2 alongtrack photon heights in Sect. 2.1 and develop a harmonic fitting procedure applied to individual IS2 beams in Sect. 2.2 (the generalized Fourier transform, GFT, method). Subsequently, we develop a multibeam, Monte Carlo method for biascorrecting alongtrack wavelengths in Sect. 3 (the Metropolis–Hastings, MH, method), which enables us to provide twodimensional wave spectra derived from alongtrack data in Sect. 3.2, and a decomposition of photon variance in Sect. 4. In Sect. 5, we discuss the limitations and assumptions of the proposed methods and conclude in Sect. 6 how they can be used to develop improved models of attenuation of waves in sea ice.
The primary aim of this analysis is to assess surface height variability in the MIZ. Hence we want to use the highest data resolution we can handle, though we are agnostic about the classification of photon returns. That is the L2level product ATL03 from Neumann et al. (2021). For comparison, we show the photon cloud data from ATL03 and the surface heights and type classifications from the higherlevel ATL07/10 product in Fig. 2 as dark blue, light blue, and orange dots (Kwok et al., 2021). By requiring 150 consecutive photons to identify a seaice segment, the ATL07 product accounts for most of the height variability from the ATL03 product. Yet it misses retrievals in the MIZ (Fig. S2, white and gray area) and within the sea ice (Fig. 2). For better resolution, the following analysis is based on the photon cloud data from ATL03.
2.1 Data preprocessing
Linearly inverting photon data require exact along and acrosstrack information about photon positions. Alongtrack photon positions are first rereferenced to the most equatorward position on the nominal ATLAS ground track (Fig. 1b, black line). The most equatorward position is evaluated from the ATL07 (Kwok et al., 2021) product and set to the beginning of the first 100 km of alongtrack data where there is an average of at least 0.02 photons per meter (defined as X=0 throughout the paper). This threshold and rereferencing are used to exclude large areas of nearly no data in the transition zone between the open ocean and MIZ. All tracks are then followed in a poleward direction, until the variance of any of the six beams exceeds a factor of 10 times the variance of the first 15 % on the equatorward side of the track (Fig. S2, dashed black lines). This avoids including impacts from coastal or land ice around the Antarctic coast. The redefined alongtrack direction x and an acrosstrack direction y are then used as the coordinate system throughout the analysis (Figs. 1b, S1 and S2).
After removing the cumulative surface height correction (dem_h
taken from the ATL07/10 dataset; Kwok et al., 2021), we bin photon measurements into 20 m stencils that overlap by 50 %. This yields a 10 m alongtrack resolution (Fig. 2 green line). Note this differs from the procedure used to form seaice surface heights in the ATL07 product, which averages height data for each set of 150 photons alongtrack. ATL07 has a constant photon count, with the tradeoff of irregularly spaced segments of varying length in the MIZ, while our approach provides more regularly spaced data with the tradeoff of having a variable photon count in each stencil and potentially including retrievals of water near sea ice. Stencils with fewer than five photons are excluded, which also leads to data gaps corresponding to no or very low photon retrievals due to seaice leads, open water, clouds, or other noise (Fig. S1 in the Supplement). These data gaps also lead to an irregularity in the data, but here each stencil mean represents the same area and hence better captures the wave phase.
A mean photon height h_{c}(x) is calculated in each 20 m stencil as the mean of the photons weighted by their inverse distance to the stencil center, using a Gaussian weighting function with a standard deviation of 10 m. We tested other data reduction methods, like using the median or mode of the stencil, finding results insensitive to the choice of the binning method (Fig. S2). The same 20 m stencil also provides an uncertainty estimate σ_{h}(x) (proportional to Fig. 2 blue area) representing the varying photon density. This uncertainty is used to define the data prior (Appendix A1) for the harmonic inversion in Sect. 2.2.
The resampled mean photon height data are used to calculate a series of alongtrack surface slopes (Fig. 2 thin blue line) by taking the alongtrack derivative and applying a spikeremoving algorithm. Using the alongtrack surface slope data focuses the time series analysis on local photon height changes rather than the magnitude of the total surface height field. The spike removal reduces peaks in the slope data, which can be from seaice height changes or especially those resulting from ice–ocean transitions. Because the slope field approximates the derivative of the height field, the spectrum of slopes ${\widehat{S}}_{c}\left(k\right)$ is readily connected to the spectrum of surface heights ${\stackrel{\mathrm{\u0303}}{S}}_{h}\left(k\right)$, with ${k}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u0303}}{S}}_{h}\left(k\right)={\stackrel{\mathrm{\u0303}}{S}}_{c}\left(k\right)$. The surface height field can therefore be directly reconstructed from the slope spectrum, as we show below. The generalized Fourier transform (Sect. 2.2) and directional estimates (Sect. 3.1) are then applied on 25 km long segments of these surface slopes, with uncertainty estimates (Fig. 2). The 25 km segments also overlap by 50 %, providing an alongtrack spectral estimate every 12.5 km.
2.2 Generalized Fourier transform (GFT)
We estimate alongtrack wave spectra using a generalized Fourier transform (GFT). The GFT is a harmonic fit of sin and cos pair coefficients, which together determine amplitude and phase at each wavenumber. The model complexity is defined by the number of resolved wave numbers, and its success depends on prior (assumed) knowledge about the data's uncertainty and model structure.
We use a GFT to overcome several disadvantages appearing when implementing a standard discrete Fourier transform (DFT) to unstructured data. While the DFT is a fast varianceconserving algorithm, it requires periodic, continuous, and equally spaced data. The DFT implies that frequency bands are harmonics over a domain or segment length L, an arbitrary limitation on the resolved frequencies. To make segments periodic, often tapering or windowing is applied to the segmented data. In addition to the data's nonperiodicity, the common presence of data gaps in IS2 retrievals requires extrapolation or gapfilling to create continuous, equally spaced data suitable for a DFT. These both lead to commonly known problems of the DFT, like energy leaks or compensation in spectral space. The data handling needed for DFT can erode the signal substantially, especially in the MIZ (Fig. 3b and c, gray and green lines).
The GFT method outlined below works on any grid, incorporates data uncertainty, and does not require periodicity. The GFT can be customized to the frequencies of interest with the additional benefit of providing a standard error in real and wavenumber space.
2.3 Harmonic wave inversion
We follow Wunsch (1996), Menke (2018), and Kachelein et al. (2022), using generalized least squares to derive spectral amplitudes in each 25 km segment X_{i} (Fig. 4). Slope data in each section X_{i} are a series of unevenly spaced meanzero data points and are expressed as a column vector $\mathit{b}={\partial}_{x}\mathit{z}$ of length N_{i}, where z is the height data. These data are then modeled as the sum of sinusoidals with wave numbers in the range of swell and wind waves:
where H is an N_{i}×2M regressor matrix of basis functions, p is the model parameter vector of length 2M, and r is the vector of the residual time series of length N_{i}. The columns of H are sines and cosines of prescribed wave numbers, ${k}_{m}^{\prime}=\mathrm{2}\mathit{\pi}/{\mathit{\lambda}}_{m}$ for wavelengths λ_{m} indexed by $m=\mathrm{1},\mathrm{2},\mathrm{\dots}M$. We use prime notation here and throughout the paper to indicate observed wave variables along the direction of each beam and unprimed notation for variables in the direction of the traveling wave (see Fig. 1c). The problem can then be written as
with model parameters
at positions
To find the most probable b given a set of model parameters p, we need to estimate the autocovariance matrices of the residual R=〈rr^{T}〉, i.e., the error of the data, and the autocovariance matrix of the model P=〈pp^{T}〉, where 〈⋅〉 is the expected value. Then the most probable estimate of the data b can be found by estimating the maximum of the posterior probability distribution P(pb). Using Bayes' theorem (Kachelein et al., 2022) or, alternatively, the matrix inversion lemma (Wunsch, 1996), given the data b, the most likely estimate of the model parameters $\widehat{\mathit{p}}$ is found as
Once the model parameters $\widehat{\mathit{p}}$ are estimated, the model can be expressed in real space using
(Fig. 3a, green lines) or as a power spectrum as ${\widehat{\stackrel{\mathrm{\u0303}}{\mathit{S}}}}_{\text{GFT}}$ (Fig. 3b and c green lines; see Sect. S1 in the Supplement for derivation). Note that ${\widehat{\stackrel{\mathrm{\u0303}}{\mathbf{S}}}}_{\text{GFT}}$ is substantially different from that of the DFT of the same data (Fig. 3b and c, gray and green lines, respectively). Gaps in the data, as well as the DFT's requirement for the data to be periodic, create artificial power in the swell's wavenumber range that leads to misleading results.
To estimate the model error, the posterior autocovariance of the difference between estimated and true model parameters is defined as the inverse Hessian:
In practice, this is calculated from the “kernel” matrix H and the (assumed) Gaussian distributed data and model priors R and P. The trace tr(⋅) of the inverse Hessian is then used to estimate the error of the model parameters:
(Fig. 3d, shown in the same units as the power spectra). The error of the fit to the modeled data is also related to the inverse Hessian:
where j is a unit vector of length 2M.
The harmonic inversion of this segment of Track 2 shows how wave spectra can be calculated from strong (gt2r) and weak (gt2l) beams, even if the data have gaps (Fig. 3a). Both beams' spectra are similar in most parts and show a maximum at k=0.03 (Fig. 3b and c). However, the strong beam (gt2r) does show a second local maximum at about k=0.06. Note that the DFT of the same track, and with tapered data, results in a different probability density function (PDF) than the harmonic inversion because of the data gaps.
2.3.1 Choice of model resolution and degrees of freedom
The quality of the GFT model depends on the degrees of freedom, the model prior P, and data priors R (Eq. 7). While the number of model parameters, 2M, remains fixed throughout the analysis, the number of data points in segment i, N_{i}, is variable and may be controlled by the segment length L, which is in our case 25 km, with a 12.5 km overlap (Fig. 4).
The number of degrees of freedom is 2M−N_{i} and depends on the number of data points in each 25 km segment. A segment with no data gaps and a 10 m resolution (Sect. 2.1) contains N_{i}=2500 data points. With 2M=1740 model parameters, this is an overdetermined problem. However, frequent data gaps reduce the data points per segment, which may result in an underdetermined problem with more model parameters than data (2M>N_{i}). The result is then a larger residual r and larger uncertainty estimate ${\widehat{\mathit{p}}}_{\text{err}}$ (Eqs. 8 and 9). Even in cases where Eq. (1) is underdetermined, we are confident in our wave spectrum estimation within a given error because P contains prior knowledge about the shape of the solution, i.e., the shape of typical surface wave spectra (Sect. 2.4). Most of the segments of the four example tracks in this study have close to 2500 data points and are overdetermined (Fig. S3 in the Supplement). Only Tracks 1 and 2 (Fig. 3) are underdetermined close to the edge of the ice cover (Fig. S3a and b).
The choice of segment length also determines the smallest resolvable wavenumber. For example, a segment length of 25 km resolves a wave with a 20 s period at an incident angle of $\pm \mathrm{75}{}^{\circ}$ about 10 times. We set the lowest observed wavenumber to be resolved to ${k}_{\mathrm{1}}^{\prime}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{3}}$ rad m^{−1}, which corresponds to a maximum observed wavelength of 2500 m (Sect. 3.1). The highest wavenumber is chosen as ${k}_{m}^{\prime}=\mathrm{0.11}$ rad m^{−1}, a typical period of wind waves of about 6 s. Using evenly spaced wavenumbers with $dk=\mathrm{1.25}\times {\mathrm{10}}^{\mathrm{4}}$ results in M=869 wavenumbers.
2.4 Iterative inversion along each beam
The GFT solution $\widehat{\mathit{b}}$ depends on prior assumptions about the wave spectrum: the model prior P. Since the GFT is iteratively applied along each beam, results from a previous segment inform the subsequent segment as illustrated in Fig. 4. Here we describe how a successive application of the GFT along the IS2 beam can lead to an efficient solution assuming that the wave's spectral shape only slowly varies between the segments.
The iterative solution along each beam is initialized at the most oceanward edge of the data with a prior P_{init} that follows a common shape of a narrow banded surface wave field: the Pierson–Moskowitz (PM) spectral slope function (based on Pierson and Moskowitz, 1964). For this segment alone, the PM function is fit to a DFT of the data, and, for the DFT only, any locations with missing data are defined with a slope of zero. This gapfilling creates artificial ringing in the DFT but is sufficient to estimate the spectrum's peak wave number and energy. The PM spectrum has, in its simplest form, only two free parameters, the peak frequency and spectral amplitude, which are fit to the DFT power spectra via an objective function that is regularized by the observed spectral peak of the smoothed data (Appendix A; Hell et al., 2019).
The initial inversion of the most equatorward segment X_{0} is performed using P_{0}=P_{init} in Eq. (7), leading to model parameters ${\widehat{\mathit{p}}}_{\mathbf{0}}$. For this first segment, a second inversion is applied on the same data, using an updated prior that is a smoothed version of ${\widehat{\mathit{p}}}_{\mathbf{0}}$ (Fig. 4, left). The smoothing uses a Lanczos running smoother in wavenumber with a stencil width of ±0.19 2π m^{−1} or 150 data points. Inversions of the successive segments ${X}_{\mathrm{1}},{X}_{\mathrm{2}},{X}_{\mathrm{3}},\mathrm{\dots}$ are then performed once, with the prior P_{i} a smoothed version of ${\widehat{\mathit{p}}}_{i\mathrm{1}}$. If missing data do not allow for a successful inversion of a segment, the algorithm is reinitiated as done to obtain P_{0} and ${\widehat{\mathit{p}}}_{\mathbf{0}}$.
This twostage inversion for segments with no preceding alongtrack segment ensures that ordinary wave spectra will be identified at the margin while still having the flexibility to allow for more complex wave signals. The effect of the PM prior and details about the derivation are shown in Appendix A.
2.5 Tracking of wave energy through sea ice
The GFT is applied to each 25 km segment with more than 250 data points, leading to wave spectral estimates along each beam. In Fig. 5a–f we show wavenumber spectra for each segment and for each of the beams of Track 3 in the Southern Hemisphere on 2 May 2019. The persegment crossbeam mean (Fig. 5h) and mean spectral error (Fig. 5i) are derived by weighting each segment by its photon density before taking the mean. We define this weighted mean and error as our best estimate of the alongtrack spectral evolution of wave energy.
The example track shows an attenuating swell signal starting at X=0 km (X is the distance from the ice edge as defined in Sect. 2.1) and a second waveenergy maximum with shorter wavelength (${k}^{\prime}\approx \mathrm{100}$ m) at about X=150 km both in the best estimate and in individual beams (Fig. 5a–f and h). (These wave events are further analyzed in Sect. 4.) The corresponding errors are often larger where the wave signals are larger (Fig. 5i). Instances with a low photon density and more frequent data gaps may fail to invert for the wave signal, resulting in a spectrogram that may not follow expected spectral shapes. These can be identified through their substantially larger error (Fig. S4g–i in the Supplement).
The estimated wave numbers are the observed alongtrack wave numbers k^{′}, which are different than the true wave number length $\left\mathit{k}\right$ along the incident wave vector k (see Fig. 1). To estimate a wavenumber spectrum along the dominant propagation direction rather than the direction it is observed, we outline a method to correct this bias in the following section.
3.1 Metropolitan estimates of the incident angle
The observed wave spectra are distorted by a misalignment between the wave's incident angle and the beam's direction (Fig. 1c). While the IS2 track orientations are wellknown, the surface waves can originate from any direction, and the angle between these two directions is θ. If θ is known, the observed wave number k^{′} or wavelength λ^{′} along the beam can be corrected using $k={k}^{\prime}\mathrm{cos}(\mathit{\theta}{)}^{\mathrm{1}}$. The same geometrical distortion will also affect estimates of the attenuation rate between X_{i} positions along the track (Fig. 1b) because the wave energy attenuates along their dominant propagation direction and not along the direction they are observed by the satellite.
We use the phase lag between weak–strong beam pairs to estimate θ from the photon data. This requires that wave crests observed in one beam are also observed in another. Using the phase lag to measure the incident angle has several limitations that have to be taken into account when designing an optimization method:

The phase lag between beam pairs can only usefully be calculated for angles that are not too oblique (Fig. S5 in the Supplement and Yu et al., 2021) and that have high enough photon densities in both beams. The angle limits in which the phase lag can be resolved depend on the chosen wavelength and the distance between the beams. Since both wavelength and distance change for each segment, we here limited the analysis to angles of $\pm \mathrm{75}{}^{\circ}$.

Acrossbeam phase lag relies on accurately measuring the distance between the beams. Uncertainty in the intrabeam distance d can further obscure angle estimates based on the phase lag (Fig. 1b). Since d varies along the track (Fig. S6 in the Supplement), using the nominal distance rather than the observed value will also bias angle estimates.

The complex generation and propagation history of waves (Kitaigorodskii, 1962; Villas Bôas and Young, 2020; Marechal and Ardhuin, 2021; Hell et al., 2021a) leads to a dynamic distortion in the incident angle. While a monochromatic plane wave would be coherent across beam pairs, estimating its direction is limited by the periodicity of the waves and observational noise. In reality, the incident swell wave energy at any given time is contained in several wavenumbers although concentrated in a narrowbanded 2D spectrum (LonguetHiggins and Deacon, 1957). A narrowbanded 2D swell field leads to wave groups in real space and limits the observable phase lag between strong and weak beam pairs. While a narrower 2D spectrum results in spatially larger coherent wave crests, a broader spectrum spans more randomphase waves, diminishing the observable phase lag between wave crests (Fig. S7 in the Supplement).
With IS2, the incident wave energy along k is observed along the beam direction for wavenumbers ${k}^{\prime}=\left\mathit{k}\right\mathrm{cos}{\mathit{\theta}}^{\mathrm{1}}$. In the case of IS2 beam pairs, we know neither θ nor the bandwidth of the incident spectrum (LonguetHiggins, 1984), and both factors limit the possible angle inversion based on beam pairs alone (Fig. S7). Despite these limitations, we describe in the following how the incidence angle θ can still be retrieved within these limitations.
As explained in Sect. 2.2, the surface wave field can be interpreted as the superposition of monochromatic plane waves. For a narrowbanded swell spectrum, the majority of the wave energy is contained in a few wavenumbers, and hence a superposition of these most energetic monochromatic waves explains most of the surface slope variability. In the following, we optimize the incident angle and phase of the most energetic monochromatic wave numbers using a Metropolis–Hastings (MH) algorithm (ForemanMackey et al., 2013). We accumulate the marginal distributions of possible incident angles across the dominant wave numbers and beam pairs, which results in the best guess of the possible incident angle. This approach leads to directional wave information similar to the maximum entropy method used in wave buoys (Lygre and Krogstad, 1986).
We focus on the 25 most energetic wavenumbers of each beam pair and segment X_{i} based on the GFT result (Sect. 2.2). To identify these wavenumbers, the beampair mean wave power is smoothed using a threewavenumber running mean to select possible wave numbers within the distorted narrowbanded spectrum (similar to the thick green lines in Fig. 3b and c). For each of these n=25 observed wavenumbers ${k}_{n}^{\prime}$, we define a monochromatic model,
in the local reference system of the segment centered around X_{i} and y=0 such that the local alongtrack coordinate is $\mathit{\eta}=\mathit{x}{X}_{i}$ and the acrosstrack coordinate is ν=y with the observable acrosstrack wavenumber ${l}^{\prime}={k}^{\prime}\mathrm{tan}\left(\mathit{\theta}\right)$ and the phase ϕ.
The monochromatic model is then used to define the objective function Φ_{n} for each wavenumber:
where b is the normalized slope data of the beam pair, β_{θ} is a hyperparameter which controls the regularization P_{θ,n} of the incidence angle θ for the nth wavenumber, and ${P}_{\mathit{\theta},n}(\mathit{\theta},k)$ is a prior estimate that we describe in Sect. 3.1.2.
The log probability of the objective function Eq. (11) is sampled for each beam pair, selected wavenumber, and alongtrack position X_{i}. To derive independent estimates of the incident angle for each nth wavenumber we use an MH scheme (Markov chain Monte Carlo, MCMC; ForemanMackey et al., 2013) by first initializing equally spaced samples of the objective function over the domain $\mathit{\theta}=[\mathrm{0.42}\mathit{\pi},\mathrm{0.42}\mathit{\pi}]$ and $\mathit{\varphi}=[\mathrm{0},\mathrm{2}\mathit{\pi}$) and advancing the ensemble of samples (ensemble of walkers) using MCMC. The MCMC method will quickly cluster walkers in the areas of lowcost or small objective function (Fig. 6b, black dots). A high density of walker positions is then interpreted as a high likelihood of an incident angle and phase for the chosen wavenumber (Fig. 6b, black dots).
We derive a sample of the joint phase and angle distribution by advancing the walkers 300 iterations, and only the last 270 iterations for each walker are used to establish the joint histogram D(θ,ϕ). The joint histogram D is then marginalized for the incident angle θ and normalized to a probability distribution function (Fig. 6c). This procedure is repeated for each selected wavenumber k_{n} and for each (available) beam pair per segment X_{i} (Fig. S8a–c in the Supplement). The best incident angle PDF ${\mathit{\theta}}_{\text{PDF}}(X,k,\text{beam pair})$ can then be derived using the weighted average across wavenumber, beam pair, or both.
An example of the resulting beam and wavenumberaveraged PDF is shown in Fig. 7b for Track 3 at X_{i}=87. Here, the individual PDFs are weighted by the mean power of the respective wavenumber and the number of data points per segment pair. The most likely incident angle is at $\mathrm{37}{}^{\circ}$, while two other angles −63 and 0^{∘} also show high likelihood. Marginal PDFs with multiple maxima are a typical result for this method and appear in many other tested sections and tracks (not shown). They come from different maxima in the joint PDFs of different wavenumbers. If one trusts the angle estimate of a single wavenumber, this result can be interpreted as a wave field with waves from multiple directions. The alternative – likely better – perspective is that the marginal PDF of a single wavenumber is not a robust estimate of the incident angle, and hence the PDF in Fig. 7b helps estimate the uncertainty of the method.
3.1.1 Robustness of the marginal PDFs
The limitations in retrieving the incident angle (Sect. 3.1) lead to a low signaltonoise regime and demand a careful evaluation of the method for sampling the objective function Φ_{n}. While larger samples may ensure convergence of the distribution estimate, a large sampler of each wavenumber, beam pair, and segment may not be necessary or computationally affordable. We decided for a systematic undersampling of each realization of the marginal PDFs $\mathit{\theta}({X}_{i},k,\text{beampair})$ and, in a second step, made a supersample from the marginal PDFs across beam pairs and wavenumber if needed.
When combining the systematic undersampled PDFs, their supersample will still provide a good estimate of the mean incident angle and its standard deviation. Each realization of the joint θ–ϕ distribution requires 6750 function evaluations for 270 iterations per walker. The walker's autocorrelation is about 20 to 30 iterations, which implies that each joint distribution may not be well established (the effective degrees of freedom per walker are about 9 to 14). Hence the marginalization of each joint distribution may misrepresent the angle uncertainty (i.e., a distribution of the walker's PDF that is too wide). To reduce uncertainty, we take a supersample of the marginal PDFs, by averaging across wavenumbers and/or the three beam pairs. The supersampling results in a statistically robust result with 5×10^{5} function evaluations per segment X_{i}, which is about 3.5–5.5×10^{3} effective degrees of freedom per segment. Longer Markov chains, i.e., more iterations, may result in a better sampling of the individual fit but may not affect the overall result since they sample from generally smooth objective functions in a low signaltonoise regime. However, in cases where a directional estimate per wavenumber or beam group is needed, the MCMC iteration length can be adjusted.
3.1.2 Constraining direction estimates with other data products
A sampling of the objective function Eq. (11) as described in Sect. 3.1 results in a joint distribution of most likely incident angles and phases per sampled wavenumber k_{n}. This joint distribution may have multiple equally likely maxima, i.e., multiple likely incident angles due to the periodicity of the wave (2π ambiguity). As illustrated in Fig. 6d (shading) this can lead to (a) maxima for positive and negative incident angles and (b) multiple maxima on both sides. To break the symmetry in the marginalized PDF of incident angles (Fig. 6e) we define a prior P_{θ}(θ,k) in the objective function using ridge regression (Appendix B). The effect of the prior on the joint and marginalized distribution is shown by comparing Fig. 6b and c with d and e. Here we inform the prior with WAVEWATCH III (WW3) global hindcast wave partitions (WW3DG, 2009, using the Integrated Ocean Waves for Geophysical and other Applications, IOWAGA, hindcast). WW3 must be treated with caution due to wind observational biases in the Southern Ocean (Belmonte Rivas and Stoffelen, 2019; Hell et al., 2020). This wave hindcast is currently the only readily available dataset for this global purpose, and priors from observational datasets would improve the quality of these data and the overall wave inversion.
The level of certainty in the WW3 prior is expressed in the hyperparameter β_{θ}, and the performance of the MCMC sampling is sensitive to its value (Eq. 11). Since validation of the WW3 prior is limited, we set β_{θ}=2. Its effect on the objective function can be seen by comparing the shading in Fig. 6b and d. The choice of β_{θ}=2 leads to the desired result in breaking the directional ambiguity while not fully determining the incident angle distribution (Fig. 7a). We tested other values of β_{θ} but found empirically that higher values tend to overfit to the prior, and lower values do not break the ambiguity well.
This method is limited to angles of about $\pm \mathrm{75}{}^{\circ}$ deviation from the nominal track direction. A more oblique, i.e., steeper, incidence angle can not be captured by this method because a steeper angle requires more coherence between wave crests. The coherence of a single wave crest is, however, limited by the curvature of the wave spectrum and not wellknown (Fig. S7). In addition, the model has a 180^{∘} ambiguity such that waves coming from the Equator side of the track (as assumed) or waves coming from the pole side (less likely) can result in the same phase lag and hence in the same incident angle, even though they come from the opposite direction.
3.2 Twodimensional spectra in alongtrack data
With the spectral and angle estimates (Sects. 2.2 and 3.1), we now can describe waves observed alongtrack in terms of their twodimensional wavenumber spectra (Fig. 8). The estimated wavenumber amplitudes $\widehat{\mathit{b}}$ (Eq. 6) are corrected by cos (θ)^{−1} using the most likely incident angle (Sect. 3.1 and Fig. 7b) resulting in the corrected wavenumber spectrogram (Fig. 8a). We use the most likely angle along the track, although the above analysis can provide angle distributions for each segment X_{i} and wavenumber k (Sect. 3.1.2 and Fig. 8b).
The corrected power spectrum and directional distribution (Fig. 8a and b) can be expressed as directional surface wave spectra every 12.5 km in the MIZ, similar to conventional surface wave buoys (Fig. 8c–e). This permits tracking the attenuation of wave energy per frequency in MIZ. In the case of Track 3 (granule 05160312), for example, we see a wave event coming from about 45^{∘} to the right of the ground track that mostly attenuates in the first 75 km from the seaice edge, while the overall attenuation rate is similar between the six beams (Fig. 5a–f). One could identify a migration of the peak wavelength from about 275 m to about 300 m within 12.5 km (Fig. 8d and e, similar to Alberello et al., 2022); we leave this analysis of the attenuation to future work.
Past the primary wave event, a second signal further into the ice with energy on scales shorter than 200 m extends from X>75 km for about 100 km or more (Fig. 8a and b) in the ice. This signal is at shorter wavelengths than the identified cutoff frequencies for the event at the track beginning. Without further information about the ice conditions, we suggest that this shortwavelength energy deeper in the sea ice is due to seaice variability itself rather than due to waves.
The estimate of the wave signal from Sect. 2.3 can be used to decompose wave and ice surface variability. Each photon retrieval is a superposition of ocean waves and seaice signals like surface roughness, floe size, and freeboard height. A decomposition of the surface variability between waves and ice can rely on the coherence of acrossbeam wave statistics, a common noise level in wavenumber space, and an approximate scale separation of the dominant wave energy and the sea ice. For the purpose of decomposing the data we define similar signaltonoise levels across beams in each 25 km segment (Fig. 5a–f).
The results of the GFT (Sect. 2.3) are used to delineate ATL03 photon heights between wave and seaice surface variability. We construct a binned wave height field along the track from the GFTderived surface spectrum by filtering out highwavenumber components that likely do not correspond to swell waves. In Fig. 9 we show the identified lowpass filters and the displacement spectrum (m^{2} k^{−1}) rather than the slope spectrum ($(\mathrm{m}/\mathrm{m}{)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{k}}^{\mathrm{1}}$, as in Fig. 8) to better separate the highfrequency noise from the lowerfrequency waves. The lowpass filter is defined by a cutoff wavenumber ${k}_{\mathrm{c}}^{\prime}$, the first wavenumber where the observed power spectrum changes slope. A change in the slope of the displacement spectrum in log–log scaling from the expected slope of surface wave spectra (${k}^{\mathrm{5}/\mathrm{2}}$ or similar; Toba, 1973) to horizontal indicates a change in the signaltonoise regime in the data (Fig. 9). Hence, horizontal slopes at high wavenumber indicate Gaussian (white) noise, while steeper slopes at lower wavenumber result from wave–wave interaction (Kitaigorodskii, 1962; Hasselmann et al., 1973). The critical wavenumber ${k}_{\mathrm{c}}^{\prime}$ between both regimes is found using piecewise regression on the weighted crossbeam log–log power spectrum (Fig. 9; Pilgrim, 2021). In cases where the piecewise regression fails to identify a separation between a steep and horizontal slope, no primary wave spectrum is identified and no lowpass filter is applied (Fig. 8b, X>87.5 km).
For illustrative purposes, we define a lowpass filter by setting ${k}_{\mathrm{c}}^{\prime}$ as the cutoff wavenumber. This lowpass filter potentially creates artificial ringing in real space, and for better results, this should be replaced by a more complex filter design. Here, wavenumbers higher than ${k}_{\mathrm{c}}^{\prime}$ are excluded from the wave height model of all beams (Fig. 10a, gray and blue lines and the gray area in the inlet) by truncating the wavenumber space of the slope model p. From this truncated slope model, we can directly construct a coefficient matrix for the waveheight model $\widehat{\mathit{z}}$ for each beam by integrating in wavenumber space. The coefficient matrix of the waveheight model $\widehat{\mathit{z}}$ is
where c_{c} and a_{c} are the model amplitudes corresponding to the cutoff frequency ${k}_{\mathrm{c}}^{\prime}$ (note the integration of the trig formula changes the order and sign of the indices). The reconstructed waveheight model $\widehat{\mathit{z}}$ can then be directly calculated from the original regressor matrix:
with ${\mathit{k}}_{\mathbf{c}}^{\mathbf{\prime}}=[{k}_{\mathrm{1}}^{\prime},{k}_{\mathrm{2}}^{\prime},\mathrm{\dots},{k}_{\mathrm{c}}^{\prime},{k}_{\mathrm{1}}^{\prime},\mathrm{\dots},{k}_{\mathrm{c}}^{\prime}{]}^{\mathrm{T}}$ as shown in Fig. 10b (blue line). The residual between the height model $\widehat{\mathit{z}}$ and the observed smoothed photon heights z, ${\mathit{z}}_{\text{free}}=\mathit{z}\widehat{\mathbf{z}}$, is an estimate of the freeboard height in the absence of the influence of waves (Fig. 10d). The residual z_{free} has similar data density to the original ATL03 photon retrievals but may reveal secondary, nonwavelike structures in the photon heights as shown in Fig. 10d. We provide additional examples in Figs. S9 and S10 in the Supplement.
Decomposing heights into wave and seaice components allows us to estimate the fraction of the total height variance that is due to neither waves nor photon variance on scales shorter than the 10 m stencil. As shown in Fig. 10e, the majority of the total variance is due to the photon variance around its 20 m stencil mean for scales smaller than the stencil (Fig. 10d, red line and black dots; Kwok et al., 2021). In this particular track, wave variance comprises then another 20 % to 50 % of total photon height variance. The remaining variance, about 5 % to 20 %, is then due to neither waves nor the photon cloud. It is from differences in the observed and modeled wave heights, $\widehat{\mathit{z}}$ and z, and we assign this to seaicerelated variability.
The distribution of the residual statistics is, by model design, approximately Gaussian (Fig. 3e), and hence nonwave signals with nonlinear imprint could contaminate the wave estimate and decomposition. If the contribution of nonwave signals to the acrossbeam average is minor, this decomposition removes waves as the dominant source of variance on scales larger than 20 m. This allows for additional analysis of the residual signal and more consistent surface height signals in waveaffected and lowseaice regions. A better filter design can further improve this separation between waves and sea ice.
ICESat2 photon data frequently show wavelike signals in sea ice, and these substantially impact the Marginal Ice Zone. In this paper, we show, for the first time, how paired laser altimeter observations can be converted into directional surface wave spectra. We describe a twopart algorithm that efficiently decomposes the IS2 photon retrievals into a surface wave signal as well as variability due to sea ice. The first (GFT) part of the algorithm is based on a linear inversion method to fit wavenumber coefficients to the ATL03 data (Sects. 2.2–2.4), and the second (MH) part uses a nonlinear inversion method that optimizes for most likely wave incident angles (Sect. 3.1).
The combined method provides a highly resolved 2D surface wave spectra every 12.5 km along each IS2 track (Fig. 8a and b) as well as an improved surface height estimate when the wave signal is removed (Fig. 10). The surface wave estimate relies on the one hand on the redundancy across beams to optimize the signaltonoise ratio in wavenumber space and on the other hand on the difference across beams for the angle inversion. The iterative solution proposed here leads to an interpretation of the IS2 track as a streak of twodimensional wave spectra, including error estimates on each variable (Fig. 8c–e).
We identify the range of wavenumbers that contain wave energy in each segment by establishing a dynamic noise level (Sect. 4). Removing the wave energy as a dominant source of variance reveals additional structure in the ATL03 photon cloud data that is not as readily present in the ATL07 or other higherlevel products (Fig. 10d, or Figs. S9 and S10), because either it is obscured by wave signals or the data are not present. Even though we do not investigate the residual photon heights further, we believe that a removal of the wave signal may have substantial benefits for understanding the seaice structure and classifying photon data for ATL07 products and above. Not removing the wave signal likely leads to an aliasing effect of the waves into the freeboard height (compare panels b and d in Fig. 10).
Our quantification of wave energy allows for an improved understanding of the observed surface elevations in seaicecovered areas. We showed that ocean surface waves have an important contribution to the variance in the MIZ. While the example tracks have a substantial amount of variance in photon height on scales smaller than 20 m (Fig. 11a, c, and e, gray area), the variance on scales longer than 20 m is clearly dominated by the effect of waves (Fig. 11a, c, and e, blue area). Especially in the MIZ, only a small fraction of this rebinned variance is due to nonwavelike features of the surface (Fig. 11a and e, grayblue area at 0–50 km from the ice edge).
The chosen examples show a clear wave signal that can be separated from highfrequency noise by using a simple lowpass filter with a cutoff frequency ${k}_{\mathrm{c}}^{\prime}$ (Sect. 4). This cutoff frequency assumes a scale separation between waves and seaice variance that generally may not exist. The identified cutoff frequency lies in the plausible range of wind waves (k=0.05 to 0.1). In cases with wind waves and complex seaice structures (Alberello et al., 2022), a separation between the wave and seaice signal is hampered. Without adding observations other than IS2, we are not able to directly differentiate between wind sea and swell.
This waveinduced variance in the photon cloud of ATL03 can, under certain conditions, also be captured by ATL07. However, we suggest that the ATL07 algorithm also can potentially capture an aliased wave signal or can fail to provide a sufficient seaice height product at all. While in the case of Track 3 (Fig. 11a and b), waves clearly affect the observation, ATL07 does not provide data in this region (see also Fig. S1) and so our inversion method using ATL03 allows for improved MIZ freeboard data. In other cases, like Track 2 (Fig. 11d), ATL07 does provide data at some but not all locations along the track, even though there is a high photon density throughout the track.
In the MIZ, ATL07 variance can exceed the variance estimates of waves, indicating potential aliasing of waveinduced signal to other scales (Fig. 11d). This aliasing can be due to the binning of data based on photon counts which result in varying bin length. Varying bin length potentially subsamples the wave's energy at scales on or around the Nyquist frequency of the dominant waves. Since both the wavelength and photon density highly vary, it is generally unknown whether or not the photocountbased measure samples these wavelengths correctly and does not negatively affect freeboard retrievals or waveenergy estimates in the MIZ.
5.1 Applying the generalized Fourier transform method for alongtrack wave spectra
We chose a generalized Fourier transform (GFT) method for the wave field inversion instead of more common methods like DFT or Lomb–Scargle (LS; Sect. 2.2; Lomb, 1976; Scargle, 1982; Wunsch, 1996; Kachelein et al., 2022). In contrast to the DFT and LS, the GFT is a varianceconserving method that can be applied to unstructured data and does not require periodicity over an (arbitrary) window length.
As with any linear inverse method, the GFT assumes Gaussian statistics, which is obeyed by linear waves but potentially violated by seaice surface variability. To minimize the effect of seaice heights on the wave inversion, we use mean surface slopes rather than heights (Sect. 2.1), which results in an approximate Gaussian residual (Fig. 3e).
The GFT is customizable to the wavenumbers of interest and additionally provides uncertainty bounds on all parameters. In turn, it comes with the requirement to a priori know what spectral resolution is needed. Given Atlas's highresolution photon cloud data, we chose to resolve the plausible wavenumber range for surface waves on a resolution about twice the one of the DFT. Other, more targeted, narrowbanded wavenumber spaces are possible, but here we choose wavenumber ranges that are the most general for surface waves. Higher resolution, especially in the long wavenumber range, allows us to provide new insights into how narrowbanded the surface wave field is – a parameter that is related to the surface's curvature and likely important for waveinduced seaice breakup (Meylan and Squire, 1994).
A major advantage of the GFT is that it can be extended to inversions of the wave field for each IS2 track by coupling neighboring or overlapping segments, similar to Kalman inversion methods. We illustrate this by simply iterative updating the data segments and model priors (Sect. 2.4 and Fig. 4). This coupling of inversions along segments advances the task on hand to the task of solving the coupled inversions consistently along each beam rather than independently for each segment. The same idea can be extended by coupling the model prior P across beams. This coupled approach ensures smoothly varying wave statistics along and acrosstrack, with the amount of smoothness tuned through the amplitude of P in Eq. (5). Future comparison with other observations will help to constrain the amplitude of P.
5.2 Applying the Metropolis–Hastings method for wave angle inversion
The inversion for the wave's incident angle is based on the crossbeam phase lag (Sect. 3). The phase lag between two beams is limited by the geometry (Fig. 1b), a difficulttoestimate property of the wave field (i.e., its “groupness”, Fig. S7; LonguetHiggins and Deacon, 1957), as well as observational noise (Sect. 3.1). We choose to approach this low signaltonoise problem with a supersample of marginal distributions derived from independent MCMC samples of monochromatic plane waves. The unweighted mean of this method across all wavenumbers is similar to the lag crosscorrelation of the beam pair. However, by focusing on a limited set of energycontaining wavenumbers, the signaltonoise ratio improves above a lagged crosscorrelation approach. Limiting the sample to the most energetic waves and using a prior raises the signaltonoise ratio and is what enables an inversion of the approximate wave angle (Fig. 7d).
The quality of the MH inversion method depends on the wavelength, wave amplitude, and curvature of the wave spectrum. The longer the wave, the better the phase lag can be observed, but those are not the most energetic. In turn, the most energetic waves have typically shorter wavelengths that are of 80–250 times the segment length (25 km), which can lead to multiple minima in the optimization due to a 2π ambiguity. Finally, the curvature of the wave spectrum characterized the length of wave groups, which in 2D erode the ability to observe the average phase lag between the two beams (second bullet point in Sect. 3.1 and Fig. S7).
The inversion is limited by the geometry of the observation (Fig. 1c). Waves coming from steep angles relative to the IS2 track cannot be resolved, such that our 2D wave field estimates are limited in range (about $\pm \mathrm{75}{}^{\circ}$). This problem may be overcome by using better, observationally based priors or enriching our WW3 priors with data from other sources.
The angle inversion method generally can sample multiple incident angles, i.e., multiple minima in the objective function, and may be able to detect multiple wave systems from different directions (Alberello et al., 2022). However, after testing the effect of wave groups and uncorrelated noise on the sampling method, we came to the conclusion that the signaltonoise ratio is too low for a frequencydependent angle correction, even after adding a prior. Hence, for the purpose of angle correction, we here assume that all wave energy comes from the same angle, i.e., the angle of the most energetic waves. Even though the true wave field may be a superposition of multiple wave systems with varying directions, the single incident angle is justifiable here because we are focused on the attenuation and propagation of the dominant wave energy.
The low signal to noise of the angle inversion requires regularization (Sect. 3.1). Since directional wave observations colocated with IS2 tracks are sparse and not readily available, we relied on WW3 hindcast models as a prior (IOWAGA; WW3DG, 2009). The wave hindcasts may perform sufficiently well in the Northern Hemisphere but are known to have limitations in the Southern Ocean MIZ, potentially due to wind biases (Belmonte Rivas and Stoffelen, 2019; Hell et al., 2020, 2021b). The lack of certainty in WW3's peak direction and frequency is expressed in the value of the hyperparameter β_{θ} (Eq. 11). A value of β_{θ}=2 leads to the desired behavior of breaking the symmetry (compare shading in Fig. 6b and d) but not imposing the optimization result through the prior (Fig. 7d, blue and orange lines).
The proposed MCMC method shares aspects with the Wavelet Directional Method (WDM; Donelan et al., 1996, 2015), which decomposes the signals of at least three stationary wave observations into wavelets for each frequency. Similar to our method, WDM uses the phase lag of the wavelets between the three stations to identify a wave incident angle per frequency. WDM could be applied to transects of the wave surface as present in our analysis. However, IS2 only provides two neighboring laser beams, and other beam pairs are too distant (about 3 km) for coherent phase analysis. In addition, the signal to noise may be substantially lower in the IS2 observations, as wave crests are potentially distorted by seaice structure. Therefore, we introduced a wave angle prior (Eq. 11) to break the ambiguity in the observed phase lag (Fig. 6b and d, shading).
We proposed and tested a method to decompose photon retrievals from the IS2 satellite into a surface wave signal and residual variance. The surface wave signal is identified using a generalized Fourier transform, and it can be expressed as a directional surface wave spectrum by adding a Metropolis–Hastings sampling method to identify the incident angle of the dominant wave energy (Sect. 3). The wave and seaice signal is separated using a simple cutoff frequency that implies a separation of scales.
Surface waves and sea ice can have complex nonlinear interactions that are important to model for improving seaice and climate projections. This method will enable us to observe the interaction between the dominant waves and sea ice by utilizing large datasets from IS2. However, IS2's nearly onedimensional snapshots of the surface height can hardly capture all possible interactions between waves and ice. Besides gaps in the observations due to clouds and varying photon densities between seaice types, a correct wave field inversion is only possible with sufficient data density and a limited range of incident angles (Sect. 3). Even though this method outlines a better, more transparent wavefield inversion than a DFT, it remains to be seen how the interaction of those limitations can be used to provide a highly resolved global waveinice product. Comparisons with other data sources, from either in situ or remote sensing observations, are needed to understand these limitations better and validate this method.
Waves and sea ice have scales ranging multiple orders of magnitude such that it is challenging to separate both in the IS2 observations. The choice of the parameters in this analysis (10 m bins, 25 km segment length, and the slopebased cutoff frequency ${k}_{\mathrm{c}}^{\prime}$) focuses on identifying swell wave events routinely created by synoptic storms (Hell et al., 2021a). However, even on these scales (80 to 350 m), a separation between wave and seaice signal may only be possible when the seaice variance is weak on those scales and the data are not too gappy, as in the chosen example tracks (Figs. 3 and 10). High levels of seaice variance or frequent data gaps will lead to systematic biases and aliasing effects in the wave spectral estimates. To identify these more complex cases, we proposed mitigation strategies that exploit the fact that swell spectra normally vary over larger scales than the segment length (12.5 km) or the separation distance of the beam pairs (3 km): by applying an iterative inversion (Fig. 4); using crossbeam average (Fig. 5); and providing error estimates in frequency, real space, and direction (Eqs. 8 and 9, Figs. 5i and 7b), the method provides ample auxiliary data to detect unusual features in the observations. Complex cases that show a large spread between beams, or large errors, can be identified and excluded from further analysis (Fig. S4).
Despite the shortcomings and limitations of individual track inversions, the method has promise, especially when applied at scale. We expect 10 %–15 % of the IS2 tracks in polar regions to be dominated by waves (Brouwer et al., 2022), which means there is already a large collection of diverse wave observations from IS2. Because IS2 can also record freeboard heights, floe sizes, and seaice types, this analysis can provide complementary seaice information to constrain dynamics in the MIZ. This will be used to statistically constrain parametrizations of wave attenuation in sea ice (Fig. 5) and can potentially leverage the waveremoved residual signal to improve ice classification (Fig. 10d).
Finally, while the heredeveloped method is customized for IS2 photon retrievals, this approach applies to any unstructured quasiinstantaneous observation of the ocean surface. In the case of IS2, crosstrack information is limited, but other future remote sensing methods may have complementary information about the surface wave field. Such an inversion could combine data from IS2, SAR, and CFOSAT to help constrain the surface wave filed in the MIZ or the open ocean (Collard et al., 2022).
A1 Data prior R
We define the data prior R based on the surface slope uncertainty as
where the stencil width is Δx=20 m, σ(b) is the standard deviation of the data b within the segment, σ_{h} is the vector of standard deviation of each data point (each stencil, Sect. 2.1), and β_{R} is a tuning parameter that determines the ratio of the model and data priors in Eq. (5) (Fig. 2, blue lines). The standard deviation, or error, of the data is divided by the stencil size to get an error in units of surface slope, and the variance of the data is then used to amplify the prior R to scale it against the model prior.
To avoid overfitting the ratio of the model prior P and data prior R has to be adjusted. For this we try different values of β_{R} and set it to 10^{2} such that the distribution of the residual r is approximately Gaussian and that
as shown in Fig. 3e. At locations with no data, this results in a model decay to zero (Fig. 3 orange lines) on scales similar to the data's autocorrelation.
A2 Model prior P
As described in Sect. 2.4, the GFT's prior for initial segments are derived from a PM spectrum that is fitted to a DFT of the segment data (Fig. A1b, dashed gray and black lines). The initial prior P_{init} is then defined as the peaknormalized PM spectrum multiplied by the data variance σ(b)^{2} plus a 10 % noise floor (Fig. A1c, dashed black line). The initial prior is used to perform a first inversion of this segment, but – to avoid overfitting to the PM spectrum – a second inversion with the same data is done, but now with a “data” prior, i.e., the smoothed result of the first inversion (Fig. A1c, black line). The power of the second (final) GFT model coefficients $\widehat{\mathit{p}}$ is then shown in green (Fig. A1c).
Segments with a successful inversion in the previous segment do not make use of the PMbased prior; instead, they use the data prior from the previous segment: for a segment i with a successful inversion in the previous segment i−1, we use a smoothed power spectrum based on ${\widehat{\mathit{p}}}_{i\mathrm{1}}$ to derive P_{i} (Fig. A1e and g, black lines, and Sect. 2.4). Note that even though the initial PM prior pushes the model to a single peak spectrogram (Fig. A1b and c, dashed line), in cases where the data do not support this shape, as in this example, the PM prior does not determine the result (Fig. A1c, e, and g, green lines). Any spectral shape, including double peaks, can appear. Cases where a PM prior improves the inversion are shown in Fig. S12 in the Supplement.
The prior in Eq. (11) uses an incident angle θ_{0}(k) with an uncertainty σ_{θ}(k) defining the prior as
Both variables in the prior have to be taken from data sources other than IS2, and here they are derived from WAVEWATCH III (WW3) global hindcast wave partitions (WW3DG, 2009, using the Integrated Ocean Waves for Geophysical and other Applications (IOWAGA) hindcast) and depend on wavenumber. The WW3 hindcast data are selected in a box around the most equatorward photons in IS2 (Fig. S11 and Sect. 2.1). In cases where this box is within the seaice mask of WW3, it is moved equatorward along the IS2 track until at least $\mathrm{2}/\mathrm{3}$ of the box's grid cells are not covered by the WW3 seaice mask.
Within each box, the mean of the peak period, peak direction, and directional spread is calculated for each of the five WW3 wavenumber partitions (Fig. 7a, black dots). These partitions are interpolated and smoothed to the k_{i} wave numbers of interest to best guess the wavenumberdependent peak direction and spread (Fig. 7, an orange line and shading). Note that it is hard to validate the WW3 partitions in the Southern Ocean due to a lack of contemporaneous in situ observations. The lack of certainty in the wave hindcast, in combination with (any) smoothing procedure, can lead to biases in the directional prior. Hence, we do not expect the direct alignment of the WW3 wave directions and those observed in IS2 observations. The WW3 incident angle is here solely used to reduce ambiguities in the objective function (Figs. 6b and d and 7) and can give a preferred incident wave angle rather than a certain estimate of the dominant wave angle (Fig. 7).
The algorithms are available through Hell (2022a, DOI: https://doi.org/10.5281/zenodo.6908645), and data are available through Hell (2022b, DOI: https://doi.org/10.5281/zenodo.6928350). The IS2 ATL03 data (https://doi.org/10.5067/ATLAS/ATL03.005; Neumann et al., 2021) and ATL07/10 data (https://doi.org/10.5067/ATLAS/ATL07.003; Kwok et al., 2021) are available through NSIDC (https://nsidc.org/data/icesat2, last access: 5 October 2022). The wave model data are available through the Integrated Ocean Waves for Geophysical and Other Applications (IOWAGA) project: https://tds3.ifremer.fr/thredds/catalog.html (Ardhuin et al., 2011). The analysis uses and modifies icesat2toolkit (https://readicesat2.readthedocs.io/, last access: 31 June 2021, https://doi.org/10.5281/zenodo.7439353, Sutterley and Siegfried, 2019).
The supplement related to this article is available online at: https://doi.org/10.5194/tc183412024supplement.
MCH developed and programmed the proposed algorithm. MCH and CH collaboratively wrote the paper.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We thank Bruce Cornuelle and BaylorFox Kemper for discussing the depths and caveats of the proposed algorithms.
Christopher Horvat and Momme C. Hell were supported by NASA (grant nos. 80NSSC20K0959 and 80NSSC23K0935) and Schmidt Futures – a philanthropic initiative that seeks to improve societal outcomes through the development of emerging science and technologies. Christopher Horvat was supported by NSF OPP2146889.
This paper was edited by Jari Haapala and reviewed by Marcello Vichi and one anonymous referee.
Alberello, A., Bennetts, L. G., Onorato, M., Vichi, M., MacHutchon, K., Eayrs, C., Ntamba, B. N., Benetazzo, A., Bergamasco, F., Nelli, F., Pattani, R., Clarke, H., Tersigni, I., and Toffoli, A.: ThreeDimensional Imaging of Waves and Floes in the Marginal Ice Zone during a Cyclone, Nat. Commun., 13, 4590, https://doi.org/10.1038/s41467022320362, 2022. a, b, c
Ardhuin, F., Hanaﬁn, J., Quilfen, Y., Chapron, B., Queﬀeulou, P., Obrebski, M., Sienkiewicz, J., and Vandemark, D.: Calibration of the “IOWAGA” Global Wave Hindcast (1991–2011) Using ECMWF and CFSR Winds, https://tds3.ifremer.fr/thredds/catalog.html (last access: 19 November 2021), 2011. a
Ardhuin, F., Stopa, J., Chapron, B., Collard, F., Smith, M., Thomson, J., Doble, M., Blomquist, B., Persson, O., Collins, C. O., and Wadhams, P.: Measuring ocean waves in sea ice using SAR imagery: A quasideterministic approach evaluated with Sentinel1 and in situ data, Remote Sens. Environ., 189, 211–222, https://doi.org/10.1016/j.rse.2016.11.024, 2017. a
Belmonte Rivas, M. and Stoffelen, A.: Characterizing ERAInterim and ERA5 surface wind biases using ASCAT, Ocean Sci., 15, 831–852, https://doi.org/10.5194/os158312019, 2019. a, b
Brouwer, J., Fraser, A. D., Murphy, D. J., Wongpan, P., Alberello, A., Kohout, A., Horvat, C., Wotherspoon, S., Massom, R. A., Cartwright, J., and Williams, G. D.: Altimetric observation of wave attenuation through the Antarctic marginal ice zone using ICESat2, The Cryosphere, 16, 2325–2353, https://doi.org/10.5194/tc1623252022, 2022. a, b
Collard, F., Marié, L., Nouguier, F., Kleinherenbrink, M., Ehlers, F., and Ardhuin, F.: WindWave Attenuation in Arctic Sea Ice: A Discussion of Remote Sensing Capabilities, J. Geophys. Res.Oceans, 127, e2022JC018654, https://doi.org/10.1029/2022JC018654, 2022. a, b, c
Donelan, M., Babanin, A., Sanina, E., and Chalikov, D.: A Comparison of Methods for Estimating Directional Spectra of Surface Waves, J. Geophys. Res.Oceans, 120, 5040–5053, https://doi.org/10.1002/2015JC010808, 2015. a
Donelan, M. A., Drennan, W. M., and Magnusson, A. K.: Nonstationary Analysis of the Directional Properties of Propagating Waves, J. Phys. Oceanogr., 26, 1901–1914, 1996. a
ForemanMackey, D., Hogg, D. W., Lang, D., and Goodman, J.: Emcee: The MCMC Hammer, Publ. Astron. Soc. Pac., 125, 306–312, https://doi.org/10.1086/670067, 2013. a, b
Hasselmann, K., Barnett, T. P., Bouws, E., Carlson, H., Cartwright, D. E., Enke, K., Ewing, J. A., Gienapp, H., Hasselmann, D. E., Kruseman, P., Meerburg, A., Müller, P., Olbers, D. J., Richter, K., Sell, W., and Walden, H.: Measurements of WindWave Growth and Swell Decay during the Joint North Sea Wave Project (JONSWAP), Ergänzung zur Deut. Hydrogr. Z., Reihe A, 12, 1–95, 1973. a
Hell, M.: Code for Directional Surface Wave Spectra And Sea Ice Structure from ICEsat2 Altimetry without Data, Zenodo [code], https://doi.org/10.5281/zenodo.6908645, 2022a. a
Hell, M.: Data for Directional Surface Wave Spectra And Sea Ice Structure from ICEsat2 Altimetry, Zenodo [data set], https://doi.org/10.5281/zenodo.6928350, 2022b. a
Hell, M. C., Cornuelle, B. D., Gille, S. T., Miller, A. J., and Bromirski, P. D.: Identifying Ocean Swell Generation Events from Ross Ice Shelf Seismic Data, J. Atmos. Ocean. Tech., 36, 2171–2189, https://doi.org/10.1175/JTECHD190093.1, 2019. a
Hell, M. C., Gille, S. T., Cornuelle, B. D., Miller, A. J., Bromirski, P. D., and Crawford, A. D.: Estimating Southern Ocean Storm Positions With Seismic Observations, J. Geophys. Res.Oceans, 125, e2019JC015898, https://doi.org/10.1029/2019JC015898, 2020. a, b
Hell, M. C., Ayet, A., and Chapron, B.: Swell Generation Under ExtraTropical Storms, J. Geophys. Res.Oceans, 126, e2021JC017637, https://doi.org/10.1029/2021JC017637, 2021a. a, b
Hell, M. C., Cornuelle, B. D., Gille, S. T., and Lutsko, N. J.: TimeVarying Empirical Probability Densities of Southern Ocean Surface Winds: Linking the Leading Mode to SAM and Quantifying Wind Product Differences, J. Climate, 34, 5497–5522, https://doi.org/10.1175/JCLID200629.1, 2021b. a
Horvat, C.: Floes, the Marginal Ice Zone, and Coupled WaveSeaIce Feedbacks, Philos. T. Roy Soc. A, 380, 20210252, https://doi.org/10.1098/rsta.2021.0252, 2022. a
Horvat, C., Roach, L. A., Tilling, R., Bitz, C. M., FoxKemper, B., Guider, C., Hill, K., Ridout, A., and Shepherd, A.: Estimating the sea ice floe size distribution using satellite altimetry: theory, climatology, and model comparison, The Cryosphere, 13, 2869–2885, https://doi.org/10.5194/tc1328692019, 2019. a
Horvat, C., BlanchardWrigglesworth, E., and Petty, A.: Observing waves in sea ice with ICESat2, Geophys. Res. Lett., 47, e2020GL087629, https://doi.org/10.1029/2020GL087629, 2020. a, b
Kachelein, L., Cornuelle, B. D., Gille, S. T., and Mazloff, M. R.: Harmonic Analysis of NonPhaseLocked Tides with Red Noise Using the Red_tide Package, J. Atmos. Ocean. Tech., 39, 1031–1051, https://doi.org/10.1175/JTECHD210034.1, 2022. a, b, c
Kitaigorodskii, S. A.: Applications of the Theory of Similarity to the Analysis of WindGenerated Wave Motion as a Stochastic Process, Izv. Geophys. Ser. Acad. Sci. USSR, 1, 105–117, 1962. a, b
Kwok, R., Petty, A., Cunningham, G., Markus, T., Hancock, D., Morison, J. H., Palm, S. P., Farrell, S. L., Ivanoff, A., Wimert, J., and ICESat2 Science Team: ATLAS/ICESat2 L3A Sea Ice Height, Version 3 [ATL07/10], NSIDC [data set], https://doi.org/10.5067/ATLAS/ATL07.003, 2021. a, b, c, d, e
Lomb, N. R.: LeastSquares Frequency Analysis of Unequally Spaced Data, Astrophys. Space Sci., 39, 447–462, https://doi.org/10.1007/BF00648343, 1976. a
LonguetHiggins, M. S.: Statistical Properties of Wave Groups in a Random Sea State, Philos. T. R. Soc., 312, 219–250, https://doi.org/10.1098/rsta.1984.0061, 1984. a
LonguetHiggins, M. S. and Deacon, G. E. R.: The Statistical Analysis of a Random, Moving Surface, Philos. T. R. Soc., 249, 321–387, https://doi.org/10.1098/rsta.1957.0002, 1957. a, b
Lygre, A. and Krogstad, H. E.: Maximum Entropy Estimation of the Directional Distribution in Ocean Wave Spectra, J. Phys. Oceanogr., 16, 2052–2060, https://doi.org/10.1175/15200485(1986)016<2052:MEEOTD>2.0.CO;2, 1986. a
Marechal, G. and Ardhuin, F.: Surface Currents and Significant Wave Height Gradients: Matching Numerical Models and HighResolution Altimeter Wave Heights in the Agulhas Current Region, J. Geophys. Res.Oceans, 126, e2020JC016564, https://doi.org/10.1029/2020JC016564, 2021. a
Menke, W.: Geophysical Data Analysis: Geophysical Data Analysis: Discrete Inverse Theory, 4th Edn., Elsevier, Amsterdam, 2018. a
Meylan, M. and Squire, V. A.: The Response of Ice Floes to Ocean Waves, J. Geophys. Res.Oceans, 99, 891–900, https://doi.org/10.1029/93JC02695, 1994. a
Neumann, T. A., Brenner, D., Hancock, D., Robbins, J., Saba, K., Harbeck, K., Gibbons, A., Lee, J., Luthcke, S. B., and Rebold, T.: ATLAS/ICESat2 L2A Global Geolocated Photon Data, Version 5. [ATL03], NSIDC [data set], https://doi.org/10.5067/ATLAS/ATL03.005, 2021. a, b
Pierson, W. J. and Moskowitz, L.: A Proposed Spectral Form for Fully Developed Wind Seas Based on the Similarity Theory of S. A. Kitaigorodskii, J. Geophys. Res., 69, 5181–5190, https://doi.org/10.1029/JZ069i024p05181, 1964. a
Pilgrim, C.: PiecewiseRegression (Aka Segmented Regression) in Python, J. Open Source Softw., 6, 3859, https://doi.org/10.21105/joss.03859, 2021. a
Rapley, C. G.: First Observations of the Interaction of Ocean Swell with Sea Ice Using Satellite Radar Altimeter Data, Nature, 307, 150–152, https://doi.org/10.1038/307150a0, 1984. a, b
Scargle, J. D.: Studies in Astronomical Time Series Analysis. II. Statistical Aspects of Spectral Analysis of Unevenly Spaced Data, Astrophys. J., 263, 835–853, https://doi.org/10.1086/160554, 1982. a
Squire, V. A.: Of ocean waves and seaice revisited, Cold Reg. Sci. Technol., 49, 110–133, https://doi.org/10.1016/j.coldregions.2007.04.007, 2007. a
Squire, V. A.: A Fresh Look at How Ocean Waves and Sea Ice Interact, Philos. T. Roy. Soc. A, 376, 20170342, https://doi.org/10.1098/rsta.2017.0342, 2018. a
Stopa, J. E., Sutherland, P., and Ardhuin, F.: Strong and highly variable push of ocean waves on Southern Ocean sea ice, P. Natl. Acad. Sci. USA, 115, 5861–5865, https://doi.org/10.1073/pnas.1802011115, 2018. a
Sutterley, T. C. and Siegfried, M.: icesat2toolkit: Python tools for obtaining and working with elevation data from the NASA ICESat2 mission, Zenodo [code], https://doi.org/10.5281/zenodo.7439353, 2019. a
The WAVEWATCH III R© Development Group (WW3DG): User manual and system documentation of WAVEWATCH III R© version 6.07. Tech. Note 333, NOAA/NWS/NCEP/MMAB, College Park, MD, USA, 465 pp. + Appendices, 2019. a, b, c
Thomson, J.: Wave propagation in the Marginal Ice Zone: connections and feedbacks within the airiceocean system, Philos. T. Roy. Soc. A, 380, 20210251, https://doi.org/10.1098/rsta.2021.0251, 2022. a, b
Thomson, J., Hošeková, L., Meylan, M. H., Kohout, A. L., and Kumar, N.: Spurious Rollover of Wave Attenuation Rates in Sea Ice Caused by Noise in Field Measurements, J. Geophys. Res.Oceans, 126, e2020JC016606, https://doi.org/10.1029/2020JC016606, 2021. a
Toba, Y.: Local Balance in the AirSea Boundary Processes, J. Oceanogr. Soc. Jpn., 29, 209–220, https://doi.org/10.1007/BF02108528, 1973. a
Tsallis, C. and Stariolo, D. A.: Generalized Simulated Annealing, Physica A, 233, 395–406, https://doi.org/10.1016/S03784371(96)002713, 1996. a
Villas Bôas, A. B. and Young, W. R.: Directional Diffusion of Surface Gravity Wave Action by Ocean Macroturbulence, J. Fluid Mech., 890, E3, https://doi.org/10.1017/jfm.2020.116, 2020. a
Wunsch, C.: The Ocean Circulation Inverse Problem, Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511629570, 1996. a, b, c
Yu, Y., Sandwell, D. T., Gille, S. T., and Villas Bôas, A. B.: Assessment of ICESat2 for the Recovery of Ocean Topography, Geophys. J. Int., 226, 456–467, https://doi.org/10.1093/gji/ggab084, 2021. a, b
 Abstract
 Introduction and problem description
 Alongtrack wave spectra from IS2
 Twodimensional wave spectra from nearly onedirectional observations
 Isolating the wave signal in the alongtrack data
 Discussion
 Conclusion
 Appendix A: GFT priors
 Appendix B: WAVEWATCH III prior
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction and problem description
 Alongtrack wave spectra from IS2
 Twodimensional wave spectra from nearly onedirectional observations
 Isolating the wave signal in the alongtrack data
 Discussion
 Conclusion
 Appendix A: GFT priors
 Appendix B: WAVEWATCH III prior
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement