Measuring ocean total surface current velocity with the KuROS and KaRADOC airborne near-nadir Doppler radars: a multi-scale analysis in preparation for the SKIM mission

Surface currents are poorly known over most of the world's oceans. Satellite-borne Doppler wave and current scatterometers (DWaCSs) are among the proposed techniques to fill this observation gap. The Sea surface KInematics Multiscale (SKIM) proposal is the first satellite concept built on a DWaCS design at near-nadir angles and was demonstrated to be technically feasible as part of the European Space Agency Earth Explorer program. This article describes preliminary results from a field experiment performed in November 2018 off the French Atlantic coast, with sea states representative of the open ocean and a well-known tide-dominated current regime, as part of the detailed design and feasibility studies for SKIM. This experiment comprised airborne measurements performed using Ku-band and Ka-band Doppler radars looking at the sea surface at near-nadir incidence in a real-aperture mode, i.e., in a geometry and mode similar to that of SKIM, as well as an extensive set of in situ instruments. The Ku-band Radar for Observation of Surfaces (KuROS) airborne radar provided simultaneous measurements of the radar backscatter and Doppler velocity in a side-looking configuration, with a horizontal resolution of about 5 to 10 m along the line of sight and integrated in the perpendicular direction over the real-aperture 3 dB footprint diameter (about 580 m). The Ka-band RADar for Ocean Current (KaRADOC) system, also operating in the side-looking configuration, had a much narrower beam, with a circular footprint only 45 m in diameter. Results are reported for two days with contrasting conditions, a strong breeze on 22 November 2018 (wind speed 11.5 m s⁻¹, Hs 2.6 m) and gentle breeze on 24 November 2018 (wind speed 5.5 m s⁻¹, Hs 1.7 m). The measured line-of-sight velocity signal is analyzed to separate a non-geophysical contribution linked to the aircraft velocity, a geophysical contribution due to the intrinsic motion of surface waves and the desired surface current contribution. The surface wave contribution is found to be well predicted by Kirchhoff scattering theory using as input parameters in situ measurements of the directional spectrum of long waves, complemented by the short wave spectrum of Elfouhaily et al. (1997). It is found to be closely aligned with the wind direction, with small corrections due to the presence of swell. Its norm is found to be weakly variable with wind speed and sea state, quite stable and close to $C_{\mathrm{0}}=\mathrm{2.0}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ at the Ka band, and more variable and close to $C_{\mathrm{0}}=\mathrm{2.4}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ at the Ku band. These values are 10 %–20 % smaller than previous theoretical estimates. The directional spread of the short gravity waves is found to have a marked influence on this surface wave contribution. Overall, the results of this study support the feasibility of near-nadir radar Doppler remote sensing of the ocean total surface current velocity (TSCV).

1 Introduction

The ocean total surface current velocity (TSCV) is defined as the Lagrangian mean velocity at the instantaneous sea surface, corresponding to an effective mass transport velocity at the surface. The TSCV is currently only reliably measured by high-frequency (HF) radars, which are only deployed in some coastal regions. Elsewhere, available estimates depend on numerical model outputs, sea level and wind measurements, and on assumptions such as the balance between the surface pressure gradient and the Coriolis force. The situation is similar regarding directional wave statistics, which are currently mainly estimated through numerical modeling.

These estimates of the TSCV are not reliable at small scales, particularly so in the tropical ocean (e.g., Sudre et al., 2013; Stopa et al., 2016), and these limitations hamper current efforts to observe and understand the fluxes of heat, fresh water, carbon and plastics as well as the coastal impacts of sea states.

Whereas new data on ocean waves are becoming available with the Surface Waves Investigation and Monitoring (SWIM) instrument carried by the China–France Ocean SATellite (CFOSAT) (Hauser et al., 2017, 2020), direct spaceborne measurements of surface current have been limited to a few regions and single projections of the current vector (Chapron et al., 2005; Rouault et al., 2010; Hansen et al., 2011). Several concepts based on SAR (synthetic-aperture radar) interferometry (Romeiser et al., 2003; Buck, 2005) or Doppler scatterometry (Rodriguez, 2018; Chelton et al., 2019) have been proposed for satellite missions aimed at mapping the ocean surface current vector (see review by Ardhuin et al., 2019a). Airborne demonstrators have also been developed in that context (Martin et al., 2018; Rodríguez et al., 2018) and are now becoming operational tools for oceanographic research.

The Doppler frequency shift (DFS) signal provided by these phase-resolving radar instruments is complex: it contains a geophysical contribution due to waves and currents, as well as a large non-geophysical contribution due to the platform motion. The platform velocity in space being of the order of 7 km s⁻¹ for low Earth orbit, it is obviously critical to have accurate knowledge of the measurement geometry to correctly estimate the non-geophysical component. The contribution due to ocean waves is, however, also an order of magnitude larger than the expected TSCV contribution (Nouguier et al., 2018) and must also be precisely estimated using an accurate sea state description.

The Sea surface KInematics Multiscale monitoring (SKIM) satellite mission has been designed to address all these requirements and provide direct global-coverage measurements of TSCV. Its main instrument payload, the SKIM Ka-band Radar (SKaR), a phase-resolved SWIM-like conically scanning radar, provides simultaneous Ka-band observations of sea state and DFS at 6 and 12^∘ incidence angles as well as state-of-the-art altimetry observations using a dedicated nadir beam. The nadir beam observations are used to control the SkaR acquisition geometry but are also processed using classical algorithms to provide sea surface elevation, significant wave height and wind speed measurements.

SKIM was preselected as one of the two candidate missions for the European Space Agency (ESA) 9th Earth Explorer. As part of the detailed design and feasibility (phase A) studies, ESA funded a dedicated measurement campaign, Drift4SKIM, which was organized from 21 to 27 November 2018 off the French Atlantic coast, in an area with sea states characteristic of the open ocean and a well-known tide-dominated current regime monitored by a two-site 12 MHz high-frequency radar system (Ardhuin et al., 2009; Sentchev et al., 2013). A range of in situ instruments (surface current drifters, drifting and moored wave-measuring buoys), as well as two airborne Doppler radars operating in the Ku (KuROS – Ku-band Radar for Observation of Surfaces) and Ka (KaRADOC – Ka-band RADar for Ocean Current) bands, were deployed. The campaign goals were to

• demonstrate how the non-geophysical contribution V_NG to the DFS can be estimated from the motion of the platform carrying the radar, the antenna diagram properties, and the azimuth and incidence angle dependencies of the radar cross section;

• explore the geophysical component V_GD and its decomposition as a sum of contributions due to currents and waves, V_CD and V_WD (Nouguier et al., 2018); and

• validate the Radar Sensing Satellite Simulator (Nouguier, 2019) and its capability to simulate airborne configurations.

Figure 1 (a) Schematic of the ATR-42 aircraft and KuROS instrument with a definition of viewing angles, azimuth φ and incidence angle θ. (b) A comparison with the SKIM viewing geometry. The unit vector e_φ is the projection on the horizontal of the line-of-sight direction vector. The variation of surface backscatter across the footprint and as a function of azimuth φ, which causes the effective mispointing δφ, is represented as gray shading. In the KuROS data, each measurement is integrated in azimuth across the antenna lobe. In the case of SKIM, the use of unfocused SAR processing allows for the separation of echoes in the azimuth direction with a resolution of dDop ≃300 m.

As highlighted in Fig. 1, the viewing geometry of an airborne system is vastly different from that of a satellite system, with a much smaller footprint and incidence angle variations at scales comparable to the wavelength of the dominant ocean waves. Another obvious difference is the stability of the platform and its velocity, 7 km s⁻¹ for low Earth orbit and around 120 m s⁻¹ for the ATR-42 aircraft used here. As a result, transposing the performance of an airborne system to a satellite system requires a thorough analysis, supplemented by carefully designed and validated simulation tools. Performing this analysis is, however, worthwhile, as it leads one to develop valuable insight into the instrument imaging principle and design trade-offs.

This article is intended to provide an overview of the Drift4SKIM campaign data and a first discussion of their implications for the emerging field of near-nadir Doppler radar observations of TSCV. It is structured as follows: the principle of the pulse-pair measurements and the different contributions to the observed DFS are detailed in Sect. 2 and Appendix A. Section 3 gives a brief account of the fieldwork performed and conditions encountered during the campaign. The results of the airborne measurements are presented in Sect. 4. Results and implications for SKIM are then discussed in Sect. 5. Conclusions and perspectives follow in Sect. 6.

2 Near-nadir radar Doppler measurements of ocean velocities: theory

Shipborne Doppler measurements of ocean currents are routinely performed using so-called vessel-mounted acoustic Doppler current profilers (VMADCPs; see, for instance, Rossby et al., 2019). Some of the data processing concepts transpose directly to the spaceborne context: the raw DFS signal contains a large non-geophysical contribution due to the platform motion, which must be estimated from ancillary sensors and compensated for. The accuracy of the final geophysical product is practically set by the accuracy of the non-geophysical velocity estimation and correction procedure. In the VMADCP context, however, the backscattering elements responsible for the production of the acoustic return signal (particulate suspended matter, zooplanktonic organisms) are passive and accurately follow the water mass. This does not carry over in the electromagnetic case: here, the return signal is produced by the interaction of the transmitted signal with the roughness elements of the sea surface, which move with respect to the water mass with an intrinsic phase velocity that is an order of magnitude larger than typical ocean currents. This effect is, for instance, well known in the ground-based HF radar current measurement context (Stewart and Joy, 1974) and must also be compensated for.

In our case, the measurement geometry is represented in Fig. 1, and the line-of-sight Doppler velocity V_LOS looking towards incidence angle θ and azimuth φ (in this paper, line-of-sight DV contributions are denoted by V, and the corresponding horizontal velocity contributions are denoted by U) is the sum of the projection of a horizontal current contribution U_CD(φ), a wave-induced contribution V_WD(θ,φ) and a non-geophysical contribution V_NG(θ,φ). The equation that permits the retrieval of the TSCV contribution U_CD(φ) from the raw measured V_LOS can be written as

$$\begin{array}{}\text{(1)}& U_{\mathrm{CD}}\left(\mathit{\phi }\right)={\displaystyle \frac{\left[V_{\mathrm{LOS}}(\mathit{\theta },\mathit{\phi })-V_{\mathrm{NG}}(\mathit{\theta },\mathit{\phi })-V_{\mathrm{WD}}(\mathit{\theta },\mathit{\phi })\right]}{\sin \mathit{\theta }}}.\end{array}$$

The aim of this section is to provide a detailed analysis of the different terms of this expression. The non-geophysical contribution V_NG is discussed in Sect. 2.1 and Appendix A. The wave Doppler contribution is discussed in Sect. 2.2. A brief summary of the measurement error budget is finally provided in Sect. 2.3.

2.1 Non-geophysical velocity V_NG

As mentioned above, the accuracy of shipborne acoustic Doppler current measurements is affected in a dominant way by the platform motion compensation process. In the spaceborne context, the platform velocity is almost 3 orders of magnitude larger (7000 m s⁻¹ vs. 10 m s⁻¹ for shipborne measurements). Thus, the accuracy requirements are tremendously exacerbated. In particular, attention must be paid to the detailed effects of the antenna radiation diagram and the sea surface normalized radar cross section (NRCS) variations with space and observation azimuth. A detailed discussion of these effects is given in Appendix A.

In summary, in the case of a sufficiently narrow radiation diagram, V_NG can be approximated as the radar carrier velocity projected on an effective look direction. This effective look direction differs from the geometric boresight direction by an effective azimuthal mispointing δφ due to the finite antenna beamwidth combined with the variations of NRCS within the radar footprint, as well as by an effective incidence angle mispointing δθ due to radar timing or surface-tracking errors.

Table 1 KuROS and KaRADOC antenna radiation diagram characteristics. All angles are in degrees. See Appendix B for the definitions of α and β.

The beamwidth at the working incidence angle is thus a very important parameter of a radar intended for TSCV measurements. Table 1 summarizes the parameters of the KuROS and KaRADOC antennas. For KuROS they have been determined following the procedure detailed in Appendix B. For KaRADOC, they are the result of anechoic chamber measurements (Appendix C). As discussed in Appendix B, these parameters describe the antenna radiation diagrams when expressed as functions of variables, α and β, which do not coincide with azimuth and incidence angle. In the case of constant-altitude flight and near-nadir observations with the antenna looking towards azimuth φ_b, one can, however, obtain a Gaussian approximation of the one-way radiation diagram as

$$\begin{array}{}\text{(2)}& G\simeq \exp \left[-{\displaystyle \frac{{\left(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{b}}\right)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)}{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}}+{\displaystyle \frac{\left({\mathit{\beta }}_{\mathrm{0}}-\tan \left(\mathit{\theta }\right)\right)\tan \left(\mathit{\theta }\right)}{{\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}}}\right]\right],\end{array}$$

where ${\mathit{\sigma }}_{\mathit{\alpha }}={\mathit{\alpha }}_{-\mathrm{3}\mathrm{dB}}/\sqrt{\mathrm{8}\log \left(\mathrm{2}\right)}$ and ${\mathit{\sigma }}_{\mathit{\beta }}={\mathit{\beta }}_{-\mathrm{3}\mathrm{dB}}/\sqrt{\mathrm{8}\log \left(\mathrm{2}\right)}.$ For 12^∘ observations the second term in the exponential can safely be neglected, and the effective azimuthal beamwidth can be estimated as

$$\begin{array}{}\text{(3)}& {\mathit{\phi }}_{-\mathrm{3}\mathrm{dB}}={\displaystyle \frac{{\mathit{\alpha }}_{-\mathrm{3}\mathrm{dB}}}{\sin \left(\mathit{\theta }\right)}}.\end{array}$$

When projected on the ground, φ_−3 dB is thus larger than α_−3 dB by a factor 1∕sin (θ), equal to 4.8 for 12^∘ measurements. Provided that the beam is not too wide, the Gaussian approximation in Eq. (A29) of G as a function of φ can then be used with the parameter

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathit{\phi }}\simeq {\mathit{\alpha }}_{-\mathrm{3}\mathrm{dB}}/\left[\sin \mathit{\theta }\sqrt{\mathrm{8}\log \left(\mathrm{2}\right)}\right].\end{array}$$

Due to the width of the azimuthal aperture, the NRCS-weighted line-of-sight azimuth φ_a can differ from the boresight azimuth φ_b by a mispointing angle δφ. Expressions for δφ are obtained in Appendix A in the two limiting cases of slow linear and fast sinusoidal variations of the ocean surface NRCS with respect to azimuth. In the slow variation case, δφ is obtained as

$$\begin{array}{}\text{(5)}& \mathit{\delta }\mathit{\phi }={\mathit{\phi }}_{\mathrm{a}}-{\mathit{\phi }}_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{0}}}{\partial \mathit{\phi }}}.\end{array}$$

Denoting by φ_t the flight track azimuth and V_p the along-track flight velocity, the spurious azimuth gradient Doppler (AGD) contribution to the DV caused by the mispointing reads

$$\begin{array}{}\text{(6)}& U_{\mathrm{AGD}}=\sin \left({\mathit{\phi }}_{\mathrm{b}}-{\mathit{\phi }}_{\mathrm{t}}\right){\displaystyle \frac{V_{\mathrm{p}}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}{{\sin }^{\mathrm{2}}\mathit{\theta }}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\sigma }}_{\mathrm{0}}}{\partial \mathit{\phi }}}.\end{array}$$

