Measuring ocean surface velocities with the KuROS and KaRADOC airborne near-nadir Doppler radars: a multi-scale analysis in preparation of the SKIM mission

Surface currents are poorly known over most of the oceans. Satellite-borne Doppler Waves and Current Scatterometers (DWCS) can be used to fill this observation gap. The Sea surface KInematics Multiscale (SKIM) proposal, is the first satellite concept built on a DWCS design at near-nadir angles, and now one of the two candidates to become the 9th mission of the European Space Agency Earth Explorer program. As part of the detailed design and feasibility studies (phase A) funded by ESA, airborne measurements were carried out with both a Ku-Band and a Ka-Band Doppler radars looking at the sea surface at 5 near nadir-incidence in a real-aperture mode, i.e. in a geometry and mode similar to that of SKIM. The airborne radar KuROS was deployed to provide simultaneous measurements of the radar backscatter and Doppler velocity, in a side-looking configuration, with an 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 KaRADOC system has a much narrower beam, with a circular footprint only 45 m in diameter. 10 The experiment took place in November 2018 off the French Atlantic coast, with sea states representative of the open ocean and a well known tide-dominated current regime. The data set is analyzed to explore the contribution of non-geophysical velocities to the measurement and how the geophysical part of the measured velocity combines wave-resolved and waveaveraged scales. We find that the measured Doppler velocity contains a characteristic wave phase speed, called here C0 that is analogous to the Bragg phase speed of coastal High Frequency radars that use a grazing measurement geometry, with little 15 variations ∆C associated to changes in sea state. The Ka-band measurements at an incidence of 12◦ are 10% lower than the theoretical estimate C0 ' 2.4 m/s for typical oceanic conditions defined by a wind speed of 7 m/s and a significant wave height of 2 m. For Ku-band the measured data is 1 https://doi.org/10.5194/os-2019-77 Preprint. Discussion started: 27 August 2019 c © Author(s) 2019. CC BY 4.0 License.

with the SKIM viewing geometry. A variation of surface backscatter across the footprint and as a function of azimuth ϕ is represented by the grey shading, and gives an effective mispointing δ. In KuROS data, each measurement is integrated in azimuth across the antenna lobe.
In the case of SKIM, the use of unfocused SAR allows the separation of echoes in the azimuth direction with a resolution dDop 300 m. ϕ, is the sum of a horizontal geophysical Doppler velocity contribution U GD (ϕ), and a non-geophysical velocity V NG . The following measurement equation is given by projections of the target and sensor velocity vectors onto the line of sight as shown in figure 1, (2) 2.1 Non-geophysical velocity V NG

