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 ﬁll this observation gap. The Sea tions 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 0 = 2 . 0ms − 1 at the Ka band, and more variable and close to C 0 = 2 . 4ms − 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 inﬂuence 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).

Abstract. 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 Kaband 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 veloc-ity 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 −1 , Hs 2.6 m) and gentle breeze on 24 November 2018 (wind speed 5.5 m s −1 , 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 . It is found to be closely aligned with the wind direction, with small correc-tions 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 0 = 2.0 m s −1 at the Ka band, and more variable and close to C 0 = 2.4 m s −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).

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(Hauser et al., , 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 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 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 −1 for low Earth orbit, it is obviously critical to have accurate knowledge of the measurement geom-etry 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  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 tidedominated 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 ; and validate the Radar Sensing Satellite Simulator (Nouguier, 2019) and its capability to simulate airborne configurations.
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 −1 for low Earth orbit and around 120 m s −1 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. The unit vector e ϕ is the projection on the horizontal of the line-ofsight 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.
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 dis-cussed 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 nongeophysical 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 groundbased 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, lineof-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 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.

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 −1 vs. 10 m s −1 for shipborne 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. 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 where σ α = α −3 dB / 8 log(2) and σ β = β −3 dB / 8 log(2). For 12 • observations the second term in the exponential can safely be neglected, and the effective azimuthal beamwidth can be estimated as 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 Figure 2. KuROS azimuth integral weight at θ = 12 • for a north-facing (ϕ b = 0 • ) antenna (black), Gaussian approximation (Eq. A29) (green) and variation of σ 0 for a typical 11 m s −1 wind from 140 • (dashed black). The peak of the σ 0 G 2 product (red) is shifted with respect to the peak of G 2 by δϕ −0.81 • . approximation in Eq. (A29) of G as a function of ϕ can then be used with the parameter Due to the width of the azimuthal aperture, the NRCSweighted 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 Denoting by ϕ t the flight track azimuth and V p the alongtrack flight velocity, the spurious azimuth gradient Doppler (AGD) contribution to the DV caused by the mispointing reads As an example, Fig. 2 shows the variations of the two-way antenna radiation diagram G 2 , its Gaussian approximation and the G 2 σ 0 (see Eq. A11) product as a function of azimuth at a 12 • incidence angle for a northward-looking KuROS antenna (ϕ b = 0 • ), using σ 0 data from the Drift4SKIM campaign on 22 November 2018. The effect of the wind-induced azimuthal gradient of σ 0 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 δϕ = ϕ a − ϕ b = −0.81 • = −15 × 10 −3 rad. For comparison (see Sect. 2.3 and Table 2), the pointing accuracy required to achieve a 15 cm s −1 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 3a shows a typical example of the azimuthal variation of σ 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 (ϕ = 320 • ), 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 p = 120 m s −1 and the Ku-band NRCS fit for both instruments (this is a reasonable assumption for orderof-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 −1 in some azimuth ranges.
Because the azimuth gradient U AGD contribution to the observed DV is proportional to V p σ 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 ∂ ϕ σ 0 /σ 0 are larger for larger incidence angles, this is more than compensated for by the 1/sin 2 θ reduction in azimuthal diversity across the footprint. This effect can thus be neglected for much higher incidence angles (Rodríguez et al., 2018).

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 backscatterweighted 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.
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 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 ρ(ξ , τ ) = η(x + ξ , t + τ )η(x, t) , the space-time covariance function of the sea surface elevation, and Q H = − 4π λ sin(θ )e ϕ and Q z = 4π λ cos(θ ), the horizontal and ver-tical components of the Ewald vector (with λ the radar wavelength), one obtains C(τ ) and ∂ τ C as The clear upwind-downwind asymmetry of σ 0 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) 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 where, in the notation of Nouguier et al. (2018), msv stands for the mean slope velocity and Mss for the mean square slope matrix, msv = mss xt mss yt , Mss = mss xx mss xy mss yx mss yy , with mss x α y β t γ = 2 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), where is the dispersion relation of gravity-capillary waves in deep water, with k the wave vector and κ M = 363.2 rad m −1 the wavenumber corresponding to the gravity-capillary regime transition. Introducing this expression in Eq. (10), and defining msv 0 as the spectral moment obtained using the dispersion relation of Eq. (15) in Eq. (13), one obtains the approximate DFS as and the corresponding V GD as While clearly oversimplified (it is, for instance, independent of the electromagnetic wavelength, which is known to have a significant influence on σ 0 ), 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 0 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 Mss −1 · msv 0 . In the rest of this article, M WD denotes the norm of this vector.
-As noted in Nouguier et al. (2018), msv 0 is equal to one-half the surface Stokes drift velocity of deepwater waves U ∞ S . 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 σ 0 as a function of incidence angle for different azimuths.
-In simple cases represented by parametric spectral forms such as the  spectrum used in this work, msv 0 and the eigenvectors of Mss are aligned with the downwind direction, and the V WD = −G D sin(θ ) e ϕ · U ∞ S relation proposed in Chapron et al. (2005) is recovered, with G D = 1 mss 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 2 for the msv 0 and Mss terms, respectively, while the parametric spectrum of , used in this work, decays as k −3 , leading to a logarithmic divergence for the Mss components and a slow convergence of the msv 0 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 0 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 lowwavenumber 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 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 Kaband 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 https://doi.org/10.5194/os-16-1399-2020 Ocean Sci., 16, 1399-1429, 2020  (Stopa et al., 2016), completed in the high-wavenumber range by  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 0 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 0 , with C 0 2.6 and C 0 2.2 m s −1 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 0 and Mss terms and can thus be expected to marginally change the numbers.

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 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 θ = 12 • .
-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 ψ = ϕ b − ϕ t , one obtains the consolidated error budget as 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.
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 −1 , ensuring a 0.4 m s −1 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-tracklooking 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.

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.

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 wavecurrent interactions (Ardhuin et al., 2009(Ardhuin et al., , 2012Guimaraes 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).
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 head-  ing 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).
-The R/V Thalia worked in the offshore area, providing continuous underway measurements of meteorolog-ical 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 conductivitytemperature-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 −3 . 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.
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.

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 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. 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 −1 from 140 • ). The 24 November flight, in contrast, took place during a steady weak southwesterly wind period (5 m s −1 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.

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 −1 , which is much larger than expected from the measurements (below aircraft speed of 120 m s −1 ). 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.

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 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 r = f − k · U/(2π ) for the moored Trèfle buoy. The red and blue arrows represent the AROME wind and MARS2D surface current vectors directions, respectively.  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 θ = 12 • 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.

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 σ 0 and DV data collected on 22 November 2018 when the wind speed was approximately 11 m s −1 , 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 geom- The 570 m scale bar applies to (f)-(i) and corresponds to the alongtrack 3 dB width of the radar beam at a 12 • incidence angle, i.e., near the middle of the swath. The mean trend of σ 0 as a function of θ has been removed from the σ 0 data. etry 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 Sutherland et al., 2018), as the swell propagates during the measurements at a phase speed of 22 m s −1 , while the aircraft moves at 120 m s −1 .
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 direc-tional 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.

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 rotatingantenna 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 σ 0 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 downwindcrosswind 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 σ 0 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 0 + a 1 cos(ϕ − ϕ σ,1 ) + a 2 cos[2(ϕ − ϕ σ,2 )]. 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 σ 0 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.

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-passfiltered 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. 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 −1 . 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 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). buoy data using the tools discussed in Sect. 2.2. The method combines the buoy spectrum up to 0.35 Hz and adds a highfrequency 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 −1 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 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 upcurrent 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 −1 instead of 11 m s −1 . 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 nearconstant value C 0 of the wave Doppler velocity magnitude, and this is the main result of the present paper.
The process leading to the estimation of the constant C 0 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 ) 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 Ocean Sci., 16, 1399-1429, 2020 https://doi.org/10.5194/os-16-1399-2020 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.
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 0 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 −1 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.

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 σ 0 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 MLMderived U WD vectors. The difference is clearly very large, reaching 2 m s −1 in places on 22 November, and smaller on 24 November as the azimuthal modulation of the radar NRCS is much weaker.
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  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 −1 on 22 November and of the order of 0.2 m s −1 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.

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 σ 0 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 σ 0 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 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. made 50 % wider in azimuth. With typical variations of σ 0 up to 1 dB over scales of the order of 1 km (e.g., Fig. 9), the variation of σ 0 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 −1 . This spurious velocity is larger than the 0.5 m s −1 significant orbital velocity of the swell. As a result the phase relation between DFS and σ 0 can change sign as a function of azimuth due to the combination of two imaging mechanisms with comparable magnitudes and possibly opposite signs. 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 σ 0 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.

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 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 wave-length of 30 m is typically enough for wind speeds higher than 7 m s −1 . 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 higherorder 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(ϕ − ϕ t )∂ ϕ log(σ 0 ) term in Eq. (6) of U AGD (as can be seen in Fig. 10, ∂ ϕ log(σ 0 ) is typically 0.1 rad −1 at a 12 • incidence angle). As the apparent mispointing due to σ 0 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.
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). 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 σ 2 ϕ V p /2 is the prefactor of sin(ϕ − ϕ t )∂ ϕ log(σ 0 ) in the expression of U AGD .

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 −1 at the Ka band and 2.4 m s −1 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 , 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 pencilbeam Doppler scatterometry exposed in Rodriguez (2018) and Rodríguez et al. (2018) to the case of near-nadir fanbeam 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 alongtrack 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) 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 = 2π/λ 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 ) 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 with σ 0 (t, x) the normalized radar backscattering cross section (NRCS) in the appropriate polarization and γ TS (|τ |) 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 1 and t 2 = t 1 + t, with t the pulse repetition interval (PRI), one obtains 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. Figure 1a and b summarize the acquisition geometry in the airborne and spaceborne settings. The antenna radiation diagram G 2 (t 1 , x) is represented as gray shading of the sea surface, while the range-point response function χ 2 (t 1 , r − r(t 1 , x)) is represented as white grating. In Eq. (A4), we have made the assumptions that G(t 1 , x) = G(t 2 , x) and χ (r − r(t 1 , x)) = χ (r − r(t 2 , x)), neglecting the effect of the spatial translation of the beam illumination pattern and rangeresolution 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 1 ) to x R (t 1 ) + 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 where the neglected terms are of the order of t the unit vector pointing from the radar location at t 1 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 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, the normalized weighting function, the average and fluctuating parts of the NRCS and borrowing the algebraic technique of "cumulant expansion" from probability theory, it is possible to express PP t as with κ n the successive "cumulants" of e(x) · (V R − v s (x)) with respect to the "density distribution" σ 0 (t 1 , r , x) W (t 1 , r , 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 As expected, the expression of κ 1 , 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 κ 2 , is a description of the impact of the variability of e(x), σ 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 the (geophysically relevant) weighted projection of the scatterer velocity in that frame on the radar line of sight, and Ocean Sci., 16, 1399-1429, 2020 https://doi.org/10.5194/os-16-1399-2020 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 Using Eq. (A13), one can obtain κ 1 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π , κ 1 is only obtained up to a multiple of λ 2 t . This effect can be neglected as long as both V NG and κ 1 remain within the unambiguous interval − λ 4 t ; λ 4 t . For larger platform velocities, care must be taken to add the right multiple of λ 2 t to κ 1 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 It is easy to see that V GD can be retrieved as 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 σ 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 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 V R = V R cos(g) cos(c)N + cos(g) sin(c)E + sin(g)D . (A22) In the NED frame, the pointing vector e can be expressed as e = sin(θ ) [cos(ϕ)N + sin(ϕ)E] + cos(θ )D. (A23) Its components in the platform-fixed frame can be determined using the fact that DCM −1 = 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 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 1 , r , x) and σ 0 (t 1 , r , x), one sees that at a 12 • incidence angle and for a platform velocity of 7000 m s −1 (spaceborne instrument), the SKIM 40 cm s −1 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 −1 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 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 1 , r , x) σ 0 (t 1 , r , 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 σ 0 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 θ = 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 surfacetracking algorithm to suffer from an error δh and detect the surface at range H − δr − δh. The data from this range bin will thus be processed using an angle of incidence different from the correct value by 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 surfacetracking 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 costconstrained 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 1 , r , ϕ) σ 0 (t 1 , r , ϕ) is well approximated by the Dirac distribution δ(ϕ − ϕ b ). In this limit, V NG (t 1 , r ) = V R cos(g) sin(θ ) cos(ϕ b − c) + sin(g) cos(θ ) . (A27) 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 σ 0 is to shift the peak of the distribution by an angle Assuming for W t 1 , t , ϕ a Gaussian approximation, in which σ ϕ (r ) is a parameter describing the width of the antenna diagram at the working incidence angle, one obtains V NG (t 1 , r ) = V R × cos(g) sin(θ ) cos(ϕ b − c + δϕ) + sin(g) cos(θ ) . (A30) with δϕ = σ 2 ϕ (r ) 2 ∂ ϕ log( σ 0 ).
Alternatively, one can choose to express V NG as the sum of V geo , the geometric approximation, plus an azimuth gradient Doppler velocity contribution, V NG t 1 , r = V geo (t 1 , r ) + V AGD t 1 , r , with V geo (t 1 , r ) = V R cos(g) sin(θ ) cos(ϕ b − c) + sin(g) cos(θ ) and V AGD (t 1 , r ) = − V R cos(g) sin(θ ) sin(ϕ b − c) 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 σ 0 and σ ϕ : At small scales, spatial gradients add to the azimuthal gradient and also induce a spurious velocity with the same expression as a function of σ 0 . Using the simple case of a single Fourier component σ 0 = ε sin [ν (ϕ − ϕ b )] allows one to evaluate the importance of different scales. The azimuthal shift can be obtained as In the slow variation limit (ν, σ ϕ → 0) and Eq. (A36) this expression coincides with Eq. (A31). For faster variations, one sees that the largest disturbance is obtained when ν ∼ √ 2/σ ϕ . 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 δϕ max = εσ ϕ e −1/2 / √ 2. 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 nongeophysical 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 α = 0 • , and the rotation axis of the antenna turntable points towards α = 0 • and β = 0 • . 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(α) = sin(θ ) sin(ϕ − ϕ b ) and tan(β) = tan(θ ) cos(ϕ − ϕ 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 constantaltitude 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 Another expression for G(α, β), more suitable for use with the half-power beamwidths α −3 dB and β −3 dB obtained from anechoic chamber measurements, is The width parameters in these equations are linked by σ α = α −3 dB / 8 log(2), σ β = β −3 dB / 8 log(2). (B3) Figure B1. Reconstructed α and β dependence of the two-way KuROS antenna diagram. For each 30 s data segment, the constantaltitude 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 β = 15 • data were excluded from the fit.
The parameter values used in this study are collected in Table 1.
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 altitudetracking error and the sea state; Chelton et al., 1989). Assuming Gaussian shapes, the squares of the width parameters add, leading to σ observed σ true 1 + σ 2 patch 2σ 2 true . (B4) 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 β 0 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 α 0 to vary as a function of β. The variation law of α 0 was determined by minimizing the average U GD over all flights for each value of β. As the variation of this quantity with respect to α 0 (β) 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.