Figure 2 KuROS azimuth integral weight at $\mathit{\theta }={\mathrm{12}}^{\circ }$ for a north-facing (${\mathit{\phi }}_{\mathrm{b}}={\mathrm{0}}^{\circ }$) antenna (black), Gaussian approximation (Eq. A29) (green) and variation of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ for a typical 11 m s⁻¹ wind from 140^∘ (dashed black). The peak of the $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}{G}^{\mathrm{2}}$ product (red) is shifted with respect to the peak of G² by $\mathit{\delta }\mathit{\phi }\simeq -{\mathrm{0.81}}^{\circ }$.

As an example, Fig. 2 shows the variations of the two-way antenna radiation diagram G², its Gaussian approximation and the $G^{\mathrm{2}}\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ (see Eq. A11) product as a function of azimuth at a 12^∘ incidence angle for a northward-looking KuROS antenna (${\mathit{\phi }}_{\mathrm{b}}={\mathrm{0}}^{\circ }$), using σ₀ data from the Drift4SKIM campaign on 22 November 2018. The effect of the wind-induced azimuthal gradient of σ₀ is to shift the effective radiation diagram towards the brighter upwind and downwind directions, with an apparent pointing azimuth φ_a. The shift induced in this case is $\mathit{\delta }\mathit{\phi }={\mathit{\phi }}_{\mathrm{a}}-{\mathit{\phi }}_{\mathrm{b}}=-{\mathrm{0.81}}^{\circ }=-\mathrm{15}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{rad}$. For comparison (see Sect. 2.3 and Table 2), the pointing accuracy required to achieve a 15 cm s⁻¹ error on the horizontal current in the airborne configuration is 1.2 mrad.

Here, it is important to note that KuROS was not specifically designed for this experiment but primarily as a calibration and validation instrument for the CFOSAT mission, which required a broad radiation diagram. Though the analysis of the KuROS data helped uncover many interesting effects relevant to Doppler observations of the sea surface, its design was not fully appropriate to validate the inversion of the geophysical velocities, for which the pencil-beam antenna diagram of KaRADOC was better suited.

Figure 3 (a) Example of azimuthal variation of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ at a 12^∘ incidence angle, corresponding to the 22 November case (11 m s⁻¹ wind from 140^∘) discussed in Sect. 3, and (b) associated spurious velocity U_AGD as a function of look azimuth φ_b in the case of a port-looking antenna mounted on a platform in constant-altitude flight at 120 m s⁻¹. For the KuROS case, the green line shows the result of the approximation in Eq. (A35), and the black line shows the result of the full azimuthal integration in Eq. (A24). The blue line represents the result of Eq. (A35) for KaRADOC using the same $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ as in the Ku-band case.

Figure 3a shows a typical example of the azimuthal variation of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ at a 12^∘ incidence angle for the Ku band. As expected for near-nadir measurements (Chapron et al., 2002; Munk, 2008; Chu et al., 2012), the NRCS is largest in the downwind look direction ($\mathit{\phi }={\mathrm{320}}^{\circ }$), has a secondary peak in the upwind direction and is weakest in the crosswind look directions. Figure 3b shows the corresponding U_AGD contribution for the KuROS and KaRADOC cases using an aircraft velocity $V_{\mathrm{p}}=\mathrm{120}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and the Ku-band NRCS fit for both instruments (this is a reasonable assumption for order-of-magnitude estimates). As detailed in Appendix A, Eqs. (5) and (6) only apply for a narrow beam when projected on the ground, which is not a very good approximation for the KuROS case, even at a 12^∘ incidence angle. As shown in Fig. 2, the Gaussian approximation for the antenna diagram as a function of φ gives a distribution that is too narrow and does not properly take into account the azimuthal integration, leading to an overestimation of U_AGD. It is clear, however, that even the more exact Eq. (A24) gives very large correction magnitudes, in excess of 1.2 m s⁻¹ in some azimuth ranges.

Because the azimuth gradient U_AGD contribution to the observed DV is proportional to $V_{\mathrm{p}}{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}$, this effect is much larger (and correcting it is correspondingly more demanding in terms of antenna characterization) for KuROS than for KaRADOC or DopplerScatt (Rodríguez et al., 2018) thanks to their narrow azimuthal beam aperture. Another remark is that the approximate expression in Eq. (A35), though it gives the appropriate dependency of U_AGD with respect to look azimuth, tends to overpredict its magnitude as the widening associated with the ground projection saturates for broad beams.

Although the relative variations ${\partial }_{\mathit{\phi }}{\mathit{\sigma }}_{\mathrm{0}}/{\mathit{\sigma }}_{\mathrm{0}}$ are larger for larger incidence angles, this is more than compensated for by the 1∕sin ²θ reduction in azimuthal diversity across the footprint. This effect can thus be neglected for much higher incidence angles (Rodríguez et al., 2018).

2.2 Geophysical velocity U_GD: waves and current Doppler velocities

The geophysical part of the DFS measured by a microwave radar over the ocean, using both along-track interferometry and Doppler centroid techniques, emerges from the average over the instrument field of view (FOV) of the backscatter-weighted line-of-sight projection of the surface velocity, as illustrated in Fig. 4.

In the well-understood case of decametric electromagnetic waves interacting with the sea surface at grazing incidence, the interaction is dominated by the Bragg coherent backscattering mechanism (Crombie, 1955), in which the backscattered field reflects the properties (amplitude, phase speed) of a very finely selected component of the sea state, namely that whose wave vector is precisely equal to the so-called Ewald vector, the difference between the wave vectors of the scattered and incident electromagnetic waves. Exploiting the deviation of the phase speed of this sea state component from its theoretical value is the principle of the HF radars operationally used to measure ocean TSCV in coastal areas (Barrick et al., 1974; Stewart and Joy, 1974).

In the case of the near-nadir interaction of microwaves with the sea surface, which is the configuration considered for SKIM and used by the AirSWOT, KuROS and KaRADOC airborne instruments, this mental picture must be adapted: the Bragg scattering mechanism is not dominant, and the main contribution comes from quasi-specular reflections on those facets of the sea surface which are normal to the Ewald vector. The backscattering cross section of the sea surface and DFS in this case do not depend on the properties of a single Fourier component of the sea state but on the probability density function of the sea surface slope, which is a complex functional of its entire directional spectrum.

Figure 4 Schematic of (a) wave and (b) wave and current contributions to Doppler velocities at the scale of elementary facets. These small-scale processes are averaged over the radar field of view, and a mean velocity signal emerges due to the correlation of surface brightness and velocities in the wave field.

As discussed by Nouguier et al. (2018), who applied it to the analysis of AirSWOT NRCS and DFS data collected during the Gulf of Mexico LASER experiment in 2016, the theoretical framework appropriate for this configuration is the Kirchhoff approximation (Beckmann and Spizzichino, 1987). In this approximation, the geophysical DFS ω_GD can be expressed as

$$\begin{array}{}\text{(7)}& {\mathit{\omega }}_{\mathrm{GD}}=-i{\left.{\displaystyle \frac{{\partial }_{\mathit{\tau }}C}{C\left(\mathit{\tau }\right)}}\right|}_{\mathit{\tau }=\mathrm{0}},\end{array}$$

where C(τ) is the temporal covariance function of the ensemble-averaged electromagnetic field backscattered in the direction of the radar.

Assuming Gaussian statistics for the sea surface, introducing $\mathit{\rho }(\mathit{\xi },\mathit{\tau })=〈\mathit{\eta }(\mathit{x}+\mathit{\xi },t+\mathit{\tau })\mathit{\eta }(\mathit{x},t)〉$, the space–time covariance function of the sea surface elevation, and ${\mathit{Q}}_{\mathbf{H}}=-\frac{\mathrm{4}\mathit{\pi }}{\mathit{\lambda }}\sin \left(\mathit{\theta }\right){\mathit{e}}_{\mathit{\phi }}$ and $Q_z=\frac{\mathrm{4}\mathit{\pi }}{\mathit{\lambda }}\cos \left(\mathit{\theta }\right)$, the horizontal and vertical components of the Ewald vector (with λ the radar wavelength), one obtains C(τ) and ∂_τC as

$$\begin{array}{}\text{(8)}& {\displaystyle }C\left(\mathit{\tau }\right)=\int e^{i{\mathit{Q}}_{\mathbf{H}}\cdot \mathit{\xi }}\left[e^{Q_z^{\mathrm{2}}\left(\mathit{\rho }(\mathit{\xi },\mathit{\tau })-\mathit{\rho }(\mathrm{0},\mathrm{0})\right)}-e^{-Q_z^{\mathrm{2}}\mathit{\rho }(\mathrm{0},\mathrm{0})}\right]\mathrm{d}\mathit{\xi },\text{(9)}& {\displaystyle }{\partial }_{\mathit{\tau }}C=Q_z^{\mathrm{2}}\int {\partial }_{\mathit{\tau }}\mathit{\rho }(\mathit{\xi },\mathit{\tau })e^{i{\mathit{Q}}_{\mathbf{H}}\cdot \mathit{\xi }}{e}^{Q_z^{\mathrm{2}}\left(\mathit{\rho }(\mathit{\xi },\mathit{\tau })-\mathit{\rho }(\mathrm{0},\mathrm{0})\right)}\mathrm{d}\mathit{\xi }.\end{array}$$

The clear upwind–downwind asymmetry of σ₀ observed in the Drift4SKIM radar observations (see Fig. 10) shows that the Gaussian assumption, which is unable to describe such skewness-related effects, is clearly questionable. It is, however, the only practical option, as going further would require prescriptions for the higher-order statistics of the sea surface, which are at present not available.

The occurrence of ρ as the argument of an exponential in these integrals renders further analytical progress difficult (see, however, Nouguier et al., 2011). Approximate expressions can, however, be obtained by performing a Taylor expansion of ρ in the neighborhood of the origin. This results in a Gaussian approximation of the integrand. The integrals can be readily evaluated, yielding (denoting by “⋅” the usual matrix product)

$$\begin{array}{}\text{(10)}& {\mathit{\omega }}_{\mathrm{GD}}\simeq -{{\mathit{Q}}_{\mathbf{H}}}^T\cdot {\left[{\mathrm{\nabla }}_{\mathit{\xi }\mathit{\xi }}\mathit{\rho }\right]}^{-\mathrm{1}}\cdot {\partial }_{\mathit{\tau }}{\mathbf{\nabla }}_{\mathit{\xi }}\mathit{\rho }.\end{array}$$

The derivatives of ρ are taken at ξ=0 and τ=0, and they can be expressed as moments of the directional sea state spectrum S_d(k) as

$$\begin{array}{}\text{(11)}& {\partial }_{\mathit{\tau }}{\mathbf{\nabla }}_{\mathit{\xi }}\mathit{\rho }=\mathit{m}\mathit{s}\mathit{v},{\mathrm{\nabla }}_{\mathit{\xi }\mathit{\xi }}\mathit{\rho }=-\mathbf{Mss},\end{array}$$

where, in the notation of Nouguier et al. (2018), msv stands for the mean slope velocity and Mss for the mean square slope matrix,

$$\begin{array}{}\text{(12)}& \mathit{m}\mathit{s}\mathit{v}=\left[\begin{array}{c}{\mathrm{mss}}_{xt}\\ {\mathrm{mss}}_{yt}\end{array}\right],\mathbf{Mss}=\left[\begin{array}{cc}{\mathrm{mss}}_{xx}& {\mathrm{mss}}_{xy}\\ {\mathrm{mss}}_{yx}& {\mathrm{mss}}_{yy}\end{array}\right],\end{array}$$

with

$$\begin{array}{}\text{(13)}& {\mathrm{mss}}_{x^{\mathit{\alpha }}{y}^{\mathit{\beta }}{t}^{\mathit{\gamma }}}=\mathrm{2}\underset{{\mathbb{R}}^{\mathrm{2}}}{\int }{k}_x^{\mathit{\alpha }}{k}_y^{\mathit{\beta }}{\mathit{\omega }}^{\mathit{\gamma }}{S}_{\mathrm{d}}\left(\mathit{k}\right)\mathrm{d}\mathit{k}.\end{array}$$

The surface current enters through its effect on the dispersion relation of surface waves ω(k). In the presence of a vertically homogeneous current U (a detailed discussion of the effect of shear can be found in Kirby and Chen, 1989),

$$\begin{array}{}\text{(14)}& \mathit{\omega }\left(\mathit{k}\right)=\mathit{k}\cdot \mathit{U}+{\mathit{\omega }}_{\mathrm{0}}\left(\right|\mathit{k}\left|\right),\end{array}$$

where

$$\begin{array}{}\text{(15)}& {\mathit{\omega }}_{\mathrm{0}}\left(\right|\mathit{k}\left|\right)=\sqrt{g\left|\mathit{k}\right|\left(\mathrm{1}+|\mathit{k}{|}^{\mathrm{2}}/{\mathit{\kappa }}_M^{\mathrm{2}}\right)}\end{array}$$

is the dispersion relation of gravity–capillary waves in deep water, with k the wave vector and ${\mathit{\kappa }}_M=\mathrm{363.2}\mathrm{rad}{\mathrm{m}}^{-\mathrm{1}}$ the wavenumber corresponding to the gravity–capillary regime transition. Introducing this expression in Eq. (10), and defining msv₀ as the spectral moment obtained using the dispersion relation of Eq. (15) in Eq. (13), one obtains the approximate DFS as

$$\begin{array}{}\text{(16)}& {\mathit{\omega }}_{\mathrm{GD}}={{\mathit{Q}}_{\mathbf{H}}}^T\cdot \left[{\mathbf{Mss}}^{-\mathrm{1}}\cdot \mathit{m}\mathit{s}{\mathit{v}}_{\mathrm{0}}+\mathit{U}\right]\end{array}$$

and the corresponding V_GD as

$$\begin{array}{}\text{(17)}& V_{\mathrm{GD}}=-\sin \left(\mathit{\theta }\right){\mathit{e}}_{\mathit{\phi }}\cdot \left[{\mathbf{Mss}}^{-\mathrm{1}}\cdot \mathit{m}\mathit{s}{\mathit{v}}_{\mathrm{0}}+\mathit{U}\right].\end{array}$$

While clearly oversimplified (it is, for instance, independent of the electromagnetic wavelength, which is known to have a significant influence on σ₀), this expression has a definite pedagogical interest, as it allows one to distinguish a number of interesting features.

• The raw velocity projection V_GD accessible to Doppler radar instruments is composed of a “genuine” current Doppler (CD) component V_CD, equal to the projection of the TSCV along the radar line of sight, plus a wave Doppler (WD) component V_WD induced by the natural motion of the sea surface.

• This V_WD component involves sea surface statistics of two different natures: the mean slope velocity vector msv₀ and the mean squared slope matrix Mss. To this order of approximation it can be seen as the projection along the radar line of sight of the constant vector ${\mathbf{Mss}}^{-\mathrm{1}}\cdot \mathit{m}\mathit{s}{\mathit{v}}_{\mathrm{0}}$. In the rest of this article, M_WD denotes the norm of this vector.

• As noted in Nouguier et al. (2018), msv₀ is equal to one-half the surface Stokes drift velocity of deepwater waves ${\mathit{U}}_S^{\mathrm{\infty }}$. As noted in Nouguier et al. (2016), the effective mean squared slope matrix Mss components (mss_shape), accounting for the electromagnetic filtering effect and part of the non-Gaussianity of the sea surface statistics, can be obtained from the derivatives of σ₀ as a function of incidence angle for different azimuths.