5
In practice, V NG is the radar velocity projected onto the effective look direction, that includes an apparent azimuth mispointing includes both spatial gradients and azimuthal gradients. As a result, the beamwidth as a function of incidence is a very important parameter of the radar.  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. In the case of level flight and low incidence observations, one can however obtain a Gaussian approximation to the 1-way radiation diagram as: G exp − ϕ 2 2 sin 2 (θ) σ 2 α + (β 0 − tan(θ)) tan(θ) where σ α = α −3dB / 8 log(2), σ β = β −3dB / 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, ϕ −3dB is thus larger than α −3dB by a factor 1/ sin(θ), equal to 4.8 for 12°measurements.
Provided that the beam does not grow too wide, the Gaussian approximation eq. (A27) of G as a function of ϕ can thus be used, with parameter 5 σ ϕ α −3dB / sin θ 8 log(2) .
When the sea surface NRCS is variable, this finite radar aperture gives an apparent mispointing that is the difference between the apparent azimuth ϕ a and boresight azimuth ϕ b , as detailed in Appendix A If not corrected for, this apparent mispointing gives a spurious velocity that is the projection of the platform velocity onto the 10 apparent line of sight where V p is the along-track velocity of the platform carrying the radar in the frame of reference of the solid Earth.
For our experimental KuROS configuration σ α 6.36 • = 0.11 rad at θ = 12 • . As shown in section 3, ∂ ϕ σ 0 /σ 0 is of the order of 0.10 rad −1 for a uniform wind speed of 11 m/s. This gives an apparent mispointing of the order of δ = 0.8 • = 14×10 −3 rad, which would correspond to a 1.7 m/s error on U GD . Here, it is important to note that KuROS was not specifically designed for this experiment, but primarily as a Calibration/Validation instrument for the CFOSAT mission, which required a 5 broad radiation diagram. Though the analysis of the KuROS data helped uncover many interesting effects relevant to Doppler observations of the sea surface (for instance on the speckle noise dependence on observation direction) 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.
As detailed in Appendix A, eq. (4) only applies for a narrow beam when projected on the ground, which is not a very 10 good approximation for the KuROS case, even at an incidence of 12 • . As shown in figure 2, the Gaussian approximation for the antenna pattern as a function of ϕ gives a too narrow distribution and does not take properly into account the azimuthal integration, leading to an overestimation of U AGD .
As an example, figure 2 shows the variations of the two-way antenna radiation diagram G 2 , of its Gaussian approximation, and of the G 2 σ 0 product as a function of azimuth at 12 • incidence, for a northward-looking KuROS antenna (ϕ b = 0 • ), using 20 Figure 3.A shows a typical azimuthal variation of σ 0 for an incidence of 12 • using a Ku-band radar. As expected for nearnadir measurements (Chapron et al., 2002;Munk, 2008), the NRCS is largest in the downwind look direction (ϕ = −40 • ), has a secondary peak in the upwind direction, and is weakest in the crosswind look directions. Figure 3.B shows the expected spurious contribution U AGD to the geophysical velocity U GD , if the apparent mispointing δ is not corrected for. This uses an aircraft velocity V p = 120 m s −1 , for the KuROS and KaRADOC cases. The Ku-band NRCS fit has been used for the Ka-band 25 instrument as well. This is a reasonable assumption for order-of-magnitude estimates.
Because both the azimuth gradient Doppler U AGD and the spatial gradient Doppler U SGD are proportional to V p σ 2 ϕ it is clear that the broad KuROS antenna pattern requires a very accurate estimation to correct for U AGD , which is almost negligible for KaRADOC or DopplerScatt (Rodríguez et al., 2018), thanks to their narrow azimuthal beam aperture. Another remark is that the approximate expression eq. (A35), though it gives the appropriate dependency of U AGD with respect to look azimuth, tends 30 to over-predict its magnitude, as the widening associated to the ground projection saturates for broad beams.
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 . KuROS azimuth integral weight at θ = 12 • , north-facing antenna (black), Gaussian approximation (eq. A27) (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 • .
importance of different scales. The azimuthal shift can be obtained as In the slow-variation limit ν, σ ϕ → 0, and eq. (6) this expression coincides with eq. (A30). 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 5 the brightness contrast ε. δ in this case is given by  Although the relative variations ∂σ 0 (ϕ)/∂ϕ/σ 0 are larger for larger incidence angles, this is more than compensated by the 1/ sin 2 θ reduction in azimuthal diversity across the footprint. This is why this effect can be neglected for much higher incidence angles (Rodríguez et al., 2018).

Geophysical velocity U GD : Waves and Current Doppler
The geophysical part of the Doppler shift measured by a microwave radar over the ocean, using both Along-Track-Interferometry and Doppler centroid techniques is caused by the backscatter-weighted average of the surface velocities along the line of sight, as illustrated in figure 4.
For a perfect sine wave of period T propagating over deep water, the phase speed of the wave is where U is the current speed, ϕ w is the wave propagation azimuth direction, and ϕ U is the current direction. Measuring the phase speed deviation from the theoretical value is the principle of the coastal HF radars (Barrick et al., 1974;Stewart and Joy, 1974), for which the grazing angles coherent Bragg back-scattering mechanism selects very effectively the sine wave components of the sea state which interact with the radio waves.  In the case of KuROS, SKIM, or other systems for which the backscattering mechanism is not selective, a superposition of waves contributes. We get a compound mixture of C's, and the different terms become Yurovsky et al., 2019), where the current U sampled at depth k/4π (Stewart and Joy, 1974) for a monochromatic sine wave, is now U CD , a weighted average of the currents at different depths, where the weighting function is determined by the wave slope spectrum. Indeed, the velocity U GD integrates the velocity of the tilted-facets on the ocean surface. Within facets, quasi-specular and specular points are selected, further modulated by the local directional tilts ∇η, leading to a modulated averaged intensity and a weighted 10 velocity as with σ(∇η) an individual local radar cross section, corresponding to a tilted facet, with a local line-of-sight velocity V LOS , and P (∇η) the facet-tilt probability distribution.

15
As the incidence angle is increased beyond 25 • , the backscatter is dominated by Bragg scattering and the phase speeds that contribute are the Doppler-shifted Bragg waves (Rodríguez et al., 2018). In our case, Bragg scattering is generally negligible except for the lowest wind speeds, and the Ka-band resonant Bragg scattering scale at 12 • is about 2 cm, around the capillarygravity wave transition, corresponding to the minimum phase velocity of about 23 cm/s. A general analysis valid for all incidence angles is presented by Yurovsky et al. (2019). Here we focus on incidence angles from 4 to 20 • where the backscatter 20 modulation is dominated by tilt effects (Kudryavtsev et al., 2017).
The correlation of surface slope and line-of-sight velocity defines the mean slope velocity in direction ϕ, msv(ϕ) . For linear ocean waves, this equals the correlation of vertical velocity gradients and displacements, equal to half the Stokes drift in direction ϕ, U S cos(ϕ−ϕ S ). The surface Stokes drift magnitude U S and direction ϕ S can be computed from wave buoy measurements (Kenyon, 1969;Ardhuin et al., 2009). Even though U S is highly correlated with the wind speed, with 25 a Pearson's linear correlation coefficient of 0.85 or so, it has a strong variability with the sea state as illustrated in Fig. 5.a.
Using a Kirchoff approximation, the variation of wave Doppler with observation direction ϕ can be computed from the wave spectrum (eq. 16 in Nouguier et al., 2018, see also Appendix C), where U S is the magnitude of the surface Stokes drift vector, which has a direction ϕ S , and the direction ϕ WD is found to be 30 within a few degrees of ϕ S . The Doppler imaging factor G D is a weakly varying function of the radar frequency and incidence angle, but also of sea state ) For Ka-band, we may use a typical variation of σ 0 of the form, (Walsh et al., 2008;Nouguier et al., 2018), where non-Gaussian corrections are, B 0.5676A 1.332 and the modification of σ 0 max related to the mss shape of Nouguier et al.

5
(2016), and A = 1/mss shape For Ku-band the mss shape is generally smaller than the Ka-band value, in particular for wind speeds over 5 m/s, as shown in Average values of the mss shape give typical values G D 25, slowly decreasing with increasing θ. Therefore, both Stokes 10 drift magnitude and G D grow with wind speed. It may thus be more practical to express U WD in the following form where M WD varies very little for most of the sea state conditions, around C 0 = 2.2 m/s in Ka-band and C 0 = 2.4 m/s in Ku band. In other words, most of the variability of U WD is controlled by the directionality effect and the magnitude M WD is a 15 weakly varying function of the wave age and presence of swell (see also Yurovsky et al., 2019;LOPS, 2019).
3 Campaign overview

Organization of the campaign
Given our objectives of demonstrating the sensitivity of airborne Doppler measurements to the geophysical contributions of currents and waves, it was important to have commonly accepted reference measurements for these parameters. Also, in order 20 to verify the limited effect of wave development on the geophysical Doppler velocity U GD , we decided to go to an oceanic environment, open to offshore swells, which makes our experiment different from previous Doppler airborne campaigns (Martin et al., 2016;Rodríguez et al., 2018).
Field work was focused in two "boxes" (Figure 6) named the "offshore box", around the "Trefle" buoy (see below), and the "Keller race box" to the north of the Ushant island. Both boxes are in the range of coverage of a WERA-type High-Frequency 25 radar (Gurgel et al., 1999;Gurgel and Barbin, 2008), operated by Service Hydrographique et Oceanographique de la Marine (Shom) and already used for several studies, in particular related to wave-current interactions (Ardhuin et al., 2009(Ardhuin et al., , 2012Guimaraes et al., 2018). The "Keller race box" is an area with very strong horizontal gradients of the current. Although it is easy to show a strong effect of the current on a measured Doppler, 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 box", 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 box" are presented in this paper.

5
The week around spring tides of November 2018 was targeted, in order to allow for a wide range of current speeds (Fig. 7a), as well as to accommodate plane availability constraints.
Ground truth measurements comprised two permanent operational systems: the 12 MHz WERA-type HF radar mentioned previously, with expected depth of measurement around 1 m (Stewart and Joy, 1974), and a wave-measuring buoy 'Pierres Noires', also known by its World Meteorological Organization number 62069.

10
Dedicated instrumentation was also deployed for the campaign: -the "Trefle" buoy was moored at 5°25' W, 48°25' N, in the middle of the offshore box. This buoy monitored the surface motion (Sutherland et al., 2016) and provided directional wave spectra (Fig. 8).

15
-the R/V Thalia worked in the offshore box, provided continuous underway measurements of meteorological parameters using a Météo-France "BATOS" operational system comprising a Vaisala WXT-series sonic anemometer located approximately 10 m above sea surface. The ship also carries a SBE21 thermosalinograph.
In the summer, the so-called "Ushant tidal front" (Le Boyer et al., 2009) has a strong influence on the current and conditions in the offshore box. CTD casts were performed from R/V Thalia during the campaign, that showed the water column to be very 20 well mixed, 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 second plane (a Piper PA-23 also operated by SAFIRE) which surveyed the offshore box in a "lawn mowing" pattern, flying under the clouds from an altitude of 500 m to 1000 m. Interesting small-scale surface signatures could be observed on calm days, but it is clear that no density-associated mesoscale structure was present.

25
The aircraft measurement geometry over the "offshore box" consists of relatively long (12 km) and straight tracks with differ-   on 11/26 from 09:40 to 11:00. In this paper, we focus on data acquired on 11/22 and 11/24 as the geophysical conditions were interesting and complementary (see below) and data were acquired with the largest azimuth diversity on these two days.

Geophysical conditions
The November 22 flight took place at the end of an interesting steady southeasterly wind episode (13 m/s from 140°). The November 24 flight in contrast took place during a steady weak south-westerly wind period (5 m/s from 270°) (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 to 0.9 m from the 21 to 24 November, with a peak frequency increasing from 0.07 Hz to 0.1 Hz, and a mean direction gradually veering from west to north-west. This swell contributes little to the Stokes drift (about 10% of the windsea contribution on 22/11).
The main environmental conditions at the time of these star-pattern flights with a fixed antenna are summarized in Table 2.  (Lygre and Krogstad, 1986), and Doppler shifted with fr = f − k · U/(2π) for the moored Trefle buoy.

KuROS set-up
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 pre-launch studies. Of the two antennas, the Low Incidence antenna is nominally centered on 14°incidence, while the Medium Incidence antenna is nominally centered on 40°incidence. Only the LI antenna, which was the more relevant for SKIM, was used during the campaign. This antenna uses a HH polarization. A 5 comprehensive description of the system can be found in Caudal et al. (2014), with an important modification consisting of a new antenna with characteristics given in Table 1.
The radar transmits a frequency-modulated pulse (chirp) with a 100 MHz frequency band, achieving a 1.5 m range resolution and an effective resolution around 7 m in elevation once projected on the ground (at 12°). The 1-way 3-dB footprint in azimuth is 580 m wide at 12°and 3000 m 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 on the Doppler velocity measurement (see section A2.1 in the appendices) is about 126 m/s, which is much larger than expected from the measurements (below aircraft speed of 120 m/s).
As discussed in Appendix A, accuracy requirements on observation geometry are much less stringent for cross-track than for up/down-track Doppler observations. The Doppler 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 sensor
The Ka-band RADar for Ocean Current monitoring (KaRADOC) airborne radar sensor has been 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 Ka-band. Further details on the system are given in Appendix C.
KaRADOC was mounted under the ATR42 aircraft, in a port-looking configuration. The two-way 3 dB footprint from 5 3000 m 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 incidence is selected by varying the working frequency. Data was acquired at different incidence angles, from 6 to 14 • , corresponding to a change of frequency from 32.5 to 38.2 GHz. Here we only report on θ = 12 • observations. This radar does not incorporate a range-resolution scheme: the transmitted pulses last several µs, 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 have been varied during the acquisitions. Though they have a strong impact on NRCS estimates quality and Doppler estimates noise, we have found the low-pass filtered Doppler signal to be robust.

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

Measurements
4.1 Qualitative description of the KuROS data 10 We will start with the σ 0 data because it is a more generally understood property of the sea surface and the analysis of its variation in azimuth ϕ and elevation θ is necessary to interpret the Doppler mean value and modulations. However, it is interesting to first look at the rawest available data, which can be displayed as images of backscatter and Doppler. Figure 9 shows pieces of the first few tracks acquired on November 22, when the wind speed was around 11 m/s, after removing the aircraft velocity (see Appendix B and C for details).

15
The first remark is that the back-scatter, after correcting for an incidence-varying mean trend, is smooth with a typical amplitude of variation of 1 dB ( Fig. 9.A). This smoothness comes from the large footprint and radiometric quality of the data.
The KuROS data clearly reveals the presence of the swell from the west, with a peak frequency of 0.07 Hz corresponding to a wavelength L = 320 m. This is particularly visible in the north-south oriented flight tracks number 6 and 5 in Fig. 9.B and Fig. 9.E-F for a zoom on track 6. The apparent swell direction (dashed lines in figure 9.B) differs from its true direction due to 20 scanning distortion Sutherland and Brozena, 2018), as the swell propagates during the measurement at a .E-F for a zoom on track 6. The "long-crested" appearance of short waves in (D) and (E) is an artefact of the wavefrontmatching observation (Jackson et al., 1985), with all other directions averaged out by the large azimuth width of the radar beam. If purely geophysical, the velocity associated to the Doppler shifts 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 30 eastward velocities toward the radar (blue). This will be discussed in further details below. Finally (H) and (I) exhibit chevron patterns with crests facing both north-east and nort-west directions, whereas the waves from the south-west are expected to be much longer than those from the south-east, which is not apparent in the KuROS data.

Ku-band backscatter
Given the expected large influence of the wind, wave and current directions, and the high noise on the aircraft measurement geometry estimation, data have been averaged along 30-s track portions. Data have been averaged in incidence angle (across track) every 1 • with one point every 30-s. As highlighted before, full tracks are straight and relatively long (12 km) in a homogeneous ocean region. Each track has a different aircraft heading, forming a "star" pattern, sampling a wide range of  The variation of σ 0 in Ku-band is shown for 22 November in Fig. 10, with one dot for each 30-s long record. These measurements show the expected azimuthal modulation of 0.8 to 0.9 dB with a downwind-crosswind contrast that increases from low to high incidence angle. This contrast is larger for the higher winds of 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). Discarding these azimuths (shaded in grey in figure 10.C), we fitted a functional variation of the form a 0 +a 1 cos(ϕ−ϕ σ,1 )+ a 2 cos[2(ϕ − ϕ σ,2 )]. As explained in section 2.1, this azimuthal variation is critical for the interpretation of the mean Doppler 5 due to the spurious azimuth gradient Doppler velocity. 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 November 24, the σ 0 (Fig. 11) was much more uniform with azimuths, due to the much lower wind speed.