• In simple cases represented by parametric spectral forms such as the Elfouhaily et al. (1997) spectrum used in this work, msv₀ and the eigenvectors of Mss are aligned with the downwind direction, and the $V_{\mathrm{WD}}=-G_{\mathrm{D}}\sin \left(\mathit{\theta }\right){\mathit{e}}_{\mathit{\phi }}\cdot {\mathit{U}}_S^{\mathrm{\infty }}$ relation proposed in Chapron et al. (2005) is recovered, with $G_{\mathrm{D}}=\frac{\mathrm{1}}{{\mathrm{mss}}_{\mathrm{shape}}}$.

• Both these statistics are, however, known to be influenced by waves at all scales. The asymptotic behaviors of the weighting factors as functions of the surface wave wavenumber in the gravity wave range are k^3∕2 and k² for the msv₀ and Mss terms, respectively, while the parametric spectrum of Elfouhaily et al. (1997), used in this work, decays as k⁻³, leading to a logarithmic divergence for the Mss components and a slow convergence of the msv₀ components at high wavenumbers.

• The Mss components are sensitive to the detailed shape of the spectrum up until the short capillary wave roll-off or to the electromagnetic cutoff, whichever is reached first.

• Estimating these terms requires knowledge of all the components of the sea state: the long gravity wave range can be measured (either in situ, as during the Drift4SKIM campaign, or using the radar measurements themselves as intended in the SKIM context), but the high-wavenumber range cannot be neglected, and its effect must be accounted for, for instance, through the use of a parametric spectral form.

• The msv₀ vector appears as a multiplicative factor, to which the inverse of the Mss matrix is applied. These terms thus have opposing influences on the final result: modifications of the sea spectrum, which tend to increase the weight of small-scale components, increase the mean slope velocity but also, and rather more, the mean squared slope by which it is divided. A certain degree of stability of the end result is thus likely.

• On a similarly reassuring note, whereas the low-wavenumber part of the spectrum is affected by swell systems of remote origin that have arbitrary orientations, the short waves represented by the parametric tail of the spectrum are known to be aligned with the wind direction and to depend on local variables only (wind strength and direction, fetch).

Figure 5 Computed variability of the Stokes drift velocity, the diffraction-effective mean square slope mss_shape and the wave Doppler velocity magnitude M_WD. (a) 2010–2017 statistics of Stokes drift magnitude at Ocean Station Papa, computed using the buoy wind speed data and wave data from the nearby WMO buoy 46 246, maintained by the University of Washington (Thomson et al., 2013). (b) mss_shape estimated from GPM satellite backscatter using modeled co-located wind speed and wave height; reproduced from Nouguier et al. (2018). (c, d) Statistics of the Ka- and Ku-band M_WD, computed using the theoretical model of Nouguier et al. (2011) for ocean wave spectra modeled over the global ocean using the WAVEWATCH III model (Stopa et al., 2016) and plotted as a function of the wind speed. The colored curves show the median value for different classes of wave height for a given wind speed; each curve is separated by 0.5 m in panels (a) and (b) and by 1 m in panels (c) and (d). In (c) and (d), the gray shading represents the histogram of the computed M_WD values in the global simulation.

Figure 5 gives orders of magnitude for the natural range of variability of the different factors thus isolated. Figure 5a shows the variability of the Stokes drift velocity estimated following Kenyon (1969) and Ardhuin et al. (2009) using wind and directional wave measurements collected from 2010 to 2017 at Ocean Station Papa. Even though U_S is highly correlated with the wind speed with a Pearson's linear correlation coefficient of 0.85, a strong dependence on the long-wavelength part of the spectrum, for which Hs is a proxy, definitely has to be accounted for.

Figure 5b, taken from Nouguier et al. (2016), instead shows the dependence on wind speed and Hs of Ku- and Ka-band effective mean squared slope mss_shape retrieved from the GPM (Global Precipitation Mission) satellite measurements. The variability is even more strongly dominated by the dependence on wind speed, the variability due to the long-wavelength part of the spectrum being much smaller. These measurements very clearly show the filtering effect of the electromagnetic wavelength and are a clear warning that Eq. (17), suggestive though it is, should be considered with caution.

Finally, Fig. 5c and d show the magnitude of the horizontal U_WD component as a function of wind speed and Hs, estimated by numerically evaluating the integrals of Eqs. (8) and (9) for C(τ) and ∂_τC using the numerical tools of Nouguier et al. (2011) on the basis of long-wavelength spectra extracted from global runs of the WAVEWATCH III model (Stopa et al., 2016), completed in the high-wavenumber range by Elfouhaily et al. (1997) spectral tails. The shading in the background represents the histogram of the different (wind speed, M_WD) pairs. As could be hoped for, the opposing influences of the wind speed on msv₀ and Mss tend to counteract each other, greatly reducing the range of variability of M_WD with wind. This effect appears stronger in the Ka rather than the Ku band, possibly due to the saturation of the Ku-band mss_shape at high winds. These figures show a strong remaining impact of the long-wavelength waves, which clearly must be accounted for. As wind speed and significant height are highly correlated variables, the frequently encountered situations fall in a quite narrow interval M_WD≃C₀, with C₀≃2.6 and C₀≃2.2 m s⁻¹ in the Ku and Ka band, respectively. In other words, most of the variability of U_WD is controlled by the directionality effect, and the magnitude M_WD is a weakly varying function of the wind, the wave age and the presence of swell (see also Yurovsky et al., 2019; Ardhuin et al., 2019b).

A final remark is that, though these general patterns can probably be assumed to be robust, the precise numerical values depend on the parametric spectral shapes which have been used to fill the high-wavenumber range of the spectra. Changing, for instance, the high-wavenumber azimuthal spreading functions, which are for the moment not very well constrained observationally, has different impacts on the msv₀ and Mss terms and can thus be expected to marginally change the numbers.

2.3 Error budget

Considering the errors on the different terms to be independent, developing Eq. (1) allows one to derive the error variance of Doppler radar measurements of the TSCV as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathrm{Var}\left(\mathit{\delta }{U}_{\mathrm{CD}}\right)=& {\left({\displaystyle \frac{U_{\mathrm{CD}}}{\tan \mathit{\theta }}}\right)}^{\mathrm{2}}\mathrm{Var}\left(\mathit{\delta }\mathit{\theta }\right)+{\displaystyle \frac{\mathrm{Var}\left(\mathit{\delta }{V}_{\mathrm{LOS}}\right)}{{\sin }^{\mathrm{2}}\mathit{\theta }}}\\ & +{\displaystyle \frac{\mathrm{Var}\left(\mathit{\delta }{V}_{\mathrm{NG}}\right)}{{\sin }^{\mathrm{2}}\mathit{\theta }}}+\mathrm{Var}\left(\mathit{\delta }{U}_{\mathrm{WD}}\right).\end{array}\end{array}$$

As a first step, four contributions to the uncertainty on U_CD can thus be isolated, with different origins.

• A first part corresponds merely to the error caused by imperfect knowledge of the projection angle between the TSCV and the line of sight. Its order of magnitude is controlled by the TSCV, and it is thus negligible with respect to similar terms that involve the platform velocity.

• The second term corresponds to the random error in the DFS measurements and subsumes the dependence on the signal-to-noise ratio, antenna beamwidth, orientation of the boresight with respect to the platform velocity vector and algorithmic choices. A very thorough analysis of this term can be found in Rodriguez (2018). The standard deviation of the raw DV signal carries over to the end result, multiplied by a 1∕sin θ factor of the order of 5 for $\mathit{\theta }={\mathrm{12}}^{\circ }$.

• The third term corresponds to the error caused by mismatches between the actual platform motion contribution to V_LOS and the estimate computed from the ancillary sensors. The order of magnitude of this term is set by the (very large) platform velocity. It is by far the largest.

• The fourth and final term corresponds to the uncertainty on the wave Doppler removal stage. Errors in the U_WD model carry directly over to the U_CD estimates.

The third term dominates the overall error budget and must be further analyzed. It is convenient for that purpose to start from Eq. (A23), which gives the expression of the beam direction vector, and use the platform velocity components in the local (northward–eastward–downward, NED) frame at the observation point. Neglecting terms involving the vertical velocity of the platform and introducing the difference between the boresight and flight track azimuths $\mathit{\psi }={\mathit{\phi }}_{\mathrm{b}}-{\mathit{\phi }}_{\mathrm{t}}$, one obtains the consolidated error budget as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\mathrm{Var}\left(\mathit{\delta }{U}_{\mathrm{CD}}\right)& ={\displaystyle \frac{\mathrm{Var}\left(\mathit{\delta }{V}_{\mathrm{LOS}}\right)}{{\sin }^{\mathrm{2}}\mathit{\theta }}}+\mathrm{Var}\left(\mathit{\delta }{U}_{\mathrm{WD}}\right)+{\displaystyle \frac{\mathrm{Var}\left(\mathit{\delta }{V}_D\right)}{{\tan }^{\mathrm{2}}\left(\mathit{\theta }\right)}}\\ & +\mathrm{Var}\left(\mathit{\delta }{V}_N\right){\cos }^{\mathrm{2}}\left({\mathit{\varphi }}_{\text{b}}\right)+\mathrm{Var}\left(\mathit{\delta }{V}_E\right){\sin }^{\mathrm{2}}\left({\mathit{\varphi }}_{\text{b}}\right)\\ & +V_{\mathrm{P}}^{\mathrm{2}}\left[{\displaystyle \frac{{\cos }^{\mathrm{2}}\left(\mathit{\psi }\right)}{{\tan }^{\mathrm{2}}\left(\mathit{\theta }\right)}}\mathrm{Var}\left(\mathit{\delta }\mathit{\theta }\right)+{\sin }^{\mathrm{2}}\left(\mathit{\psi }\right)\mathrm{Var}\left(\mathit{\delta }\mathit{\phi }\right)\right].\end{array}\end{array}$$

This equation summarizes the dependence of the overall U_CD error on the errors introduced by the Doppler measurements, the U_WD model, the individual platform velocity components, and the incidence angle and azimuth mispointing errors.

Table 2 Standard deviations of the different error terms in Eq. (19) necessary to achieve a 0.40 m s⁻¹ standard deviation for U_CD.

As an illustration, Table 2 summarizes the requirements that have to be met to keep the standard deviation of each of the seven terms below 0.15 m s⁻¹, ensuring a 0.4 m s⁻¹ standard deviation for U_CD. The requirement for θ is translated to the corresponding altitude-tracking accuracy requirement for the KuROS and SKIM configurations. The requirements for linear velocity components are stringent but can be reached using current-day technology. The requirement for altitude accuracy is easily within the specifications of the SKIM nadir altimeter payload but definitely out of reach of KuROS. The KuROS data could, however, be analyzed in the cross-track-looking configurations for which this requirement does not apply. The requirement for azimuthal pointing accuracy is by far the most stringent. In the airborne case, it is met for the antenna boresight by the plane IMU (inertial measurement unit), allowing a straightforward analysis of the KaRADOC data. In the KuROS case, however, it is exceeded by a factor of 10 by the mispointing induced by the azimuthal gradients of sea surface, which required the development of a specific data correction procedure. Finally, in the spaceborne case, it seems only achievable using a combination of high-end inertial measurements and data-driven analysis techniques.

3 Campaign overview

This section provides a general overview of the campaign. The location, timing and overall organization are described in Sect. 3.1, the environmental conditions encountered during the campaign are described in Sect. 3.2, and the two main instruments, the KuROS and KaRADOC airborne radars, are described in Sect. 3.3 and 3.4, respectively.

3.1 Campaign organization

The Drift4SKIM experiment differs from previous airborne Doppler radar campaigns (Martin et al., 2016; Rodríguez et al., 2018) in two important respects: in order to observe the effect of wave development on the geophysical Doppler velocity U_GD, it was performed in a midlatitude, eastern basin oceanic environment open to offshore swells. Also, given the campaign objectives of demonstrating the sensitivity of airborne radar Doppler measurements to the geophysical contributions of currents and waves, it comprised an extensive in situ component designed to have commonly accepted reference measurements for these parameters.

Fieldwork was performed in two areas (denoted by square boxes in Fig. 6) named the “offshore” area, centered on the Trèfle buoy (see below), and the “Keller Race” area to the north of the island of Ushant. Both locations are in the range of coverage of a two-site WERA (Gurgel et al., 1999) high-frequency radar system, operated by Service Hydrographique et Oceanographique de la Marine (Shom) and already used for several studies, in particular related to wave–current interactions (Ardhuin et al., 2009, 2012; Guimaraes et al., 2018).

Keller Race is an area with very strong horizontal gradients of the current (Sentchev et al., 2013). Although it is easy to show a strong effect of the current on the measured DFS, the spatial variability of the sea state is difficult to measure in situ, introducing uncertainties when combining U_CD+U_WD in a forward model or using U_WD estimates when retrieving U_CD from the measured U_GD. The offshore area, on the other hand, was chosen for its spatial uniformity, being located far enough from the islands and with a near-uniform depth of 110 m. Only airborne data acquired over the offshore area are presented in this paper.

The week around spring tides in November 2018 was selected in order to allow for a wide range of current speeds (Fig. 7a).

Figure 6 Location of the measurement campaign and in situ assets, including a map of the KaRADOC measurements of the geophysical Doppler velocity U_GD acquired on 22 November 2018.

The KuROS and KaRADOC radars were installed on an ATR-42 plane operated by the French institutional scientific flight facility, SAFIRE, which is equipped with an AIRINS™ GNSS-FOG INS providing position, pitch, roll and heading information with stated tolerances of a few centimeters, 0.005, 0.005 and 0.01^∘, respectively.

Ground truth measurements comprised two permanent operational systems: the HF radar system mentioned previously, with an expected depth of measurement around 1 m (Stewart and Joy, 1974), and the Pierres Noires (WMO no. 62069) wave-measuring buoy. Dedicated instrumentation was also deployed for the campaign.

• The Trèfle buoy was moored at 5^∘15^′ W, 48^∘15^′ N at the center of the offshore area. This buoy monitored the surface current (Sutherland et al., 2016) and provided directional wave spectra (Fig. 8).

• Several types of drifting buoys, including CARTHE drifters (Novelli et al., 2017) drogued around 40 cm, SVP drifters (Niiler and Paduan, 1995) drogued at 15 m and Spotter wave-measuring buoys (Raghukumar et al., 2019), were deployed in the measurement areas.

• The R/V Thalia worked in the offshore area, providing continuous underway measurements of meteorological parameters using a Météo-France BATOS operational system comprising a Vaisala WXT series sonic anemometer located approximately 10 m above the sea surface. The ship also carried a SBE21 thermosalinograph.

In the summer, the so-called Ushant tidal front has a strong influence on the surface currents, as well as hydrographic (Le Boyer et al., 2009) and atmospheric (Redelsperger et al., 2019) conditions in the offshore area. This seasonal feature typically disappears in October, and conductivity–temperature–depth (CTD) casts were performed from R/V Thalia to confirm that it had indeed vanished when the campaign took place. The water column was found to be very well mixed, with surface-to-bottom potential density anomalies being smaller than 0.002 kg m⁻³. The spatial homogeneity was also checked using the ship thermosalinograph and an infrared camera mounted on a Piper PA-23 plane, which surveyed the offshore area in a “lawn-mowing” pattern, flying under the clouds at an altitude of 500 to 1000 m. While small-scale surface features were observed on calm days, it is clear that no density-associated mesoscale structures were present.