Mean Doppler from KaRADOC
We now discuss quantitatively the measured Doppler signal, in order to assess the agreement of our theory of the wave-induced 10 Doppler U WD with the measurements.
Starting with KaRADOC, easier to interpret than KuROS thanks to its narrower beam, we first examine the consistency of all the data acquired at an incidence of 12 • . Figure   Looking in details, we find that the model gives larger values of U GD than what is observed for both days. An ad hoc reduction of U WD by 10% gives the best agreement between model and data. With such a reduction, the surface current vector 10 is accurately inverted from the data, when U WD is subtracted off the fitted U GD , as shown in figure 14.
The reason for the 10% model overestimation is not clear, and it may be due to a particular processing of the airborne data or assumptions in the model. This second hypothesis was tested by varying the method for estimating U WD . Table 3 gives a subset of model tests with varying the exact input spectrum.
The largest differences, of the order of 10%, were found when changing the transition frequency from 0.35 to 0.5 Hz 15 with a general increase of the simulated U WD values. This can be interpreted as the effect of the directional distribution in the Elfouhaily (1997) spectrum that is probably slightly too narrow for intermediate wavelengths 2-10 m. This effect was discussed for very specific cases by Peureux et al. (2018) and it is not yet clear if it is specific to the very young wind seas they observed, although it could explain some properties of L-band backscatter (Yueh et al., 2013). Other effects, in particular the non-linearity of the waves (e.g. Nouguier et al., 2009Nouguier et al., , 2015 may also contribute a few percent to the deviation of mean Doppler.

20
While testing different model settings for U WD , we also estimated the importance of the swell. Although the swell has a limited contribution to the mean square slope, in the present case it contains 6% of the zonal Stokes drift component V s . When projected on the Doppler velocity direction, this explains the increase by 4% of U WD when the swell is removed in table 3 because the swell contribution is in the direction opposite to the wind sea and thus reduces U s and U WD (see also Yurovsky et al., 2019). This is more easily understood by looking at the contribution of the different spectral components to U s . In 25 practice, U s (ϕ) is the projection of the horizontal vector (U s , V s ) in direction ϕ and can be obtained from the heave spectrum E(k) of the sea surface and the first directional moments a 1 (k) and b 1 (k), or the associated mean direction ϕ 1 (k), which is the direction of the vector (a 1 (k), b 1 (k)) and spread σ 1 (k) (Kuik et al., 1988), that can be measured from drifting wave buoys, Using the dispersion relation for linear surface gravity waves in deep water, σ 2 = gk with g the acceleration of gravity, the 30 Stokes drift contribution at each wavenumber is given by the following spectrum (Kenyon, 1969),  Figure 14. Retrieval of surface current vector UCD in blue obtained by reducing UWD by 10%, compared to in situ measurements by HF-radar, CARTHE and SVP drifter (three shades of grey). which integrates to Table 3. Modeled wave Doppler amplitude MWD and direction ϕWD using the same input spectra but varying the transition frequency ft between the buoy spectrum and a Elfouhaily (1997) spectrum for the high frequencies, and the wave age Ω for this high frequency spectrum.
Removing the swell from the input spectrum was also tested.  Figure 15 shows example of these spectra in the case of November 22, with a 2 m swell from the west, giving positive contributions to F x for frequencies under 0.1 Hz (wavelengths larger than 150 m), and a wind-sea from the south-east, giving positive contributions to F y and negative contributions to F x .
Although the wind-sea energy is much less than that of the swell, its contribution, weighted by k 1.5 ∝ f 3 is much larger. a bi-harmonic form of σ 0 fitted to the measurements (Fig. 10) to compute U AGD . As shown in Fig. 16.A,B, removing U AGD reduces the maximum average velocities (from 4 to 3 m/s) that correspond to the saturated red patches in Fig. 9.D.

Mean Doppler from KuROS
Still the measured Doppler velocities do not fit very well the forward model, in particular for the azimuths 210 to 300 degrees. The fit is still poor but relatively better at higher incidence angles, with 18 • shown in Fig. 16.C. Given the lower importance of antenna lobe effects for higher incidences, it is quite possible that the mean Doppler deviations from the model Another issue is the general magnitude of the measured U GD that is smaller with KuROS compared to KaRADOC, whereas theory would suggest the contrary. In particular, the wave Doppler, in a first approximation is ratio of the wave Stokes drift and the effective mean square slope mss shape as defined by Nouguier et al. (2016). With the known reduction of mss shape from Ka to Ku band, this should lead to a higher value of U WD and thus U GD . The KuROS data thus appear anomalous.
As a result, any attempt to retrieve surface current vectors from KuROS is not very successful for most incidence angles, as The amplitude of the Doppler velocity are not much reduced, in spite of a more than halved current and Stokes drift. This is consistent with the expected near-constant value C 0 of the wave Doppler magnitude. We also note that the upwind/downwind asymmetry of the Doppler is much reduced.
A simplified vector interpretation of the measured Doppler is shown in Fig. 18. It shows large differences, of the order of 50 10 to 100 cm/s between the KuROS estimates and the in situ measurements.

Observed Doppler modulations
In spite of issues related to the insufficient knowledge of the antenna pattern for a retrieval of the mean Doppler velocity, the KuROS measurement are very interesting for the analysis of the small scale velocity gradients. In particular, Caudal et al. (2014), with a different antenna (slightly narrower beam) had successfully used the cross-spectra of Doppler and NRCS to 15 resolve the 180 • ambiguity in the wave propagation direction. Also, an alternative to the theoretical model for U WD used above, it may be possible to use the resolved part of the correlation between σ 0 and Doppler to estimate the unresolved part and the full U GD .
In practice the Doppler modulations are also caused by the gradients of σ 0 and the speed of the aircraft, just like the mean spurious U AGD velocity. These spurious modulations are enhanced by 70% when the antenna pattern is made 50% wider in 20 azimuth, as illustrated in Figure 19. 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   Figure 19. Qualitative validation of radar imaging mechanism in R3S simulations (Nouguier, 2019). Both the real data and simulation contain the geophysical modulation of velocities associated to surface velocities and slopes in the look direction (part of UGD ) and aircraft velocities and slopes in the flight direction (part of UAGD). Note that the wave phases in the R3S simulation are random and cannot be expected to match those in the data or between the two simulations. This effect will be weaker for shorter (wind sea) components as soon as the wavelength and crest length becomes much shorter than the KuROS footprint L y , as given by eq. (6): 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 radar, in Ka and Ku band, using the same pulse-pair technique but very different antenna patterns has 5 provided important insight for the preparation of the SKIM mission.
In terms of radar measurement, the DRIF4SKIM campaign clearly demonstrated the robustness of the retrieval approach proposed for SKIM ESA, 2019), with the use of detailed wave spectrum measurements to estimate the wave Doppler contribution U WD , associated to the wave intrinsic phase speed sensed through the surface slope spectrum.
Accurate wave spectral measurements of the first directional moments, such as provided by a buoy, is sufficient to estimate U WD 10 and resolving wavelengths of 15 m (a frequency of 0.32 Hz) is generally sufficient to estimate the full spectral contribution, appending a parametric spectral shape of the unresolved shorter waves. In fact it is most important to resolve the peak of the windsea, and a resolved wavelength of 30 m is typically enough for wind speeds higher than 7 m/s. Subtracting off U WD from the measured geophysical Doppler provides an estimate of the total surface current velocity U CD .
This velocity is expected to be representative of the top meter of the ocean as given by the depth of measurement for each monochromatic wave train (Stewart and Joy, 1974), weighted by the slope spectrum. This velocity should contain most of the Stokes drift (Barrick and Weber, 1977;Broche et al., 1983;Ardhuin et al., 2009). As a result the velocity from any near-nadir Doppler Wave and Current Scatterometer, such as SKIM, is expected to give values between those of CARTHE drifters and a 5 12 MHz HF radar, consistent with the results shown in figure 18.
The campaign also stressed the necessity of a very good knowledge of the measurement geometry, including antenna pattern, and the spatial and azimuthal variation of the radar cross-section. In this respect, the main characteristics of the present campaign and planned SKIM satellite mission are recalled in table 4. Table 4. Main differences between airborne KaRADOC and KuROS system as used in the present paper, and the SKIM system as presented by ESA (2019). The factor σ 2 ϕ sin θVp/2 is the common factor that, multiplied by sin ϕ∂(log(σ 0 ))/∂ϕ gives the spurious velocity UAGD, given by eq. (A35) A true mispointing in azimuth δ ϕ gives a spurious velocity V p δ ϕ sin ϕ, and aircraft velocities lead to much less strict require-10 ments than satellite systems for which a pointing accuracy of a few microradians cannot be achieved by attitude measurements alone (gyroscopes and star-trackers) but uses 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 pattern (Rodríguez et al., 2018, see also Appendix D).
The apparent mispointing due to σ 0 gradients in azimuth or space is proportional to the beam width squared, and leads to 15 larger non-geophysical velocities for KuROS at 12 • incidence, than for SKIM, even at 6 • . That effect however is practically negligible for KaRADOC at 12 • , as shown in figure 2.
Finally we recall that the incidence angle is estimated from the range measurement in the case of KuROS and from the highly directive antenna pattern for KaRADOC. Translated to SKIM, which uses range measurements like KuROS, it requires a knowledge of the local slope of the ocean as the nearest ocean target for the nadir beam is not exactly at nadir but slightly 20 offset due to the geoid slope, up to 300 microradians (Sandwell and Smith, 2014) and more importantly the sea surface height difference between the nadir and the oblique beams.
Other radar system constraints or optimizations for satellite systems are discussed by Rodriguez (2018) and (ESA, 2019, chapter 5), with sampling issues further analyzed by Chelton et al. (2019).
Although it was the first deployment of the KaRADOC instrument and only a limited dataset could be acquired due to the necessary adjustment process, the DRIF4SKIM campaign clearly demonstrated that surface geophysical velocities can be measured by a pulse-pair method with a Ka-band measuring at 12 • incidence. The campaign data are consistent, with a 10% bias, with a Geophysical Model Function (GMF) that expresses the geophysical Doppler as the sum of the range component For 11 m/s winds, the difference in drift velocity between CARTHE and SVP drifters and HF radar, of the order of 20 cm/s, is not due to a large scale vertical shear, and may be associated to either a particular behaviour of CARTHE drifters in wave 10 motions, or some shear in the top meter. This will be further investigated elsewhere using the 30 cm depth measurements of the Trèfle buoy (which capsized on November 21 in the present experiment).
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.
Finally, the test of near-nadir satellite measurements is limited by the very different viewing geometry due to the difference 15 in altitude. Airborne measurement footprints are at most 500 m or so, and thus cannot reproduce the averaging properties of a much wider footprint from a satellite. 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. This limited azimuthal resolution is probably sufficient for estimating parameters such as the Stokes drift given by eq. (22) and does not require a full turn of the antenna to observe all waves, thereby 20 making possible a higher spatial resolution of the Stokes drift vector.
Future airborne systems may ideally combine higher incidence angles such as used on DopplerScatt (Rodríguez et al., 2018) and OSCAR/Wavemill , with near-nadir angles that allow unambiguous wave measurements. In that case, the large azimutal footprint of KuROS is probably not necessary, and a narrower beam like KaRADOC can be used, greatly simplifying the analysis of antenna beam patterns.

25
Code and data availability. Data and numerical model results presented in this article are available via ftp at the following address: http://tinyurl.com/SKIMftp for a single pulse is random and uniformly distributed over the unit circle. The radar returns of successive pulses transmitted at short intervals are however correlated, and the time history of the phase can be used to measure the relative motion of the radar and the scatterers. SARs make use of this property to refine the along-track resolution of backscattering cross-section 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 10 can be expressed as where the integral is performed over the sea surface; A(r ) is a time-independent weakly-dependent function of range, unim-15 portant 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 pattern; χ(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.

20
As mentioned by Rodríguez et al. (2018), the thermal noise contribution, though it plays a major role in setting the quality of the measurements, poses no great conceptual difficulty, and can be safely considered as δ-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 a much richer physics. It is affected by electromagnetic phenomena, by the geometry and kinematics of the sea surface itself, and its statistics are further complicated by the so-called "speckle" 25 phenomenon. As stated by Rodríguez et al. (2018) the correlation function of this coefficient as a function of time and space separation, averaging over speckle realizations, can be modelled as with σ 0 (t, x) the Normalized Radar backscattering Cross Section (NRCS) in the appropriate polarization, and γ T S (|τ |) a function describing its time decorrelation at a fixed location, due to the life history of individual scattering patches. 30 Combining expressions (A1) and (A2) to compute the speckle-averaged product of the return signals for two consecutive radar pulses sent at t 1 and t 2 = t 1 + ∆t, with ∆t the pulse repetition interval (PRI), one obtains P P ∆t (t 1 , r ) = E 2 (t 2 = t 1 + ∆t, r )E 1 (t 1 , r ) * as The phase of eq. (A3) 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 and this method is the so-called "Pulse-Pair" technique of Zrnic (1977).