Figure 7 Time series at the location of the Trèfle buoy (5^∘15^′ W, 48^∘15^′ N) in the offshore zone of (a) ocean surface current speed from the MARS2D numerical model run at LOPS (Lazure and Dumas, 2008). (b) Wind speed (black) and direction (blue) from the AROME regional operational model run by Météo-France. (c) Total (blue) and swell (black) significant wave height and wave peak frequency (red) from the WAVEWATCH III numerical wave model run at LOPS (Roland and Ardhuin, 2014). The four time periods shaded in gray correspond to the times of fixed-antenna KuROS measurements. The corresponding observed environmental parameters are detailed in Table 3.

The airborne radar measurements geometry over the offshore area consisted of relatively long (12 km) and straight tracks with different aircraft headings, forming a star pattern, as for the 22 November 2018 flight shown in Fig. 6. Tracks were flown every 12, 22.5 or 45^∘ in azimuth, depending on flight duration constraints. The KaRADOC antenna was fixed relative to the aircraft and looking to port, while the KuROS antenna could either be fixed in the up-track or port cross-track directions, or it could rotate in the clockwise sense relative to the flight track. The KuROS Doppler data presented in this paper were acquired in the port-looking configuration.

Figure 8 Directional wave spectra E(f_r,θ), as functions of the relative wave frequency f_r and incoming wave azimuth θ, estimated from the motions of the Trèfle buoy on 22 November at 13:00 UTC and Spotter buoy number 10 on 24 November at 12:00 UTC. The measured directional moments were transformed with the maximum entropy method (Lygre and Krogstad, 1986) and Doppler-shifted with $f_{\mathrm{r}}=f-\mathit{k}\cdot U/\left(\mathrm{2}\mathit{\pi }\right)$ for the moored Trèfle buoy. The red and blue arrows represent the AROME wind and MARS2D surface current vectors directions, respectively.

Table 3 Surface current velocity, Stokes drift and wind speed measured or estimated near position 48^∘15^′ N, 5^∘15^′ W. For each table entry, the parenthesized pair contains the (eastward, northward) components of the vector (cm s⁻¹) for current or Stokes drift estimates and (m s⁻¹) for wind speed estimates. Please note that the Stokes drift is only integrated up to 0.5 Hz. Stokes drift buoy data correspond to the Trèfle buoy for 22 November and Spotter buoy number 10 for 24 November.

3.2 Geophysical conditions

A wide range of geophysical conditions were encountered during the 1-week-long campaign. Four flights were performed over the offshore area on 21 November from 13:50 to 15:50 UTC, on 22 November from 12:15 to 15:00 UTC, on 24 November from 11:20 to 13:20 UTC, and finally on 26 November from 09:40 to 11:00 UTC. In this paper, we focus on data acquired on 22 and 24 November as the geophysical conditions were interesting and complementary (see below), and data were acquired with the largest azimuth diversity on these two days.

The 22 November flight took place at the end of a steady southeasterly wind episode (13 m s⁻¹ from 140^∘). The 24 November flight, in contrast, took place during a steady weak southwesterly wind period (5 m s⁻¹ from 225^∘) (Fig. 7b).

The wave height during the campaign was dominated by the presence of two swell systems from North Atlantic remote storms. The swell height decreased from 2.5 m on 21 November to 0.9 m on 24 November, with a peak frequency increasing from 0.07 to 0.1 Hz and a mean direction gradually veering from northwest to west. This swell has a small contribution to the Stokes drift of the order of 10 % of the wind-sea contribution on 22 November.

The main environmental conditions at the time of these star-pattern flights are summarized in Table 3.

3.3 KuROS instrument

KuROS is a Ku-band (13.5 GHz) pulse-pair Doppler radar with a dual antennae system and azimuthal scanning possibility, which was developed in the framework of the CFOSAT prelaunch studies. Of the two antennas, the low-incidence (LI) antenna is nominally centered on a 14^∘ incidence angle, while the medium-incidence (MI) antenna is nominally centered on a 40^∘ incidence angle. Only the LI antenna, which was the more relevant for SKIM, was used during the campaign. This antenna uses a horizontal transmit–receive polarization (HH) polarization. A comprehensive description of the system can be found in Caudal et al. (2014). A new antenna was used for the Drift4SKIM campaign, with characteristics given in Table 1.

The radar transmits a frequency-modulated pulse (chirp) with a 100 MHz bandwidth, achieving a 1.5 m range resolution and an effective ground-projected resolution of approximately 7 m (at 12^∘). The one-way 3 dB footprint in azimuth is 580 m wide at 12^∘ and 3000 m of flight altitude. The pulse repetition frequency (PRF = 1 ∕ PRI) depends on the altitude and is 23 kHz when the aircraft flies at 3000 m. The ambiguity of the Doppler velocity measurement (see Sect. A1.4 in the Appendix) is about 126 m s⁻¹, which is much larger than expected from the measurements (below aircraft speed of 120 m s⁻¹). In order to reduce the thermal noise contribution, the range-resolved pulse-pair signal is coherently averaged in the instrument over 1 ms, corresponding to 22 pulse pairs per instrument sample. For the purpose of this article, this was further coherently averaged per blocks of 15 samples.

As discussed in Appendix A, accuracy requirements for observation geometry are much less stringent for cross-track than for up-track and down-track Doppler velocity observations. The Doppler velocity data discussed in this article were all collected with the KuROS antenna in the port-looking orientation. This configuration also ensures an overlap with the KaRADOC footprint.

3.4 KaRADOC instrument

The Ka-band RADar for Ocean Current (KaRADOC) airborne radar sensor was developed for the Drift4SKIM campaign. KaRADOC is derived from the Still WAter Low Incidence Scattering (SWALIS) instrument, developed for the measurement of the NRCS of inland water surfaces in the Ka band. Further details on the system are given in Appendix C.

KaRADOC was mounted under the ATR-42 aircraft in a port-looking configuration. The two-way 3 dB footprint from 3000 m of altitude over a flat sea surface is an ellipse with diameters 45 and 60 m in the cross-track and along-track directions, respectively. The antenna is of the slotted-waveguide type and allows steering of the beam in elevation (incidence angle) by varying the working frequency. Data were acquired at different incidence angles from 6 to 14^∘, corresponding to a range of frequencies from 32.5 to 38.2 GHz. This article focuses on the observations collected at $\mathit{\theta }={\mathrm{12}}^{\circ }$ at 33.7 GHz.

KaRADOC does not implement a range-resolution scheme: the transmitted pulses last several microseconds, and the whole FOV is illuminated simultaneously. The demodulated return signal is sampled at 15 MHz and archived. It is essentially constant while the electromagnetic wave is actually interacting with the sea surface. The useful signal segment is selected and its average is computed in order to reduce the thermal noise contribution, yielding one complex amplitude for each pulse. Several hundred pulses are sent at 4 kHz PRF for each burst of measurements, with a burst repetition frequency of the order of 5 to 10 Hz, depending on the number of incidence angles in the scanning sequence. These parameters were varied during the acquisitions. Though they have a strong impact on NRCS and DFS estimate quality, we have found the low-pass-filtered DFS signal to be robust.

The pulse-pair complex signal is averaged for each burst in order to reduce the effect of coherent speckle. One complex pulse-pair sample is thus obtained per burst. Even at the lowest burst repeat frequency of 5 Hz, the plane moves by less than a third of the FOV along-track extension between bursts.

The impact of the acquisition parameters on the KaRADOC measurement normalization is not yet fully understood, and the NRCS measurements could not be exploited in the scope of this study. The noise-filtered DFS measurements are, however, not affected by these normalization changes and are valid.

Figure 9 (a, c) Mosaics of KuROS backscattering intensity and Doppler velocity data acquired on 22 November with fixed port-looking antenna. (b) Overlay of the Doppler velocity and backscattering intensity. Flight tracks are numbered 1 to 6; black arrows indicate the flight direction, and white arrows point in the radar look direction. The long dashed lines represent the apparent direction of swell crests. Panels (d)–(i) show close-up views of selected tracks from (a) and (c). The tracks shown in (d) and (e) are out of the frame of (a) and (c). Panels (f)–(i) show close-ups of flight tracks 6 and 2. The 570 m scale bar applies to (f)–(i) and corresponds to the along-track 3 dB width of the radar beam at a 12^∘ incidence angle, i.e., near the middle of the swath. The mean trend of σ₀ as a function of θ has been removed from the σ₀ data.

4 Measurements

4.1 KuROS NRCS–DFS imagery

The KuROS NRCS–DFS imagery reveals a host of interesting features, modulations and dependencies. An in-depth analysis of all these processes is clearly out of the scope of this paper and will be the subject of forthcoming contributions from the Drift4SKIM team. This section thus only provides a cursory description of a few segments of σ₀ and DV data collected on 22 November 2018 when the wind speed was approximately 11 m s⁻¹, which are displayed in Fig. 9.

A first remark is that the NRCS is smooth, with a typical modulation depth of 1 dB after removing its mean trend as a function of incidence angle (Fig. 9a). This smoothness is in part due to the large footprint, but it also shows that the radiometric quality of the data and the coherent averaging performed are sufficient to control the thermal noise. Speckle noise is, however, still present, with different statistics depending on the radar look direction and the variable considered (not shown). The cross-track observation geometry leads to the best speckle noise reduction for the NRCS but to the worst-case speckle noise statistics for the DV.

The KuROS data clearly show a modulation in both NRCS and U_GD associated with the northwesterly swell observed by the Trèfle buoy with a peak frequency of 0.07 Hz, corresponding to a wavelength L=320 m (Fig. 8). This is particularly visible on the north–south-oriented flight tracks numbered 5 and 6 in Fig. 9b (see also Fig. 9f–g for a zoom-in on track 6). The apparent swell crest direction (dashed lines in Fig. 9b) differs from the true direction due to the scanning distortion effect (Walsh et al., 1989; Sutherland et al., 2018), as the swell propagates during the measurements at a phase speed of 22 m s⁻¹, while the aircraft moves at 120 m s⁻¹.

The shorter waves measured by the Trèfle buoy (Fig. 8) occupy a wide range of directions from a narrow wind-sea peak from the south at 0.16 Hz (L=60 m) to a broad directional distribution at 0.22 Hz (L=30 m), with a mean direction of 130^∘ and a half-width (spread) of 45^∘, hence covering directions from 85 to 175^∘. These shorter components are present in the data from flight tracks 5 and 6 in the form of very narrow stripes with orientations shown by short dashed lines in Fig. 9b (see also Fig. 9f–g for a zoom-in on track 6). The “long-crested” appearance of the short waves in (d) and (e) is an artifact due to the wavefront-matching observation geometry (Jackson et al., 1985), with all other directions averaged out by the large azimuth width of the radar beam. If purely geophysical, the phase relationship between the DV and NRCS modulations is expected to give the wave propagation direction. For flight track 6 in (f) and (g), the long swell propagates towards the radar, and the brighter slopes (white) correspond to eastward velocities toward the radar (blue). This will be discussed in further detail below. Finally, (h) and (i) exhibit chevron patterns with crests facing both northeast and northwest. Whereas the waves from the southwest are expected to be much longer than those from the southeast, this is not apparent in the KuROS data.

Figure 10 (a, b) Variations as a function of azimuth φ of σ₀ for an incidence angle θ of 6 to 18^∘ on 22 November for the port-looking (12:13–13:38 UTC) and rotating-antenna (13:41–13:58 UTC) flights, respectively. (c, d) Variations of σ₀ normalized by its azimuthal average for the fixed- and rotating-antenna data, respectively.

4.2 Ku-band NRCS

In this section we discuss the dependence of the Ku-band NRCS as a function of azimuth and incidence angle for the 22 and 24 November cases. The fixed-antenna and rotating-antenna data are presented. In order to reduce the dispersion introduced by the short-scale modulating processes discussed above, the data were averaged per 1^∘ incidence angle and azimuth bins. As mentioned before, full tracks are straight and relatively long (12 km), and they view a mainly homogeneous ocean region. For the fixed-antenna observations, azimuthal diversity is obtained by performing tracks in different flight directions, forming a “star” pattern.

The variations of the Ku-band σ₀ are shown for 22 November in Fig. 10. These measurements show the expected modulation of 0.8 to 0.9 dB with azimuth, with a downwind–crosswind contrast that increases with the incidence angle. This contrast is larger for the higher winds on 22 November. The upwind–downwind asymmetry is expected from the behavior of the surface slope probability density function (Chapron et al., 2002; Walsh et al., 2008; Munk, 2008). The exception are the σ₀ values for the flight tracks with a fixed antenna around the azimuths 90 and 270 (Fig. 10a), which have anomalous normalized values between 1 and 1.3 instead of expected values much closer to 1. We have no explanation for this anomaly, which is genuine. No such anomaly was found for the rotating-antenna data collected later on the same day (Fig. 10b).

Discarding these azimuth ranges (shaded in gray in Fig. 10c), the data could be well fitted with a functional form $a_{\mathrm{0}}+a_{\mathrm{1}}\cos (\mathit{\phi }-{\mathit{\phi }}_{\mathit{\sigma },\mathrm{1}})+a_{\mathrm{2}}\cos \left[\mathrm{2}\right(\mathit{\phi }-{\mathit{\phi }}_{\mathit{\sigma },\mathrm{2}}\left)\right]$. As explained in Sect. 2.1, measuring this azimuthal variation is critical for the interpretation of the mean Doppler velocity due to the spurious azimuth gradient contribution. As expected, the fitted directions φ_σ,1 and φ_σ,2 are very close to the wind direction, except for the lowest incidence angles for which the contrast is less than 0.05 dB.

On 24 November, the σ₀ azimuthal contrast was much weaker (Fig. 11) due to the much lower wind speed and was actually not aligned with the wind direction when the measurements were performed.

Figure 11 Same as Fig. 10c but using port-looking antenna data collected on 24 November at 11:22–13:03 UTC.

4.3 Mean Doppler velocity from KaRADOC

We now quantitatively discuss the measured Doppler velocity signal in order to assess the agreement of our theory of the wave-induced contribution U_WD with the measurements. This section is focused on the KaRADOC data, which are easier to interpret than the KuROS data due to the narrower radar beam of the instrument.

We present in Fig. 6 the low-pass-filtered U_GD estimates retrieved from the 12^∘ incidence angle KaRADOC data collected on 22 November between 12:13 and 12:59 UTC.

This representation is misleading, as much of the observed variability is in fact due to the effect of the flight track orientation. For instance, the largest contrast can be observed between the northeastward- and southwestward-directed flight tracks, even though to a first approximation a mere change in observation direction has occurred.

Another representation of the same data is proposed in Fig. 12. In this figure the U_GD data are represented as blue lines shifted to the right of the plane ground track (in black) by an amount proportional to the instantaneous low-pass-filtered U_GD value. This representation removes the trivial effect of observation direction changes and allows subtler effects to be better appreciated. For instance, noise-free observations of a constant vector U_GD would appear as straight lines parallel to the flight tracks, all crossing at the tip of the vector. Deviations from this behavior, such as can be observed in Fig. 12, are indicative of measurement noise, geophysical variability or geophysical phenomena not accounted for by our theory.

For 22 November, 16 flight tracks are available collected from 12:13 to 12:59 UTC, and for 24 November 17 tracks were collected from 11:27 to 13:13 UTC.