A2 Measurement geometry
Figures 1A and 1B summarize the acquisition geometry in the airborne and space-borne settings. The influence of the antenna 10 radiation diagram G 2 (t 1 , x) is represented as a grey shading of the sea surface, while the influence of the range point-response function χ 2 (t 1 , r − r(t 1 , x)) is represented as the grating in white. In eq. (A3), 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 range-resolution weighting distribution on the sea surface.
This is a very good approximation for airborne pulse-pair radar observations, and a quite good one for spaceborne obser-15 vations. 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 wavelength to avoid phase ambiguity. For space-borne 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 as well as of the range point-response.

A2.1 Pulse-pair signal approximation
Returning to expression (A3), 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 scatterers originally located at x have moved to x + v s ∆t (specifying the reference frame is not yet necessary since only relative separations are important at this stage). The radar-to-scatterers vector has thus changed by [v s (x) − V P ] ∆t. The distance change can be approximated by where neglected are of the order of ∆t 2 ||v s − v R || 2 /||x − x R (t 1 )|| 2 . Introducing 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 exponential contributions to an integral. Obtaining an equivalent representation as the exponential of a sum of weighted integrals would be desirable. Introducing the effective illuminated surface 10 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 P P ∆t as with κ n the successive "cumulants" of e(x) · (V R − v s (x)) with respect to the "density distribution" σ 0 (t 1 , x)W (t 1 , r , x). As 20 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 , 25 κ 1 (t 1 , r ) = (A13) 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.
A2.2 Pulse-pair signal phase approximation 5 Selecting a frame of reference, we define the (geophysically relevant) weighted projection of the scatterers velocity in that frame on the radar line-of-sight and the (geophysically irrelevant) projection of the radar velocity. (Our conventions are such that V GD is positive when the scatterers 10 move towards the radar, and that V N G is positive when the radar moves towards the footprint, in keeping with everyday intuition).
With these conventions, one sees that: Using equation (A12), one can obtain κ 1 approximately as 1/(2k∆t) times the argument of the complex pulse-pair signal. To do so, one must however consider a bit carefully the ambiguity that is inherent in phase measurements. 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 arg P P ∆t (t 1 , r ) .
At this stage, even a coarse approximation of V NG can be used, as long as it is sufficient to resolve the phase ambiguity.