Figure 12 Plots of the Ka-band Doppler velocity signal on (a) 22 and (b) 24 November 2018. The flight tracks are marked as thin black lines. For each flight track, a thick blue line shifted to the right of the flight path by an amount proportional to the instantaneous low-pass-filtered Doppler velocity represents the projection of the U_GD vector along the instrument line of sight. At the beginning of each track data were discarded until the plane stabilized. The green arrow represents the maximum likelihood estimate of the U_GD vector using the whole data set. The (almost indistinguishable) red arrow shows the result of the least-squares sinusoidal fits shown in Fig. 13a and b. The 1 standard deviation error ellipse on the maximum likelihood estimate is represented in green.

Overall, the assumption of a constant vector is good to within 0.3 m s⁻¹. It is particularly striking that the three horizontal lines in Fig. 12a are almost perfectly aligned, corresponding to two flight tracks looking into azimuth 0^∘ and one flight into azimuth 180^∘. On 22 November, the largest dispersion is for the 315 and 135^∘ azimuths for which a total of four tracks are available with very different values that are, however, consistent along each track.

Using the average values from the different tracks, we compare the measured Doppler velocity to the forward model given by Eq. (1), with U_WD estimated from the in situ wave buoy data using the tools discussed in Sect. 2.2. The method combines the buoy spectrum up to 0.35 Hz and adds a high-frequency tail based on the Elfouhaily (1997) spectrum, then computes numerically the integrals of Eqs. (8) and (9) to obtain the DFS estimate. The TSCV contribution, U_CD, is taken to be the drift velocity of the nearest CARTHE drifter, which is uniform to within 3 cm s⁻¹ in the offshore area (interactive animations of all deployments and trajectories can be found at https://odl.bzh/eVRHv1TE, last access: 25 October 2020).

Figure 13 KaRADOC Doppler velocity (red circles) for the star-pattern flights on (a) 22 and (b) 24 November. Cosine function fits to the data (red lines). Modeled geophysical Doppler velocity U_GD using the MEM (MLM) estimate of the directional wave spectrum (green and darker green). The modeled U_GD is the sum of the CARTHE drifter velocity U_CD (blue) and the wave Doppler velocity estimated from the measured spectra, U_WD (midnight blue dashes).

Figure 13 shows the measured mean Doppler velocity and standard deviation for each track (the standard deviation is representative of the order of magnitude of the short-scale modulations due to waves, not of the error bar for the mean DV). On 22 November (Fig. 13a), the current vector accounts for less than half of the observed magnitude of U_GD, and it is interesting that the maximum Doppler velocity is from azimuth 147^∘, between the wind direction (128^∘) and the up-current direction (183^∘). The directions of the modeled and measured U_GD are within 5^∘ of each other.

Compared to the relatively high wind condition on 22 November, it is interesting to discuss the results for 24 November (Fig. 13b), with a wind speed of 5.5 m s⁻¹ instead of 11 m s⁻¹. The amplitude of the Doppler velocity is not much reduced, in spite of more than halved current and Stokes drift. This is consistent with the expected near-constant value C₀ of the wave Doppler velocity magnitude, and this is the main result of the present paper.

Figure 14 TSCV retrieval U_CD (in blue) obtained by subtracting the U_WD from the MLM-processed Trèfle buoy data (in midnight blue) from the U_GD vector determined from the KaRADOC measurements (in red) compared to field measurements by HF radar, CARTHE and SVP drifters (shades of gray) on (a) 22 and (b) 24 November.

The process leading to the estimation of the constant C₀ is, however, dependent on a number of assumptions: the directional wave spectrum must be evaluated from the buoy data, then matched to a parametric spectral shape before the necessary numerical integrations can be performed. The (Elfouhaily et al., 1997) spectral shape we have used depends on the wind speed and direction, but also on a wave age parameter, Ω, equal to 0.84 for equilibrium seas.

Table 4 summarizes a subset of the extensive tests we have performed to check the sensitivities of this process. It is clear from this table that drastically changing the wind speed, as occurred between the two days, affects the magnitude of the computed U_WD more at the Ku band than the Ka band, but not in a catastrophic way, and that the wave age parameter Ω can also be varied over its meaningful range quite freely. We have also checked that the transition frequency at which the spectral tail is matched to the observational data is not a very sensitive parameter, provided it is taken low enough for the buoy data to be of good quality where they are kept.

Extracting directional wave spectra from buoy data, however, is a quite an intricate and subjective step. Several methods have been developed over the years to this end, each with pros and cons (see Benoit et al., 1997, for a review). Two of the best-established methods are the maximum entropy method (MEM) and the maximum likelihood method (MLM). The MEM is a parametric method which assumes a specific form of the directional spreading function. In each frequency band, the parameters of the spreading function are chosen such that the first moments of the azimuthal Fourier spectrum match the buoy-derived ones. The MLM is a nonparametric method akin to the Capon beamformer. In terms of directional moments measured by buoys, the MEM estimates provide spectra that exactly fit the measured moments, while the MLM produces spectra that have directional spreads larger than those obtained directly from the measured moments. However, it is not clear how they compare on other properties of the spectrum that may be relevant to the mean slope velocity. Comparing results obtained with these two methods was thus a convenient way to test the sensitivity of U_WD to the sea state directional spread.

Table 4 Modeled wave Doppler velocity amplitude M_WD and direction φ_WD at the Ku and Ka band with directional wave spectra produced using the maximum entropy or maximum likelihood methods and varying the wave age parameter Ω of the Elfouhaily (1997) high-frequency spectral tail. In all cases the transition frequency between the wave data and the high-frequency spectrum is f_t=0.35 Hz. All values are estimated for $\mathit{\theta }={\mathrm{12}}^{\circ }$.

As Table 4 shows, using one technique or the other to estimate the resolved part of the wave directional spectrum does induce significant differences in the simulated U_WD values, showing that the azimuthal width of the spectrum, which is currently not very well constrained observationally, is a sensitive factor. The values obtained using the broader MLM spectrum are consistently smaller than those obtained using the MEM spectrum. A broader azimuthal distribution only redistributes the weight between the different Mss components but reduces the contributions composing the msv₀ vector.

The results obtained using both methods are shown in Fig. 13 as light green and dark green lines. It appears that the MLM processing of the buoy data gives the best fit to the radar U_GD, showing that the directional spread of the sea state should not be taken too low. The possibility that the directional distribution of the Elfouhaily (1997) spectrum could be slightly too narrow for intermediate wavelengths of 2–10 m was, for instance, discussed in specific cases by Peureux et al. (2018). It is, however, not yet clear if it is specific to the very young wind seas they observed, although it could also explain some properties of L-band backscatter (Yueh et al., 2013). The MLM was used to process the Trèfle data in the rest of this study.

Conversely, Fig. 14 illustrates the use of the DV data for the retrieval of the surface current vector by subtraction of U_WD from the fitted U_GD. The norm of the difference between the in situ measured and remotely sensed U_CD vectors is less than 20 cm s⁻¹ on both days, which is significant but quite satisfying at such an early stage of the technique, especially taking into account the fact that geophysical variability due to time variations of wind and tidal current occurred over the several hours of the flight.

Figure 15 Ku-band Doppler measurements performed on (a) 22 November and (b) 24 November with the KuROS radar in the port-looking antenna configuration at a $\mathit{\theta }={\mathrm{12}}^{\circ }$ incidence angle. The graphical conventions are identical for the two plots. The red dots and error bars represent the average and ±1 standard deviation interval of the platform-motion-corrected DV measurements along the different tracks. The magenta dots mark the mean values after correction of the U_AGD contribution. The magenta dotted line is the cosine fit to the corrected data. The blue and midnight blue lines respectively represent the projection along the line-of-sight azimuth of the CARTHE current measurements and the U_WD vector computed from the MLM-processed Trèfle buoy data. The green line represents the sum of these two contributions and should agree with the magenta dotted line.

4.4 Mean Doppler velocity from KuROS

Due to the much broader radiation diagram of the KuROS antenna, analyzing the Ku-band data requires significantly more effort, as the U_AGD spurious velocity contribution due to the azimuthal variation of σ₀ across the FOV discussed in Sect. 2.1 must be compensated for. The DV measurements, corrected for the U_NG platform motion contribution but not for the U_AGD contribution, are represented in Fig. 15 as red dots, while the green line represents the projection along the line-of-sight azimuth of the sum of the TSCV and the MLM-derived U_WD vectors. The difference is clearly very large, reaching 2 m s⁻¹ in places on 22 November, and smaller on 24 November as the azimuthal modulation of the radar NRCS is much weaker.

Figure 16 Comparison of KuROS-derived Doppler velocity, corrected for the U_AGD and U_WD wave contributions, with in situ (CARTHE, SVP drifters and HF radar) current measurements on (a) 22 and (b) 24 November.

Introducing the fits to the Ku-band NRCS data discussed in Sect. 4.2 in Eq. (A24) allows one to produce the corrected data represented by the magenta dots, which are in much better agreement with the green line (though a constant offset is apparent in the 24 November data, which is rejected by the cosine-fit procedure). Figure 16 summarizes the U_CD retrieval operation in vector form (the magenta arrow represents only the first azimuthal harmonic component of the U_AGD correction). The norm of the difference between the in situ (in gray) and remotely sensed (in blue) estimates of the TSCV is of the order of 0.5 m s⁻¹ on 22 November and of the order of 0.2 m s⁻¹ on 24 November. Again, these numbers, though admittedly not small, can be considered encouraging given the number of very large corrections applied to the data and the fact that the instrument had definitely not been designed for this purpose.

4.5 Observed Doppler velocity modulations

The range-resolution scheme implemented in KuROS makes it a very interesting instrument for the analysis of DFS and NRCS modulations. In particular, Caudal et al. (2014), with a different antenna (slightly narrower beam), have attempted to use the cross-spectrum of the DFS and NRCS to resolve the 180^∘ ambiguity in the wave propagation direction. In the SKIM context, analyzing the contribution of the resolved scales to the correlation between σ₀ and the DFS could permit the development of empirical methods to estimate the unresolved part and provide estimates of U_WD.

In practice, with the antenna used for the Drift4SKIM flights, another contribution to the DFS modulations is also caused by the gradients of σ₀ and the speed of the aircraft, just like the mean spurious U_AGD velocity. Brighter areas in the field of view tend to strongly influence the DFS signal towards positive values if they are located to the front of the aircraft and negative values if they are located aft of the aircraft.

As a test of this, simulations were performed with the Radar Sensing Satellite Simulator (Nouguier, 2019), which are illustrated in Fig. 17. The amplitude of the spurious modulations is enhanced by 70 % when the antenna diagram is made 50 % wider in azimuth. With typical variations of σ₀ up to 1 dB over scales of the order of 1 km (e.g., Fig. 9), the variation of σ₀ with azimuth φ is roughly proportional to 1∕sin θ, giving a U_AGD that does not vary much with θ, of the order of 1.5 m s⁻¹. This spurious velocity is larger than the 0.5 m s⁻¹ significant orbital velocity of the swell. As a result the phase relation between DFS and σ₀ can change sign as a function of azimuth due to the combination of two imaging mechanisms with comparable magnitudes and possibly opposite signs.

Figure 17 Qualitative validation of the R3S simulations of the radar imaging mechanism (Nouguier, 2019). Both the real data and simulation contain the geophysical modulation of velocities associated with surface velocities and slopes in the look direction (part of U_GD) as well as aircraft velocities and slopes in the flight direction (part of U_AGD). Note that the wave phases in the R3S simulation are random and cannot be expected to match those in the data or between the two simulations.

This effect will be weaker for shorter (wind-sea) components as soon as the wavelength and crest length become much shorter than the KuROS footprint L_y, as given by Eq. (A36): for a given σ₀ contrast, the gradient increases linearly as the scale L is reduced, but the U_AGD for a given gradient is reduced exponentially in $-L_y/L$.

5 Implications for SKIM

The use of two Doppler radars, in the Ka and Ku band with the same pulse-pair technique but antennas with very different radiation diagrams, has provided important insight for the preparation of the SKIM mission.

Regarding radar measurements, the Drift4SKIM campaign clearly demonstrated the feasibility of the TSCV retrieval approach proposed for SKIM (Ardhuin et al., 2018; ESA, 2019) based on the use of the SKIM wave spectrum measurements (here replaced by in situ buoy measurements) to estimate the wave Doppler velocity contribution U_WD associated with the wave intrinsic phase speed. Measuring the first directional moments (on which the buoy estimates are based) is sufficient to estimate U_WD, and resolving wavelengths of 15 m (a frequency of 0.32 Hz) is sufficient to estimate the full spectral contribution, appending a parametric spectral shape for the unresolved shorter waves. In fact, it is most important to resolve the peak of the wind sea, and a resolved wavelength of 30 m is typically enough for wind speeds higher than 7 m s⁻¹. As this article has shown, however, the angular distribution of the directional spectrum is a sensitive element in both the resolved and parameterized wavelength ranges. Work is still needed to improve the spectral parameterization and to determine whether the accuracy of the sea state restitution algorithms intended for SKIM will be sufficient to solve this issue.

This experiment has also increased confidence in the use of forward models based on the Kirchhoff approximation, such as the R3S of Nouguier (2019), for the study of higher-order effects on the measured DFS. A subject of particular interest is, for instance, the effect of shear in the surface layer on the SKIM DFS, a key to the determination of the effective SKIM measurement depth.

The campaign also stressed the necessity of very good knowledge of the measurement geometry, including the antenna radiation diagram, and the spatial and azimuthal variation of the radar cross section. In this respect, the main characteristics of the instruments used for the present campaign and for the planned SKIM satellite mission are recalled in Table 5, together with the value of the prefactor of the $\sin (\mathit{\phi }-{\mathit{\phi }}_{\mathrm{t}}){\partial }_{\mathit{\phi }}\log \left({\mathit{\sigma }}_{\mathrm{0}}\right)$ term in Eq. (6) of U_AGD (as can be seen in Fig. 10, ∂_φlog (σ₀) is typically 0.1 rad⁻¹ at a 12^∘ incidence angle). As the apparent mispointing due to σ₀ gradients in azimuth or space is proportional to the beamwidth squared, the non-geophysical velocities caused by this effect for SKIM, though non negligible, are actually much smaller than for KuROS, even at a 6^∘ incidence angle.

Table 5 Main differences between the KaRADOC and KuROS airborne radars used in the present article and the SKIM system as presented by the ESA (2019). The factor ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}{V}_{\mathrm{p}}/\mathrm{2}$ is the prefactor of $\sin (\mathit{\phi }-{\mathit{\phi }}_{\mathrm{t}}){\partial }_{\mathit{\phi }}\log \left({\mathit{\sigma }}_{\mathrm{0}}\right)$ in the expression of U_AGD.

As discussed in Sect. 2.3, due to the much reduced platform velocity, the pointing requirements for airborne systems are much easier to reach than for satellite systems for which a pointing accuracy of a few microradians cannot be achieved by attitude measurements alone (gyroscopes and star trackers) but must use a separation of the geophysical and non-geophysical patterns in the data (ESA, 2019). This data-driven approach is also used in airborne systems for correcting phase biases in the antenna diagram (Rodríguez et al., 2018).

Finally, as discussed in Sect. 2.3, we recall that the incidence angle is estimated from the range measurements in the cases of KuROS and SKIM and estimated directly from the platform attitude for the pencil-beam case of KaRADOC. In the spaceborne context, the local slope of the ocean has to be taken into account, as it can induce a mispointing of the nadir beam of up to 300 µrad (Sandwell and Smith, 2014) and induce a correction in the elevation angle at the observation point.

Other radar system constraints or optimizations for satellite systems are discussed by Rodriguez (2018) and the ESA (2019, chap. 5), with sampling issues further analyzed by Chelton et al. (2019).

6 Conclusions and perspectives

The Drift4SKIM campaign clearly demonstrated that surface geophysical velocities can be measured by microwave Doppler radars implementing the pulse-pair method at the Ka band at a 12^∘ incidence angle. The Ku-band measurements, though less easy to interpret due to the large antenna beamwidth of the instrument, also supported this view. The campaign data are consistent with a geophysical model function (GMF) that expresses the geophysical DFS as the sum of the range component of the total surface current velocity and a wave DFS that is a weakly varying function of the sea state of the order of 2.0 m s⁻¹ at the Ka band and 2.4 m s⁻¹ at the Ku band. This wave DFS integrates contributions of all wavenumbers and directions, weighted by the surface slope spectrum. It can be well estimated from the sea surface elevation directional spectrum using the Kirchhoff approximation framework.

The campaign highlighted the importance of very good knowledge of the platform motion and orientation as well as the radar line-of-sight direction vector. The Ku-band NRCS–DFS imagery, though not very successful in that respect, observed a large number of interesting modulation phenomena, which will be analyzed in more detail in forthcoming contributions.

In general, the robustness of the theoretical GMF and its possible empirical adaptation will require the acquisition of more data in a wider range of wind and wave conditions. An in-depth investigation of the angular width of the sea state directional spectrum in the short gravity wave regimes seems of particularly high interest in this respect. Also, obtaining a description of the scale-resolved statistics of sea surface slope skewness would open the path to a Kirchhoff approximation study of the upwind–downwind asymmetry of the radar NRCS and DFS, which is currently lacking.

Finally, the test of near-nadir satellite measurements is limited by the very different viewing geometry due to the difference in altitude. Airborne measurement footprints are at most 500 m or so and thus cannot reproduce the averaging properties of the much wider footprint of a satellite instrument. Still, this medium-sized footprint is comparable to the unfocused SAR resolution that will be obtained with SKIM and provides some practical application with a similar azimuthal averaging that has a limited directional resolution for swell spectrum measurement.

Future airborne systems may ideally combine higher incidence angles, such as that used on DopplerScatt (Rodríguez et al., 2018), OSCAR and Wavemill (Martin et al., 2018), with near-nadir angles that allow for unambiguous wave measurements. In that case, the large azimuthal footprint of KuROS is probably not necessary, and a narrower beam like KaRADOC can be used, greatly simplifying the analysis.

Appendix A: Doppler scatterometry theory

This Appendix proposes an extension of the theory of pencil-beam Doppler scatterometry exposed in Rodriguez (2018) and Rodríguez et al. (2018) to the case of near-nadir fan-beam instruments such as SKaR and KuROS. It compiles a number of processing steps or concepts that had to be developed for the analysis of the Drift4SKIM KuROS data. In each section the differences from and similarities to the spaceborne SKIM context are highlighted.

A1 Pulse-pair theory

A1.1 Radar pulse-pair measurements

A radar instrument works by sending microwave pulses into the environment and recording the echo from its field of view. Usual scatterometers consider only the intensity of the return signal. Coherent instruments, such as SARs, measure both the amplitude of the return signal and its phase with respect to the transmitted carrier as a function of range. Over the ocean, the phase of the return signal for a single pulse is random and uniformly distributed over the unit circle. The radar returns of successive pulses transmitted at short intervals are, however, correlated, and the time history of the phase can be used to measure the relative motion of the radar and the scatterers. SARs make use of this property to refine the along-track resolution of backscattering cross-sectional measurements. SKIM and the other proposed Doppler missions aim to use it to obtain direct surface current measurements.

As explained by Rodriguez (2018, Appendix A), the complex amplitude of the return signal of a pulse transmitted at time t_i can be expressed as (here and in the following, the × symbol between scalar operands denotes the usual scalar multiplication; there is no occurrence of the vector product in this article)

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{E}_i(t_i,r^{\prime })=& n(t_i,r^{\prime })+{\displaystyle \frac{A\left(r^{\prime }\right)}{{r^{\prime }}^{\mathrm{2}}}}\\ \times & \int G(t_i,\mathit{x})\mathit{\chi }(r^{\prime }-r(t_i,\mathit{x}\left)\right)\\ & \times \exp \left[-\mathrm{2}ikr(t_i,\mathit{x})\right]s(t_i,\mathit{x})\text{d}S,\end{array}\end{array}$$

where the integral is performed over the sea surface, A(r^′) is a time-independent weakly dependent function of range, unimportant for our purposes here (corresponding in particular to the effects of transmitted signal amplitude, receiver and processing gain, and attenuation losses), G(x) is the one-way antenna diagram, χ(r) is the range-point-target response of the instrument, r^′ is the nominal pixel range in the time sampled signal, $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the radar wavenumber, r(t_i,x) is the range from the radar to the observation point x at time t_i, $n(t_i,r^{\prime })$ is the thermal noise contribution, and s(t_i,x) is the complex reflection coefficient of the sea surface at instant t_i and location x.

As mentioned by Rodríguez et al. (2018), the thermal noise contribution, though it plays a major role in the quality of the measurements, is conceptually simple and can be safely considered δ-correlated in time and characterized by a single quantity, its average power N. The reflection coefficient s(t_i,x), on the other hand, emerges from the interaction of the electromagnetic waves with the ocean surface and has much richer physics. It is affected by electromagnetic phenomena as well as the geometry and kinematics of the sea surface itself, and its statistics are further complicated by the so-called “speckle” phenomenon. As stated by Rodríguez et al. (2018), the correlation function of this coefficient as a function of time and space separation, averaged over speckle realizations, can be modeled as

$$\begin{array}{}\text{(A2)}& {〈s(t,\mathit{x})s^{\ast }\left(t^{\prime },{\mathit{x}}^{\mathbf{\prime }}\right)〉}_S=\mathit{\delta }\left(\mathit{x}-{\mathit{x}}^{\mathbf{\prime }}\right){\mathit{\sigma }}_{\mathrm{0}}(t,\mathit{x}){\mathit{\gamma }}_{\mathrm{TS}}\left(|t-t^{\prime }|\right),\end{array}$$

with σ₀(t,x) the normalized radar backscattering cross section (NRCS) in the appropriate polarization and ${\mathit{\gamma }}_{\mathrm{TS}}\left(\right|\mathit{\tau }\left|\right)$ a function describing its time decorrelation at a fixed location due to the life history of individual scattering patches.

The so-called pulse-pair technique of Zrnic (1977) relies on the properties of the product of the return signals from consecutive radar pulses. Combining Eqs. (A1) and (A2) to compute the speckle-averaged product of the return signals for two radar pulses sent at t₁ and $t_{\mathrm{2}}=t_{\mathrm{1}}+\mathrm{\Delta }t$, with Δt the pulse repetition interval (PRI), one obtains

$$\begin{array}{}\text{(A3)}& {\mathrm{PP}}_{\mathrm{\Delta }t}\left(t_{\mathrm{1}},r^{\prime }\right)={〈E_{\mathrm{2}}\left(t_{\mathrm{2}}=t_{\mathrm{1}}+\mathrm{\Delta }t,r^{\prime }\right)E_{\mathrm{1}}{\left(t_{\mathrm{1}},r^{\prime }\right)}^{\ast }〉}_S\end{array}$$

as

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\text{PP}}_{\mathrm{\Delta }t}(t_{\mathrm{1}},r^{\prime })=& {\displaystyle \frac{A^{\mathrm{2}}\left(r^{\prime }\right)}{{r^{\prime }}^{\mathrm{4}}}}{\mathit{\gamma }}_{\text{TS}}\left(\right|\mathrm{\Delta }t\left|\right)\\ \times & \int G^{\mathrm{2}}(t_{\mathrm{1}},\mathit{x}){\mathit{\chi }}^{\mathrm{2}}(r^{\prime }-r(t_{\mathrm{1}},\mathit{x}\left)\right){\mathit{\sigma }}_{\mathrm{0}}(t_{\mathrm{1}},\mathit{x})\\ & \times \exp \left[-\mathrm{2}ik\left[r(t_{\mathrm{1}}+\mathrm{\Delta }t,\mathit{x})-r(t_{\mathrm{1}},\mathit{x})\right]\right]\text{d}S.\end{array}\end{array}$$

As can be seen in this equation, the phase of the pulse-pair signal contains a weighted average of the time rate of change of the distance separating the radar from the scattering elements in its instantaneous footprint. This rate of change can be interpreted as a velocity.

A1.2 Measurement geometry

Figure 1a and b summarize the acquisition geometry in the airborne and spaceborne settings. The antenna radiation diagram G²(t₁,x) is represented as gray shading of the sea surface, while the range-point response function ${\mathit{\chi }}^{\mathrm{2}}(t_{\mathrm{1}},r^{\prime }-r(t_{\mathrm{1}},\mathit{x}\left)\right)$ is represented as white grating. In Eq. (A4), we have made the assumptions that $G(t_{\mathrm{1}},\mathit{x})=G(t_{\mathrm{2}},\mathit{x})$ and $\mathit{\chi }(r^{\prime }-r(t_{\mathrm{1}},\mathit{x}\left)\right)=\mathit{\chi }(r^{\prime }-r(t_{\mathrm{2}},\mathit{x}\left)\right)$, neglecting the effect of the spatial translation of the beam illumination pattern and range-resolution weighting distribution on the sea surface.

This is a very good approximation for airborne pulse-pair radar observations and a quite good one for spaceborne observations. For airborne instruments, the PRI is usually chosen such that the line-of-sight projection of the platform movement over a PRI is smaller than one-half the carrier wavelength to avoid phase ambiguity. For spaceborne instruments, avoiding ambiguity is not practical due to the much larger platform velocity, but the PRI is constrained by other considerations, and the platform displacement over a PRI is much smaller than the characteristic scales of the antenna radiation diagram and of the range-point response.

A1.3 Pulse-pair signal approximation

Returning to Eq. (A4), we see that over the time interval separating the two radar pulses, the radar has moved from its original position x_R(t₁) to x_R(t₁)+V_PΔt, and the scatterer located at x has moved to x+v_sΔt. The radar-to-scatterer vector has thus changed by [v_s(x)−V_P]Δt. The distance change can be approximated by

$$\begin{array}{}\text{(A5)}& r(t_{\mathrm{1}}+\mathrm{\Delta }t,\mathit{x})-r(t_{\mathrm{1}},\mathit{x})=\mathrm{\Delta }t{\displaystyle \frac{\mathit{x}-{\mathit{x}}_{\mathrm{R}}\left(t_{\mathrm{1}}\right)}{\left|\right|\mathit{x}-{\mathit{x}}_{\mathrm{R}}\left(t_{\mathrm{1}}\right)\left|\right|}}\cdot \left({\mathit{v}}_{\mathrm{s}}-{\mathit{V}}_{\mathrm{R}}\right),\end{array}$$

where the neglected terms are of the order of $\mathrm{\Delta }{t}^{\mathrm{2}}\left|\right|{\mathit{v}}_{\mathrm{s}}-{\mathit{v}}_{\mathrm{R}}|{|}^{\mathrm{2}}/||\mathit{x}-{\mathit{x}}_{\mathrm{R}}\left(t_{\mathrm{1}}\right)|{|}^{\mathrm{2}}$. Introducing

$$\begin{array}{}\text{(A6)}& \mathit{e}\left(\mathit{x}\right)={\displaystyle \frac{\mathit{x}-{\mathit{x}}_{\mathrm{R}}\left(t_{\mathrm{1}}\right)}{|\mathit{x}-{\mathit{x}}_{\mathrm{R}}\left(t_{\mathrm{1}}\right)|}},\end{array}$$

the unit vector pointing from the radar location at t₁ to the observation point (choosing either time instant is equivalent, as the difference is of the same order of magnitude as the neglected terms), the pulse-pair signal can be expressed as

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\text{PP}}_{\mathrm{\Delta }t}(t_{\mathrm{1}},r^{\prime })=& {\displaystyle \frac{A^{\mathrm{2}}\left(r^{\prime }\right)}{{r^{\prime }}^{\mathrm{4}}}}{\mathit{\gamma }}_{\text{TS}}\left(\mathrm{\Delta }t\right)\\ \times & \int G^{\mathrm{2}}(t_{\mathrm{1}},\mathit{x}){\mathit{\chi }}^{\mathrm{2}}(r^{\prime }-r(t_{\mathrm{1}},\mathit{x}\left)\right){\mathit{\sigma }}_{\mathrm{0}}(t_{\mathrm{1}},\mathit{x})\\ & \times \exp \left[\mathrm{2}ik\mathrm{\Delta }t\mathit{e}\left(\mathit{x}\right)\cdot \left({\mathit{V}}_{\text{R}}-{\mathit{v}}_{\text{s}}\left(\mathit{x}\right)\right)\right]\text{d}S.\end{array}\end{array}$$

This equation is not very practical, as the relative motion of the scatterers with respect to the radar enters as the argument of an exponential integrand. Obtaining an equivalent representation as the exponential of a sum of weighted integrals would be desirable. Introducing the effective illuminated surface,

$$\begin{array}{}\text{(A8)}& S\left(t_{\mathrm{1}},r^{\prime }\right)=\int G^{\mathrm{2}}\left(t_{\mathrm{1}},\mathit{x}\right){\mathit{\chi }}^{\mathrm{2}}\left(r^{\prime }-r(t_{\mathrm{1}},\mathit{x})\right)\mathrm{d}S,\end{array}$$

the normalized weighting function,

$$\begin{array}{}\text{(A9)}& W\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right)={\displaystyle \frac{G^{\mathrm{2}}\left(t_{\mathrm{1}},\mathit{x}\right){\mathit{\chi }}^{\mathrm{2}}\left(r^{\prime }-r(t_{\mathrm{1}},\mathit{x})\right)}{S(t_{\mathrm{1}},r^{\prime })}},\end{array}$$

the average and fluctuating parts of the NRCS

$$\begin{array}{}\text{(A10)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(t_{\mathrm{1}},r^{\prime }\right)=\int W\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right){\mathit{\sigma }}_{\mathrm{0}}\left(t_{\mathrm{1}},\mathit{x}\right)\mathrm{d}S,\text{(A11)}& {\displaystyle }\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{0}}\left(t_{\mathrm{1}},\mathit{x}\right)}{\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(t_{\mathrm{1}},r^{\prime }\right)}},\end{array}$$

and borrowing the algebraic technique of “cumulant expansion” from probability theory, it is possible to express PP_Δt as

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\text{PP}}_{\mathrm{\Delta }t}(t_{\mathrm{1}},r^{\prime })=& {\displaystyle \frac{A^{\mathrm{2}}\left(r^{\prime }\right)}{{r^{\prime }}^{\mathrm{4}}}}{\mathit{\gamma }}_{\text{TS}}\left(\mathrm{\Delta }t\right)\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime })S(t_{\mathrm{1}},r^{\prime })\\ & \times \exp \left[\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{(\mathrm{i}\mathrm{2}k\mathrm{\Delta }t{)}^n}{n\mathrm{!}}}{\mathit{\kappa }}_n\right],\end{array}\end{array}$$