15
This is important in particular for the onboard processors of satellite instruments, which have to rely on limited quality position/velocity/pointing information and typically can not use the σ 0 distribution information ground segment processors can retrieve from the signal. Care must however be taken to take account of the correction applied by the onboard processor in later processing stages.
A3 Non-geophysical Doppler V NG

A3.1 Overview
This section will describe the different contributions to Non-geophysical Doppler (V NG ). At first order the signature is dominated by the platform velocity and its pointing knowledge. A second point is the determination of the range knowledge which is directly linked to the timing and altitude accuracy. This point is critical for pulse-limited instruments, as KuROS and SKaR,5 for which the pointing is determined for each range bin as function of the altitude. Last, asymmetric variation in NRCS within the antenna azimuth for a given range bin generate a bias in Doppler. All these three elements need to be accurately corrected in order to have an accurate geophysical Doppler measurement.

A3.2 Pointing knowledge
At this point of the discussion, we will leave the general setting and will start orienting our convention choices towards the 10 description of the airborne case of the Drift4SKIM campaign. Problems specific to the space-borne case, such as the need to work in global-scope reference frames or to account for Earth rotation and sphericity, or structural details such as the parabolic reflector featured by the SKaR instrument, do not change the deep nature of the issues, but do introduce a heavy notational burden which tends to obscure the discussion.
From now on, we will work in the simplified setting of the flat-Earth approximation, in which the elevation and incidence 15 angles γ and θ are equal. We will use a platform-fixed reference frame, the origin of which is located at the antenna phase center of the instrument, with x-vector pointing to the geometric front of the platform, y-vector pointing to starboard, and z-vector pointing to the floor, and a local geographic North/East/Down 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 20 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 (N, E, D) 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 level flight towards the North. In the above expression we have used the transparent notation c p → cos(pitch), s r → sin(roll), c h → 25 cos(heading), etc... Other quantities worth introducing are the course c and glide angle g such that the plane velocity vector in the NED frame is In the NED frame, the pointing vector e can be expressed as e = sin(θ) [cos(ϕ)N + sin(ϕ)E] + cos(θ)D.
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.

5
With these notations, and using eq. (A17), one can express V N G as Level flight corresponds to g 0. We thus concentrate on the impact of errors in the first term of the RHS of this equation.