with κ_n the successive “cumulants” of $\mathit{e}\left(\mathit{x}\right)\cdot ({\mathit{V}}_{\mathrm{R}}-{\mathit{v}}_{\mathrm{s}}(\mathit{x}\left)\right)$ with respect to the “density distribution” $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})W(t_{\mathrm{1}},r^{\prime },\mathit{x})$. As all the κ_n are real, we see that odd-n terms contribute to the argument of the pulse-pair signal, while even-n terms contribute to its magnitude. Keeping only the first two terms in the sum, one obtains

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\text{PP}}_{\mathrm{\Delta }t}(t_{\mathrm{1}},r^{\prime })=& {\displaystyle \frac{A^{\mathrm{2}}\left(r^{\prime }\right)}{{r^{\prime }}^{\mathrm{4}}}}{\mathit{\gamma }}_{\text{TS}}\left(\mathrm{\Delta }t\right)\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime })S(t_{\mathrm{1}},r^{\prime })\\ & \times \exp \left[\mathrm{i}\mathrm{2}k\mathrm{\Delta }t{\mathit{\kappa }}_{\mathrm{1}}\right]\exp \left[-\mathrm{2}(k\mathrm{\Delta }t{)}^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{2}}\right].\end{array}\end{array}$$

As expected, the expression of κ₁,

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{\mathit{\kappa }}_{\mathrm{1}}(t_{\mathrm{1}},r^{\prime })=\int & W(t_{\mathrm{1}},r^{\prime },\mathit{x})\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})\\ & \times \mathit{e}\left(\mathit{x}\right)\cdot \left({\mathit{V}}_{\text{R}}-{\mathit{v}}_{\text{s}}\left(\mathit{x}\right)\right)\text{d}S,\end{array}\end{array}$$

shows that to first order the argument of the pulse-pair signal gives access to the integral over the footprint of the relative velocity of the scatterers with respect to the radar. The expression of κ₂,

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}{\mathit{\kappa }}_{\mathrm{2}}(t_{\mathrm{1}},r^{\prime })=\int & W(t_{\mathrm{1}},r^{\prime },\mathit{x})\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})\\ & \times {\left[\mathit{e}\left(\mathit{x}\right)\cdot ({\mathit{V}}_{\text{R}}-{\mathit{v}}_{\text{s}}(\mathit{x}\left)\right)-{\mathit{\kappa }}_{\mathrm{1}}\right]}^{\mathrm{2}}\text{d}S,\end{array}\end{array}$$

is a description of the impact of the variability of e(x), $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ and v_s inside the footprint on the pulse-pair signal magnitude.

A1.4 Pulse-pair signal phase approximation

Working now in the Earth-fixed reference frame at the observation point, we define

$$\begin{array}{}\text{(A16)}& V_{\mathrm{GD}}=-\int W\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right)\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right)\mathit{e}\left(\mathit{x}\right)\cdot {\mathit{v}}_{\mathrm{s}}\left(\mathit{x}\right)\mathrm{d}S,\end{array}$$

the (geophysically relevant) weighted projection of the scatterer velocity in that frame on the radar line of sight, and

$$\begin{array}{}\text{(A17)}& V_{\mathrm{NG}}\left(t_{\mathrm{1}},r^{\prime }\right)={\mathit{V}}_{\mathrm{R}}\cdot \int W(t_{\mathrm{1}},r^{\prime },\mathit{x})\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(t_{\mathrm{1}},r^{\prime },\mathit{x}\right)\mathit{e}\left(\mathit{x}\right)\mathrm{d}S,\end{array}$$

the (non-geophysical) projection of the radar velocity (our conventions are such that V_GD is positive when the scatterers move towards the radar and that V_NG is positive when the radar moves towards the footprint, in keeping with everyday intuition).

With these conventions, one sees that

$$\begin{array}{}\text{(A18)}& V_{\mathrm{GD}}\left(t_{\mathrm{1}},r^{\prime }\right)={\mathit{\kappa }}_{\mathrm{1}}\left(t_{\mathrm{1}},r^{\prime }\right)-V_{\mathrm{NG}}\left(t_{\mathrm{1}},r^{\prime }\right).\end{array}$$

Using Eq. (A13), one can obtain κ₁ approximately as 1∕(2kΔt) times the argument of the complex pulse-pair signal. At this stage, one must, however, consider a bit carefully the ambiguity that is inherent in phase measurements. As the phase of a complex number is only known up to a multiple of 2π, κ₁ is only obtained up to a multiple of $\frac{\mathit{\lambda }}{\mathrm{2}\mathrm{\Delta }t}$. This effect can be neglected as long as both V_NG and κ₁ remain within the unambiguous interval $\left[-\frac{\mathit{\lambda }}{\mathrm{4}\mathrm{\Delta }t};\frac{\mathit{\lambda }}{\mathrm{4}\mathrm{\Delta }t}\right]$. For larger platform velocities, care must be taken to add the right multiple of $\frac{\mathit{\lambda }}{\mathrm{2}\mathrm{\Delta }t}$ to κ₁ before subtracting V_NG. For airborne instruments it is usually feasible to select a small enough PRI to avoid ambiguity altogether. For satellite instruments, one approach is to select a solid Earth-fixed reference frame, in which v_s is small, and to work on the phase-migrated pulse-pair signal

$$\begin{array}{}\text{(A19)}& {\stackrel{\mathrm{̃}}{\mathrm{PP}}}_{\mathrm{\Delta }t}\left(t_{\mathrm{1}},r^{\prime }\right)=\exp \left[-\mathrm{i}\mathrm{2}k\mathrm{\Delta }t{V}_{\mathrm{NG}}\right]\times {\mathrm{PP}}_{\mathrm{\Delta }t}\left(t_{\mathrm{1}},r^{\prime }\right).\end{array}$$

It is easy to see that V_GD can be retrieved as

$$\begin{array}{}\text{(A20)}& V_{\mathrm{GD}}\left(t_{\mathrm{1}},r^{\prime }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}k\mathrm{\Delta }t}}\arg \left({\stackrel{\mathrm{̃}}{\mathrm{PP}}}_{\mathrm{\Delta }t}\left(t_{\mathrm{1}},r^{\prime }\right)\right).\end{array}$$

At this stage, even a coarse approximation of V_NG can be used, as long as it is sufficient to resolve the phase ambiguity. This is important in particular for the onboard processors of satellite instruments, which have to rely on a limited quality of information in terms of position, velocity and beam pointing and typically cannot use the $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ distribution information that ground segment processors can retrieve from the signal. The correction applied by the onboard processor must, however, be accounted for in later processing stages.

A2 Non-geophysical contribution V_NG

The non-geophysical contribution V_NG must be estimated from the platform velocity and radar beam pointing. For pulse-limited instruments such as KuROS and SKaR, the incidence angle is determined in each range bin as a function of the altitude. The accuracies of the range-resolution and altitude determination processes are then critical. Last, asymmetric azimuthal variations of the sea surface NRCS within a given range bin tend to bias the effective observation azimuth towards the brighter part of the instrument FOV. This section discusses these different aspects.

A2.1 Beam pointing accuracy

We work in the simplified setting of the flat-Earth approximation, in which the elevation and incidence angles γ and θ are equal. We use a platform-fixed reference frame, the origin of which is located at the antenna phase center of the instrument, with the x vector pointing to the geometric front of the platform, the y vector pointing to starboard, the z vector pointing to the floor and a local geographic NED reference frame, the origin of which is fixed to the solid Earth and located at a suitable point of the campaign area.

The orientation of the platform-fixed reference frame with respect to the local geographic frame is provided by the platform IMU as (roll, pitch, heading) Euler angles, from which one can construct the direction cosine matrix

$$\begin{array}{}\text{(A21)}& \text{DCM}=\left[\begin{array}{ccc}{c}_{\mathrm{p}}{c}_{\mathrm{h}}& s_{\mathrm{r}}{s}_{\mathrm{p}}{c}_{\mathrm{h}}-c_{\mathrm{r}}{s}_{\mathrm{h}}& c_{\mathrm{r}}{s}_{\mathrm{p}}{c}_{\mathrm{h}}+s_{\mathrm{r}}{s}_{\mathrm{h}}\\ c_{\mathrm{p}}{s}_{\mathrm{h}}& s_{\mathrm{r}}{s}_{\mathrm{p}}{s}_{\mathrm{h}}+c_{\mathrm{r}}{c}_{\mathrm{h}}& c_{\mathrm{r}}{s}_{\mathrm{p}}{s}_{\mathrm{h}}-s_{\mathrm{r}}{c}_{\mathrm{h}}\\ -s_{\mathrm{p}}& s_{\mathrm{r}}{c}_{\mathrm{p}}& c_{\mathrm{r}}{c}_{\mathrm{p}}\end{array}\right],\end{array}$$

allowing one to express the components of a vector in the NED frame from its ($x,y,z$) components in the platform-fixed frame. The two reference frames are consistent in the sense that the frame vectors coincide when the platform is in constant-altitude flight towards the north. In the above expression we have used the transparent notation c_p→cos (pitch), c_r→cos (roll), c_h→cos (heading),
s_p→sin (pitch), s_r→sin (roll) and s_h→sin (heading). Other quantities worth introducing are the course c and glide angle g such that the plane velocity vector in the NED frame is

$$\begin{array}{}\text{(A22)}& {\mathit{V}}_{\mathrm{R}}=V_{\mathrm{R}}\left[\cos \left(g\right)\cos \left(c\right)\mathit{N}+\cos \left(g\right)\sin \left(c\right)\mathit{E}+\sin \left(g\right)\mathit{D}\right].\end{array}$$

In the NED frame, the pointing vector e can be expressed as

$$\begin{array}{}\text{(A23)}& \mathit{e}=\sin \left(\mathit{\theta }\right)\left[\cos \left(\mathit{\phi }\right)\mathit{N}+\sin \left(\mathit{\phi }\right)\mathit{E}\right]+\cos \left(\mathit{\theta }\right)\mathit{D}.\end{array}$$

Its components in the platform-fixed frame can be determined using the fact that ${\mathrm{DCM}}^{-\mathrm{1}}={\mathrm{DCM}}^T$. The corresponding antenna azimuth and elevation angles φ and γ, in terms of which the radiation diagram is specified, can then be expressed using the platform-fixed to antenna-fixed reference frame transformation matrix.

With these notations and using Eq. (A17), one can express V_NG as

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}{V}_{\text{NG}}& (t_{\mathrm{1}},r^{\prime })=V_{\text{R}}\int W(t_{\mathrm{1}},r^{\prime },\mathit{x})\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})\\ & \times \left[\cos \left(g\right)\sin \left(\mathit{\theta }\right)\cos (\mathit{\phi }-c)+\sin \left(g\right)\cos \left(\mathit{\theta }\right)\right]\text{d}S.\end{array}\end{array}$$

Constant-altitude flight corresponds to g≃0. We thus concentrate on the impact of errors in the first term on the right-hand side of this equation. Quite clearly, the impact of errors in sin (θ) is largest when the instrument views the area where cos (φ−c) is large, i.e., in the up-track and down-track directions, while the impact of errors in the azimuthal direction is largest when the instrument looks cross-track (i.e., where the derivative of cos (φ−c) is close to 1).

Leaving aside for the moment the effects of uncertainties on $W(t_{\mathrm{1}},r^{\prime },\mathit{x})$ and $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})$, one sees that at a 12^∘ incidence angle and for a platform velocity of 7000 m s⁻¹ (spaceborne instrument), the SKIM 40 cm s⁻¹ error budget for horizontal velocity measurements translates to pointing accuracies of 4.5 and 21 µrad in incidence angle and azimuth, respectively (see the discussion in Sect. 2.3). In the airborne case at 120 m s⁻¹ platform velocity and 3000 m of altitude, the corresponding numbers are 0.26 and 1.25 mrad for incidence angle and azimuth pointing accuracy for KuROS, respectively. In the cross-track viewing geometry of KaRADOC, only the comparatively mild (but still quite demanding) 1.25 mrad azimuth pointing accuracy requirement applies.

Figure A1 (a) True pointing. The attitude drift changes the antenna footprint direction and shifts the DFS centroid. Here γ and φ are the expected coordinates of the antenna gain ground projection, while $\stackrel{\mathrm{̃}}{\mathit{\gamma }}$ and $\stackrel{\mathrm{̃}}{\mathit{\phi }}$ are the shifted versions of these coordinates by the attitude misknowledge (adapted from Ardhuin et al., 2019b). (b) Apparent pointing φ_a for the SKIM geometry. Examples of two-way antenna gain G as a function of azimuth and distortions (exaggerated 100 times) induced by σ₀ gradients on the power integrated by the radar in the azimuth direction across the antenna diagram (gray curve). Three examples of asymmetric distortions are given: a sine function with $\mathit{\nu }=\sin \mathit{\theta }/{\mathit{\sigma }}_{\mathit{\alpha }}$, a sine function varying 3 times faster and a linear trend. Such distortions induce an apparent mispointing of the beam δφ and a correction to the geometrical line-of-sight relative velocity estimate.

Figure A1a shows the measurement geometry seen from above. One can see that uncertainties on the viewing azimuth and incidence angle have different origins.

• The uncertainty in azimuth can be due to imperfect knowledge of the weighting corresponding to the $W(t_{\mathrm{1}},r^{\prime },\mathit{x})\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{x})$ term in Eq. (A24). This can of course come from imperfect platform attitude or antenna orientation information but also from an imperfect characterization of the antenna radiation diagram or the distribution of σ₀ on the sea surface.

• The uncertainty in incidence angle is due to imperfect knowledge of the radial position of the range-resolution bins (yellow striping of the footprint in Fig. A1a). This can be due to imperfect timing accuracy or to imperfect knowledge of the vertical separation between the instrument and sea surface.

A2.2 Timing and altitude accuracy

For this brief discussion of the effects of timing and altitude accuracy on incidence angle estimation, we consider a single range bin whose “true” range from the radar is r, whose altitude with respect to the radar is H and the incidence angle is θ. In this case $\mathit{\theta }=\arccos (H/r)$. If the radar now suffers from a timing error δr, the instrument will detect a false altitude H−δr but will ascribe to range bin r−δr the signal coming from r. In the meantime, we consider the surface-tracking algorithm to suffer from an error δh and detect the surface at range $H-\mathit{\delta }r-\mathit{\delta }h$. The data from this range bin will thus be processed using an angle of incidence

$$\begin{array}{}\text{(A25)}& \mathit{\theta }+\mathit{\delta }\mathit{\theta }=\arccos \left({\displaystyle \frac{H-\mathit{\delta }r-\mathit{\delta }h}{r-\mathit{\delta }r}}\right),\end{array}$$

different from the correct value by

$$\begin{array}{}\text{(A26)}& \mathit{\delta }\mathit{\theta }\simeq {\displaystyle \frac{\mathrm{1}}{H\tan \left(\mathit{\theta }\right)}}\left[\mathit{\delta }h+\mathit{\delta }r\left[\mathrm{1}-\cos \left(\mathit{\theta }\right)\right]\right].\end{array}$$

Considering δh and δr to be independent, we see that at 12^∘ the incidence angle knowledge requirements expressed above for SKIM and KuROS respectively translate to timing accuracy requirements of 36.8 and 7.7 m and to surface-tracking accuracy requirements of 80.4 and 16.9 cm.

The timing accuracy requirements are easily met in the spaceborne context but can be challenging in the cost-constrained context of an airborne instrument.

The surface-tracking algorithm, however, does not benefit from the error compensation that exists for the timing error. The requirement for SKIM is easily met by the nadir altimeter payload of SKIM. The 16.9 cm altitude tracking requirement is out of reach of the KuROS airborne instrument. Our analysis of its DFS data will thus be restricted to the side-looking configurations for which, as per Eq. (A24), the pointing requirements are much milder.

A2.3 Effective pointing and azimuth gradient DFS

As expressed in Eq. (A24), for each range-resolution cell V_NG results from an integral over azimuth with a weight that depends on the product of the antenna radiation diagram and the sea surface NRCS, which varies as a function of the horizontal position (x,y) due to the presence of waves, varying winds, currents, surfactants, sea ice and all the physical properties of the sea surface.

Even with perfect knowledge of the platform attitude and velocity, NRCS variations can thus make the effective pointing of the measurements deviate from the pure geometric estimates. Valuable insight into this effect can be gained by considering the saddle-point approximation of Eq. (A24) in the limit of a very narrow antenna diagram (which is clearly applicable for SKIM and KaRADOC, less so for KuROS).

Considering first the case of an antenna pointing towards azimuth φ_b with an infinitely narrow radiation diagram, we see that the product $W(t_{\mathrm{1}},r^{\prime },\mathit{\phi })\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}(t_{\mathrm{1}},r^{\prime },\mathit{\phi })$ is well approximated by the Dirac distribution δ(φ−φ_b). In this limit,

$$\begin{array}{}\text{(A27)}& \begin{array}{rl}& V_{\text{NG}}(t_{\mathrm{1}},r^{\prime })=\\ & V_{\text{R}}\left[\cos \left(g\right)\sin \left(\mathit{\theta }\right)\cos ({\mathit{\phi }}_{\text{b}}-c)+\sin \left(g\right)\cos \left(\mathit{\theta }\right)\right].\end{array}\end{array}$$

We recognize in this expression V_geo, the estimate of V_NG one would have derived using direct geometric arguments.

The essence of the argument is that the sharpest factor in the integral is the beam radiation diagram. If it is now not infinitely sharp, we see that the effect of a gradient of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ is to shift the peak of the distribution by an angle

$$\begin{array}{}\text{(A28)}& \mathit{\delta }\mathit{\phi }=-{\displaystyle \frac{{\left.{\partial }_{\mathit{\phi }}\log \left(\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\right)\right|}_{{\mathit{\phi }}_{\mathrm{b}}}}{{\left.{\partial }_{\mathit{\phi }\mathit{\phi }}\log \left(W\right)\right|}_{{\mathit{\phi }}_{\mathrm{b}}}}}.\end{array}$$

Assuming for $W\left(t_{\mathrm{1}},t^{\prime },\mathit{\phi }\right)$ a Gaussian approximation,

$$\begin{array}{}\text{(A29)}& W(t_{\mathrm{1}},r^{\prime },\mathit{\phi })={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathit{\pi }}{\mathit{\sigma }}_{\mathit{\phi }}\left(r^{\prime }\right)}}\exp \left[-{\displaystyle \frac{(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{b}}{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\left(r^{\prime }\right)}}\right],\end{array}$$

in which σ_φ(r^′) is a parameter describing the width of the antenna diagram at the working incidence angle, one obtains

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}& V_{\mathrm{NG}}(t_{\mathrm{1}},r^{\prime })=V_{\text{R}}\\ & \times \left[\cos \left(g\right)\sin \left(\mathit{\theta }\right)\cos ({\mathit{\phi }}_{\text{b}}-c+\mathit{\delta }\mathit{\phi })+\sin \left(g\right)\cos \left(\mathit{\theta }\right)\right].\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A31)}& \mathit{\delta }\mathit{\phi }={\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\left(r^{\prime }\right)}{\mathrm{2}}}{\partial }_{\mathit{\phi }}\log \left(\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\right).\end{array}$$

Alternatively, one can choose to express V_NG as the sum of V_geo, the geometric approximation, plus an azimuth gradient Doppler velocity contribution,

$$\begin{array}{}\text{(A32)}& V_{\mathrm{NG}}\left(t_{\mathrm{1}},r^{\prime }\right)=V_{\mathrm{geo}}(t_{\mathrm{1}},r^{\prime })+V_{\mathrm{AGD}}\left(t_{\mathrm{1}},r^{\prime }\right),\end{array}$$

with

$$\begin{array}{}\text{(A33)}& V_{\mathrm{geo}}(t_{\mathrm{1}},r^{\prime })=V_{\text{R}}\left[\cos \left(g\right)\sin \left(\mathit{\theta }\right)\cos ({\mathit{\phi }}_{\text{b}}-c)+\sin \left(g\right)\cos \left(\mathit{\theta }\right)\right]\end{array}$$

and

$$\begin{array}{}\text{(A34)}& \begin{array}{rl}{V}_{\mathrm{AGD}}(t_{\mathrm{1}},r^{\prime })=& -V_{\text{R}}\cos \left(g\right)\sin \left(\mathit{\theta }\right)\sin ({\mathit{\phi }}_{\text{b}}-c)\\ & {\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\left(r^{\prime }\right)}{\mathrm{2}}}{\partial }_{\mathit{\phi }}\log \left(\stackrel{\mathrm{̃}}{{\mathit{\sigma }}^{\mathrm{0}}}\right).\end{array}\end{array}$$

One can see from these expressions that for a given azimuthal variation of the NRCS the order of magnitude of V_AGD is set by the width of the antenna radiation diagram: instruments with a thin diagram, such as SKIM and KaRADOC, are less affected than instruments with a broader diagram, such as KuROS. Also, one sees that V_AGD is largest when the instrument looks in the cross-track direction and is zero in the up-track and down-track viewing directions. Finally, one sees that V_AGD is equivalent to the line-of-sight projection of a spurious horizontal velocity U_AGD, which varies with incidence angle only through the variations of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$ and σ_φ:

$$\begin{array}{}\text{(A35)}& U_{\mathrm{AGD}}\left(t_{\mathrm{1}},r^{\prime }\right)=-V_{\mathrm{R}}\sin \left({\mathit{\phi }}_{\mathrm{b}}-c\right){\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\left(r^{\prime }\right)}{\mathrm{2}}}{\partial }_{\mathit{\phi }}\log \left(\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\right).\end{array}$$

At small scales, spatial gradients add to the azimuthal gradient and also induce a spurious velocity with the same expression as a function of $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}$. Using the simple case of a single Fourier component $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}=\mathit{\epsilon }\sin \left[\mathit{\nu }\left(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{b}}\right)\right]$ allows one to evaluate the importance of different scales. The azimuthal shift can be obtained as

$$\begin{array}{}\text{(A36)}& \mathit{\delta }\mathit{\phi }=\mathit{\epsilon }\exp \left(-{\displaystyle \frac{\left({\mathit{\nu }}^{\mathrm{2}}+\mathrm{1}\right){\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}}{\mathrm{4}}}\right)\sinh \left({\displaystyle \frac{\mathit{\nu }{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}}{\mathrm{2}}}\right).\end{array}$$

In the slow variation limit ($\mathit{\nu },{\mathit{\sigma }}_{\mathit{\phi }}\to \mathrm{0}$) and Eq. (A36) this expression coincides with Eq. (A31). For faster variations, one sees that the largest disturbance is obtained when $\mathit{\nu }\sim \sqrt{\mathrm{2}}/{\mathit{\sigma }}_{\mathit{\phi }}$. This azimuthal wavenumber is such that the footprint can host a bright and a dark patch, one on either side of the look direction. This configuration creates the largest disturbance for a given value of the brightness contrast ε. δφ in this case is given by

$$\begin{array}{}\text{(A37)}& \mathit{\delta }{\mathit{\phi }}_{max}=\mathit{\epsilon }{\mathit{\sigma }}_{\mathit{\phi }}{\mathrm{e}}^{-\mathrm{1}/\mathrm{2}}/\sqrt{\mathrm{2}}.\end{array}$$

Appendix B: KuROS antenna diagram determination

A precise determination of the antenna diagram is necessary for any Doppler application, given the possibly large contribution of pointing error in the estimation of the non-geophysical DFS and the effect of the antenna beamwidth in the spurious azimuth gradient velocity U_AGD. A comprehensive strategy has thus been developed for estimating the one-way antenna diagram in amplitude and phase by combining anechoic chamber measurements and verification using the campaign data with a final adjustment of systematic phase shifts in the data. In this section α and β are respectively the latitude and longitude of a set of spherical coordinates centered on the antenna such that the main lobe extends in a longitudinal sector on the Equator $\mathit{\alpha }={\mathrm{0}}^{\circ }$, and the rotation axis of the antenna turntable points towards $\mathit{\alpha }={\mathrm{0}}^{\circ }$ and $\mathit{\beta }={\mathrm{0}}^{\circ }$. With this choice of coordinates the antenna diagram has separable Gaussian dependencies on α and β. In constant-altitude flight, when the antenna points towards φ_b, $\sin \left(\mathit{\alpha }\right)=\sin \left(\mathit{\theta }\right)\sin (\mathit{\phi }-{\mathit{\phi }}_{\mathrm{b}})$ and $\tan \left(\mathit{\beta }\right)=\tan \left(\mathit{\theta }\right)\cos (\mathit{\phi }-{\mathit{\phi }}_{\mathrm{b}})$.

B1 Fixed-antenna NRCS correction

The anechoic chamber measurements are very accurate for the antenna alone. However, once integrated into the plane, the antenna diagram is perturbed. This is, for instance, particularly noticeable in the NRCS measurements in rotating mode, wherein a spurious azimuthal pattern could clearly be seen, and for fixed-antenna DFS observations, wherein a “striping” pattern as a function of incidence angle is obvious.

We have thus developed a complementary method that relies on the variations of the plane attitude during maneuvers. Using the plane IMU, we identify the angular coordinates α and β of the nadir and use the measured power to map the antenna diagram (using as a reference point the constant-altitude return power values for each data segment to account for geophysical nadir NRCS variations). The combination of all the flights during the campaign gives the distribution of measured power as a function of α and β, which is shown in Figs. B1 and B2.

The measured distribution is well approximated by a Gaussian shape

$$\begin{array}{}\text{(B1)}& G(\mathit{\alpha },\mathit{\beta })=\exp \left[-{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}}}-{\displaystyle \frac{{\left(\mathit{\beta }-{\mathit{\beta }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{\beta }}^{\mathrm{2}}}}\right].\end{array}$$

Another expression for G(α,β), more suitable for use with the half-power beamwidths α_−3 dB and β_−3 dB obtained from anechoic chamber measurements, is

$$\begin{array}{}\text{(B2)}& G(\mathit{\alpha },\mathit{\beta })={\mathrm{2}}^{\left[-\frac{\mathrm{4}{\mathit{\alpha }}^{\mathrm{2}}}{{\mathit{\alpha }}_{-\mathrm{3}\mathrm{dB}}^{\mathrm{2}}}-\frac{\mathrm{4}{\left(\mathit{\beta }-{\mathit{\beta }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{{\mathit{\beta }}_{-\mathrm{3}\mathrm{dB}}^{\mathrm{2}}}\right]}.\end{array}$$

The width parameters in these equations are linked by

$$\begin{array}{}\text{(B3)}& {\mathit{\sigma }}_{\mathit{\alpha }}={\mathit{\alpha }}_{-\mathrm{3}\mathrm{dB}}/\sqrt{\mathrm{8}\log \left(\mathrm{2}\right)},{\mathit{\sigma }}_{\mathit{\beta }}={\mathit{\beta }}_{-\mathrm{3}\mathrm{dB}}/\sqrt{\mathrm{8}\log \left(\mathrm{2}\right)}.\end{array}$$

The parameter values used in this study are collected in Table 1.

Figure B1 Reconstructed α and β dependence of the two-way KuROS antenna diagram. For each 30 s data segment, the constant-altitude values, for which the nadir is at α=0 and β=0, have been used as a reference level to account for geophysical variations in nadir NRCS.

Figure B2 Reconstructed azimuth dependence of nadir return power for different incidence angles. For each incidence angle, the α=0 value has been used as a reference. The thick line shows the final Gaussian fit used in the data analysis. The $\mathit{\beta }=\mathrm{15}{}^{\circ }$ data were excluded from the fit.

One cautionary remark is that the illuminated patch at nadir is not infinitely sharp. The measured distribution is thus the convolution of the true antenna diagram by the power distribution at the nadir patch (which depends on the altitude-tracking error and the sea state; Chelton et al., 1989). Assuming Gaussian shapes, the squares of the width parameters add, leading to

$$\begin{array}{}\text{(B4)}& {\mathit{\sigma }}_{\mathrm{observed}}\simeq {\mathit{\sigma }}_{\mathrm{true}}\left[\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{patch}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{true}}^{\mathrm{2}}}}\right].\end{array}$$

The broadening of the diagram due to finite nadir patch size is thus a small correction provided the scale of the nadir patch remains smaller than the antenna diagram scales. For reasonable orders of magnitude of the altitude-tracking error and significant wave height, the patch −3 dB width is of the order of 3^∘ when viewed from 3000 m of height. This corresponds to a 3 % correction on the value of σ_α. We have chosen to neglect this correction. The values summarized in Table 1 are the parameters of the Gaussian fits to the observed distributions.

B2 Rotating-antenna NRCS correction

Using these parameters as a starting point, we have then constructed corrections for the rotating-antenna measurements of NRCS by allowing the boresight elevation β₀ to vary as a function of antenna orientation within the plane. The variation law was determined by minimizing the dependence of the rotating-antenna NRCS measurements as a function of flight direction over the offshore area for each day.

B3 Fixed-antenna DFS correction

In a similar way, we have observed that the KuROS antenna diagram is slightly “wrinkled” in that the beam boresight azimuth changes as a function of elevation. This azimuthal mispointing transposes immediately into a striping modulation of the U_GD estimates. A correction was introduced by allowing the boresight azimuth α₀ to vary as a function of β. The variation law of α₀ was determined by minimizing the average U_GD over all flights for each value of β. As the variation of this quantity with respect to α₀(β) is not trivial, this required constructing, regularizing and inverting the observation matrix.

Appendix C: KaRADOC system

KaRADOC is built around an Agilent PNA-X network analyzer, complemented by a TX power amplifier, a T/R switch, an RX low-noise amplifier and a high-gain purpose-built slotted waveguide antenna (shown in Fig. C1).

The beam can be steered in elevation by changing the instrument working frequency (see Fig. C2a), and the antenna is usually mounted on a pitch–roll stabilization platform. For the Drift4SKIM experiment, however, the antenna was rigidly mounted in a port-looking configuration centered on a 10^∘ incidence angle with a 2^∘ backward-looking tilt to compensate for the aircraft pitch in constant-altitude flight. Plane attitude variations were accounted for in the data processing. Observations were collected at 33.7 GHz, corresponding to a 12^∘ nominal incidence angle. Other angles were also scanned, but RF leakage from the TX to the RX subsystems was too strong at the corresponding frequencies, making the signal harder to analyze.

The antenna radiation diagram is very narrow, with a beamwidth less than 1.5^∘ in elevation and less than 2^∘ in azimuth (see Fig. C2b). Figure C3a represents sections across the KaRADOC main lobe in the azimuth and elevation direction at 33.7 GHz.

Figure C1 The back (a) and front (b) of the antenna.

Figure C2 Frequency dependence of the KaRADOC main lobe azimuth and elevation boresight angles (a) as well as half-power beamwidths (b).

Figure C3 KaRADOC radiation diagram at 33.7 GHz as a function of (a) azimuth at 2.09^∘ elevation and (b) elevation at −0.05^∘ azimuth.