10
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/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 , x), one sees that at 12°incidence, and for a platform velocity of 7000 m/s (space-borne instrument), the SKIM 1 cm/s error budget on horizontal velocity measurements 15 translates to pointing accuracies of 0.3 and 1.4 microradians in incidence and azimuth, respectively (these figures are obtained by allowing each of the pointing errors to consume the full error allocation. As the two terms reach their maximal amplitudes in different parts of the swath, this is not unreasonable). In the airborne case at 120 m/s platform velocity and 3000 m altitude, the corresponding numbers are 18 and 85 microradians for incidence and azimuth pointing accuracy for KuROS, respectively.
In the cross-track viewing geometry of KaRADOC, only the comparatively mild (but still quite demanding) 85 microradians 20 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 have different origins: -the uncertainty in azimuth can be due to an imperfect knowledge of the weighting corresponding to the W (t 1 , r , x) σ 0 (t 1 , x) term in eq. (A22). This can of course come from imperfect platform attitude or antenna orientation information, but also from 25 an imperfect characterization of the antenna radiation diagram or of the distribution of σ 0 on the sea surface.
-the uncertainty in incidence is due to an imperfect knowledge of the radial position of the range resolution bins (yellow striping of the footprint if fig. A1A). This can be due to an imperfect timing accuracy, or to an imperfect knowledge of the vertical separation between instrument and sea surface.
A3.3 Timing and altitude accuracy 30 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 where the incidence angle is θ. In this case θ = arccos(H/r). If now the radar 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 that the surface-tracking algorithm suffers from an error δh, and detects the surface at range H − δr − δh. The data from this range bin will thus be processed using an angle of incidence 5 different from the correct value by Considering δh and δr as independent, we see that at 12°incidence the incidence knowledge requirements expressed above for SKIM and KuROS translate respectively to timing accuracy requirements of 2.4 m and 0.5 m, and to surface-tracking accuracy requirements of 5.4 cm and 1.14 cm.

10
The timing accuracy requirements are easily met in the spaceborne context, but can be challenging in the cost-constrained context of an airborne instrument.
The surface-tracking algorithm, however, does not benefit from the error-compensation that exists for the timing error. The requirement for SKIM is stringent, but the SKIM mission comprises a state-of-the-art Ka-band nadir-looking altimeter capable of reaching this goal. The 1.14 cm altitude tracking requirement is clearly out of reach of the KuROS airborne instrument. Our 15 analysis of its Doppler data will thus be restricted to the side-looking configurations for which, as per eq. (A22), the pointing requirements are much milder.

A3.4 Effective pointing / Azimuth Gradient Doppler
As expressed in eq. (A22), for each range resolution cell V N G 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 20 (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 a 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. (A22) in the limit of a very narrow antenna diagram (which is clearly applicable for 25 SKIM and KaRADOC, less so for KuROS).
Considering first the case of an antenna pointing towards azimuth ϕ A 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 δ(ϕ − ϕ A ). In this limit We recognize in this expression V geo , the estimate of V N G 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 effet 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(ϕ A − c + δϕ) (A28) + sin(g) cos(θ)] .
10 Alternatively, one can choose to express V N G as the sum of V geo , the geometric approximation, plus an Azimuth Gradient Doppler contribution V N G (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(ϕ A − c) + sin(g) cos(θ)] (A32) 15 and V AGD (t 1 , r ) = −V R cos(g) sin(θ) sin(ϕ A − c) (A33) 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 20 than instruments with a broader pattern, 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/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 only through the variations of σ 0 and σ ϕ : and 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, and 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, β = 0. With this choice of coordinates the antenna diagram has separable Gaussian dependencies on α and β. In level flight, when the antenna points towards φ A , sin(α) = sin(θ) sin(ϕ − ϕ A ), tan(β) = tan(θ) cos(ϕ − ϕ A )).

B1 Fixed antenna NRCS correction
The anechoic chamber measurements are very accurate for the antenna alone. However, once integrated into the plane, the antenna pattern is perturbed. This is for instance particularly noticeable in the NRCS measurements in rotating mode, where a spurious azimuthal pattern could clearly be seen, or for fixed-antenna Doppler observations, where a striping pattern is obvious.
We have thus developed a complementary method that relies on the variations of the plane attitude during maneuvers. Using the Figure B1. Reconstructed α, β dependence of the 2-way KuROS antenna diagram. For each 30 s data segment, the level flight values, for which the nadir is at α = 0, β = 0, have been used as reference level to account for geophysical variations in nadir NRCS.
plane IMU, we identify the angular coordinates α and β of the nadir, and use the measured power to map the antenna pattern (using as a reference point the level flight 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 β that 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 α −3dB and β −3dB obtained from anechoic chamber measurements, is: The width parameters in these equations are linked by σ α = α −3dB / 8 log (2), σ β = β −3dB / 8 log (2) The parameter values used in this study are collected in table 1.

10
One cautionary remark is that the illuminated patch at nadir is not infinitely sharp. The measured distribution is thus the convolution of the true antenna diagram by the power distribution at the nadir patch (which depends on the altitude tracking error as well as 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 .
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 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 box for each day.  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.

B3 Fixed antenna Doppler 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 5 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, a RX low-noise amplifier, and a high-gain purpose-built slotted waveguide antenna (shown in Fig. C1). Figure C1. The back (on the left) and the front (on the right) of the antenna The beam can be steered in elevation by changing the instrument working frequency (see Fig. C2A). For the Drift4SKIM experiment, the antenna was mounted in a port-looking configuration, centered on 10°incidence, with a 2°backward-looking 5 tilt to compensate for the aircraft pitch in level flight. Observations at 12°were collected at 33.7 GHz. 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