Sea Surface Ka-Band Doppler Measurements: Analysis and Model Development

: Multi-year ﬁeld measurements of sea surface Ka-band dual-co-polarized (vertical transmit– receive polarization (VV) and horizontal transmit–receive polarization (HH)) radar Doppler characteristics from an oceanographic platform in the Black Sea are presented. The Doppler centroid (DC) estimated using the first moment of 5 min averaged spectrum, corrected for measured sea surface current, ranges between 0 and ≈ 1 m/s for incidence angles increasing from 0 to 70 ◦ . Besides the known wind-to-radar azimuth dependence, the DC can also depend on wind-to-dominant wave direction. For co-aligned wind and waves, a negative crosswind DC residual is found, ≈− 0.1 m/s, at ≈ 20 ◦ incidence angle, becoming negligible at ≈ 60 ◦ , and raising to, ≈ +0.5 m/s, at 70 ◦ . For our observations, with a rather constant dominant wave length, the DC is almost wind independent. Yet, results confirm that, besides surface currents, the DC encodes an expected wave-induced contribution. To help the interpretation, a two-scale model (KaDOP) is proposed to fit the observed DC, based on the radar modulation transfer function (MTF) previously developed for the same data set. Assuming universal spectral shape of energy containing sea surface waves, the wave-induced DC contribution is then expressed as a function of MTF, significant wave height, and wave peak frequency. The resulting KaDOP agrees well with independent DC data, except for swell-dominated cases. The swell impact is estimated using the KaDOP with a modified empirical MTF.


Introduction
Today's ocean radar scatterometry techniques are mainly based on the analysis of the sea surface backscattering intensity, or normalized radar cross-section (NRCS). Controlled by the surface roughness and look geometry, particularly local incidence angle at the surface, the NRCS can serve to retrieve the near-surface wind and/or wave characteristics, e.g., wind scatterometry [1], wave scatterometry [2], and/or be processed to provide very high-resolution wind and wave ocean surface information [3,4]. A novel approach to further extend scatterometer capabilities is to acquire and process the backscattering phase, or Doppler centroid (DC). The DC reflects the time evolution of backscattered fields, and thus encompasses the sea surface kinematics, including more direct information on underlying surface current properties. as a function of look geometry and environmental variables. Using these data, in Section 4, we develop and validate a semi-empirical DC model which is based on a theoretical DC model (DopRIM [10,11,29]) and an empirical modulation transfer function [40]. In Section 5, we compare the semi-empirical model to other published results and discuss its features and applicability. Section 6 summarizes the results.

Materials and Methods
Experiments were conducted from the Black Sea oceanographic platform (Figure 1) operated by the Marine Hydrophysical Institute of the Russian Academy of Science. This is a truss-tower bottom fixed platform that locates 600 m offshore in a 30-m deep water (44 • 23'35" N, 33 • 59'04" E). The consideration is limited to onshore wind directions (easterly to south-westerly, 100 km < fetch < 400 km) with typical peak wavelengths of ∼30-40 m at 10-15 m/s winds. Such wave conditions ensured applicability of the deep water wave approximation. The consideration was also limited to cases without a strong swell. Local currents typically varied from 0 to 0.3-0.5 m/s.

Radar
Experiments were conducted with a Ka-band continuous-wave dual-co-polarized Doppler scatterometer (Figure 1c). It consisted of two co-aligned transmit/receive horn antennae separated by 22 cm. The radar operated at a hybrid polarization mode transmitting a wave of mixed polarization and receiving two orthogonal (vertical and horizontal) backscattering modes. This configuration implied a contamination of the co-polarized signals (vertical transmit-receive polarization (VV) and horizontal transmit-receive polarization (HH)) by cross-polarized signals, which were, however, very weak for typical wind conditions and can be disregarded (see Appendix B in [39] for more details).
Radar installation enabled us to vary the incidence angle, θ, from 0 to 70 • with various radar-to-wind azimuths, but the pure downwind azimuth was rarely used to avoid platform wind shadowing impacts. At low/high incidence angles (θ less/more than 45 • ), the radar was installed on the top/bottom deck at 12/6-m height, respectively.
The two-way beam-width at a half-power level was 5/7 • in the vertical/horizontal plane for VV/HH polarization, respectively. The radar had a real aperture horn antenna. Its ground footprint size and position depended on the incidence angle and radar height only. For example, at low incidence angle observations from the top deck, the footprint size was ≈2 m by 2 m. At large θ = 70 • observations from the bottom beck the footprint size increased up to ≈3 m by 8 m (see Appendix A in [40] for more details).
The radar employs the homodyne method to extract the Doppler shift. Along with received power calibration, a special test was conducted to evaluate the transmitter oscillator instability and check its impact on the sea surface Doppler shift measurements [41]. Detected Doppler signals (I/Q pair per polarization) were digitized by an analog-to-digital converter at 40 kHz/channel sampling frequency. The Fourier transform of the recorded I/Q data was performed over 0.2 s time intervals yielding quasi-instantaneous Doppler spectra of the raw signal. Then, these spectra were averaged over a five minute interval to estimate mean Doppler spectrum for each five minute interval.

Hydro-Meteorological Measurements
Supplementary instruments included the standard meteorological sensors installed at 23 m height mast ( Figure 1e): wind anemometer, atmosphere temperature, humidity, and pressure. The water temperature sensor was submerged at 3 m depth. The data from meteorological sensors were used to estimate wind speed, U, at the standard 10 m height using the COARE 3.0 algorithm [42]. Currents at 10-m depth were measured by a propeller sensor. Surface waves were measured by the six-wire wave gauge antenna (Figure 1b). The wires were fixed at the center and vertices of a pentagon with 25-cm length edges. Directional frequency elevation spectra were estimated using the extended maximum likelihood method implemented in the DIWASP package [43].
Because Doppler scatterometers were designed to retrieve the surface currents, the standard 10 m depth propeller measurements were augmented by surface currents estimated from video records [44]. The video data were recorded by a digital video camera installed atop the radar and directed into its footprint (incidence angle 0-70 • ). The camera viewing angles were 28 • /47 • in the vertical/horizontal planes, the frame rate is 25 Hz, and the image size was 1440 by 1080 pixels. In order to retrieve the sea surface currents, the Particle Image Velocimetry (PIV) technique (bubble tracking), along with the dispersion relation analysis were used. It was found that the PIV-derived currents are = 0.75%U faster than the dispersion-derived currents, which, it turn, were = 1.3%U faster than the 10m-depth propeller currents, where U is the 10 m neutral wind speed.
For radar data analysis, the currents derived from the wave dispersion analysis were used by default. If this estimate was unavailable (e.g., unsuccessful scene illumination), we reconstructed the currents using the vertical shear, , either from PIV or propeller data. The latter data are used if video records are not available.

Results
This section presents measured Ka-band DC for various environmental conditions and look geometry.
Raw radar records, totalling ≈60 h, were partitioned into five minute samples. DC was calculated via the 1st moment of the Doppler spectrum of five minute fragment (positive Doppler frequency corresponds to approaching target). For each fragment, the wind speed and wave spectrum were estimated, while the surface currents were available only for ∼50% of radar samples. The overall statistics for the 2009-2018 data set is shown in Figure 2. The range of moderate to high incidence angles, θ > 20 • , was covered rather uniformly for winds, U < 12-15 m/s, while the range of low incidence angles, θ < 20 • , was covered for only moderate to calm sea states, U < 7 m/s. With respect to the radar-to-wave azimuth, φ wa (zero corresponds to the upwave direction), the data were mostly concentrated around 120 • at small θ < 10 • . At higher θ, the coverage was more uniform. The data set was refined by filtering out offshore wind short-fetch sea states.

Look Geometry Dependence
Generally, the variability of measured DC was primarily determined by the radar azimuth, i.e., approaching/passing waves cause positive/negative Doppler shifts, respectively ( Figure 3). The magnitude of azimuth variation was determined by the incidence angle and the sea state. Upwind and downwind DC values were approximately the same in magnitude at small θ, but become different at larger θ (upwind magnitude is larger), especially at HH polarization. Maximal observed DC was ≈1/1.2 m/s at VV/HH polarization at θ = 70 • in the upwind direction.
In the crosswind direction, the DC was unexpectedly non-zero (green symbols in Figure 3): at small θ < 30-40 • , the crosswind DC was systematically negative V ≈ −0.05-0.1 m/s while, at large θ > 60 • , it approached positive values close to those observed in the near-upwind direction, V ≈ 0.45-0.65 m/s. In contrast with Figure 3, which shows only nearly co-aligned winds and waves (|φ wi − φ wa | < 25 • , where φ wi and φ wa are the radar-to-wind and radar-to-wave azimuth, respectively) to outline the incidence angle behavior of the DC, Figure 4 illustrates azimuth DC dependencies as a function of radar-to-wind direction for all observed wind-wave situations. To first order, the data followed a regular, cosine-like dependency for co-aligned winds and waves (black symbols). For differing wind and wave directions, the DC points stand out of the regular dependency, indicating that the DC is a complex function of both wind and wave directions. Figure 4a,b confirm that crosswind DC is somewhat negative at small incidence angles, θ = 20 • , while it becomes positive (even for co-aligned wind and waves) at θ = 70 • .

Sea State Dependence
Previous studies have shown an apparent DC dependence on the wind speed for C-band measurements at θ = 40 • [14], and much less pronounced for Ka-band at θ = 56 • for wind speeds U = 4-12 m/s [21]. Plotting our measurements against wind velocity projected on the radar look direction indicates a clear dependence on the sea state ( Figure 5). Only VV data are shown, as HH data look similar. The DC sensitivity to the sea surface current is remarkably dependent on the incidence angle (indicated by colors in Figure 5a). Although the DC followed the sea surface current, a positive/negative bias for upwind/downwind direction, respectively, is also present. Its magnitude increased with decreasing incidence angle. This bias is to be attributed to the wave-induced DC contribution. It is demonstrated more clearly by Figure 5b,c, showing DC dependency on wind speed and characteristic magnitude of wave orbital velocity projection on the the radar incidence plane. Note that u = ω p H s was used as a characteristic magnitude of wave orbital velocity, where ω p is the wave peak frequency, and H s is the significant wave height. Hereafter this projection, u cos φ wa , is referred to as "wave range velocity". The scatter in these plots was caused by (i) θ dependence of the wave-induced contribution, (ii) differences between wave and wind azimuths. This is revealed by bottom panels in Figure 5 showing the same data clouds as in the top panels, but color-scaled versus the deviation between wind and wave range velocity components. Like for waves, the wind range velocity is defined as the projection of wind velocity on the radar incidence plane, U cos φ wi . As a measure of the above deviation, the distance (along wave range x-axis) between an individual data point and a linear fit to the whole data cloud was used (Figure 5h). It is clearly seen in Figure 5f,g that strong deviations of the DC from its regular wind/wave behavior are coincident with strong deviations between wind and wave range velocity components. At first glance, the comparison in Figure 5 confirms the presence of sea state dependent DC with magnitudes close to those predicted by CDOP [14]. In fact, the DC was governed by the range projection of wind/wave velocity, ∼cos φ, rather than the velocity magnitude. This is highlighted in Figure 6, showing the DC as a function of wind speed, U, for various azimuths (specified by different colors). In line with the DopplerScatt observations [21], no significant wind dependency can be detected in Figure 6.
The main objective of this study is to present our Ka-band platform observations in an easy usable form. Normally, e.g., for NRCS analysis, this is achieved by fitting observation data as a function of look geometry and environmental variables, also known as the GMF. As discussed in the above data overview, a plain polynomial fit based solely on the look geometry parameters (such as incidence angle, azimuth) cannot be effectively applied because, among others, the fitting function should involve the sea state parameters. In general, all environmental parameters, including wind speed and direction, significant (dominant) wave height, period and direction should be accounted for. But given the data sparsity, a fit to a large number of variables cannot be well constrained. Thus, for practical purposes, a semi-empirical approach is adopted instead of purely empirical GMF. It is based on a functional form of the sea surface radar echo suggested by a theoretical model (DopRIM [10,11,29]) with scaling parameters fitted by our observations. Below, the DopRIM is briefly described. It is further simplified based on observed Doppler signal properties, available wave parameterizations, and environmental condition limitations of our platform measurements. The DC measurements presented above are used for validations of the proposed semi-empirical model.

Background (DopRIM Approach)
Our approach follows the Doppler radar imaging model (DopRIM) [10,11,29]. It accounts for small-scale waves corresponding to individual scatterers, modulated by large-scale tilting waves. The sea surface DC is an equivalent current, V = π f D /k r , corresponding to the Doppler frequency shift, f D , averaged over large temporal and/or spatial scales, where v dr is the drift speed, φ dr is the radar-to-current azimuth angle, v sc is the inherent scatter velocity projection on the radar incidence plane, u is the line-of-sight projection of the wave orbital velocity, σ = σ + σ is the sea surface NRCS, represented as a sum of the mean level and variation, θ is the incidence angle, k r is the radar wavenumber. The first term in Equation (1) describes the sea surface current contribution, and the second one is the mean scatterer velocity relative to the sea surface. The third term describes contribution arising from the correlation between large-scale surface (orbital) velocities and (tilt-and hydro-) modulations of the scatterer NRCS.
In the DopRIM, the DC contribution includes resonant and non-resonant scattering mechanisms, and the total DC reads: where n denotes the scattering mechanism (n = 1-Bragg, n = 2-quasi-specular reflection from regular non-breaking surface, n = 3-scattering from breaking waves), P n = σ n / ∑ σ is the partial NRCS contribution by the corresponding mechanism. At moderate and large incidence angles, θ > 20-30 • , the Bragg-resonant and wave breaking mechanisms both contribute, while, at small θ < 20-30 • , quasi-specular reflections dominate. Below we describe all DopRIM contributions separately and consider how they can be simplified.

Moderate and Large Incidence Angles
The inherent velocity corresponding to the resonant scattering mechanism corresponds to the Bragg wave phase velocity, where g is the gravity acceleration, γ w is the surface tension, and k br = 2k r sin θ is the Bragg wavenumber. The backscattered intensity depends on the Bragg wave spectrum, thus the Bragg scatter velocity also depends on the angular distribution of wave spectrum in the direction towards and away from a radar, where B is the directional curvature spectrum, and φ wi is the radar-to-wind azimuth.
In the DopRIM, radar scattering from rough surface patches associated with breaking of waves with wavenumbers k < min(k wb , k r /10), where k wb = 2π/0.3 rad/m, significantly contribute to the surface NRCS. Doppler velocity of breaker-facets is equal to the mean phase velocity, c, of breaking waves weighted with their contribution to the fraction of the sea surface covered by breaking zones: where Λ(k) ∝ (u * /c) 2 k −1 B(k) is the Phillips Λ-distribution [45] for the total length of breaking crests per unit area associated with wavenumber, k, and wave phase velocity, c, u * is the wind friction velocity. If B = const then, according to Equation (5), the inherent velocity of breaking facets is c wb = 2c(k wb ), or about c wb ≈ 1.4 m/s. However, as found in [40], the observed sea surface Doppler shift does not reveal strong spikes at θ < 65 • . Also, wave spectra derived from Doppler velocity agree well with in situ wave gauge spectra. It suggests that, although breakers significantly contribute to spike-like NRCS, their impact on instantaneous Doppler shift is weak. As suggested in [40], enhanced surface roughness patches on crests of breaking waves are "embedded" into the water, and thus their Doppler velocity corresponds to the orbital velocity of breaking waves. Joint analysis of collocated radar/video measurements of Doppler/optical velocity of breakers [46] further revealed that Doppler velocity of breakers is slower than the advance velocity of breakers by a factor of ≈4 (traced by whitecap movement). If one assumes that the 1/4 factor corresponds to the characteristic steepness of breaking waves, then the measured Doppler velocity of breakers corresponds to their orbital velocity. Following, one may conclude that the inherent Doppler velocity of breakers defined by Equation (5) is to be scaled by a factor of 1/4. This gives c wb = 2c(k wb )/4 or c wb ≈ 0.35 m/s, that is comparable with Bragg phase velocity.
Notice that this c wb estimate is valid for moderate incidence angles only, θ < 65 • . At larger incidence angles, especially at HH polarization, the Doppler shift demonstrates quite different features including strong spikes correlated with breaking wave events [40]. However, this study is limited to moderate incidence angles, θ < 65 • , and specific Doppler shift features observed at larger incidence angles are not considered.
The last term in Equation (1) is caused by correlated variations of facet NRCS, σ , with its long wave orbital velocity, u. This so-called wave-induced Doppler velocity can be expressed in terms of a complex modulation transfer function (MTF) [47,48], where Ω is the long wave angular frequency, C is the long wave phase velocity, G is the geometrical coefficient projecting wave orbital velocity onto the radar line-of-sight, M is the complex MTF consisting of tilt-and hydro-MTF, and φ wa is the radar-to-wave azimuth angle.
In the DopRIM, the tilt and hydro-MTF are evaluated separately for Bragg and wave breaking scattering mechanisms (see details in [29,49,50]). Here, we rely on an empirical MTF developed in [40] using the same data set for modulating wave frequency from 0.2 Hz to 0.8 Hz. In contrast to the DopRIM, the empirical MTF combines all mechanisms and there is no need to separate their contributions. Thus, Equation (6) with empirical MTF, M, accounts for both tilt-and hydro-wave-induced contributions of Bragg and wave breaking scattering mechanisms at moderate incidence angles, θ > 20-30 • .

Small Incidence Angles
At small incidence angles, θ < 20-30 • , quasi-specular reflections provide dominant contribution to the surface NRCS. The mean Doppler velocity of specular points is derived by Longuet-Higgings [51], and its spectral form is given in [24,29], where φ is the relative to wind azimuth of the spectral component associated with wave vector, dk is the mean square slope (MSS) components of the modulating waves in the upwind and crosswind direction, respectively, and integral in Equation (8) is taken over the large-scale surface waves, i.e., over spectral domain satisfying the condition: k < k r /4. Equation (8) represents the inherent mean velocity of specular points. However, given a small radar footprint, apparent modulations of NRCS and Doppler velocity by dominant surface waves were observed [40] at low incidence angles. These observations presume that the two-scale concept, with differing scales of waves providing radar backscattered intensity and modulations, can be adopted. This can also be corroborated by the evaluation of integrals in Equation (8). For a characteristic curvature spectrum, B(k) ≈ const, the MSS is dominated by the high frequency tail of spectrum, while the integral involving the phase velocity is mainly governed by the low frequency part of the spectrum (dominant waves). Following this concept and using the known NRCS approximation for low incidence angles [24,52,53], along with the definition of tilt-MTF, Equation (8) is rewritten for a single dominant long wave as, which takes a form similar to Equation (6) describing the wave-induced Doppler velocity. Such similarity is not that surprising. Equation (8) describes the mean velocity of specular points, resulting from the linear superposition of random waves of different scales [51]. The majority of specular points is linked to large local slopes of short scale waves, which are advected by fast, but small slope, dominant surface waves. As such, the mean specular point velocity is linked to parameters describing these dominant waves, as predicted by Equation (11).

The Semi-Empirical Model
The similarity between Doppler velocity relationships at low, Equation (11), and moderate, Equation (6), incidence angles helps suggest a more unified relationship, valid for all incidence angles. This is especially tempting because the empirical MTF developed in [40] applies to all incidence angles and accounts for all mechanisms responsible for radar backscattering, including quasi-specular reflections at low incidence angles, and Bragg/non-Bragg scattering at larger incidence angles.
First, the inherent Doppler velocity of Bragg waves and breakers can be merged. As discussed above, the inherent velocity of breakers (about 0.35 m/s) is close to the Bragg wave phase velocity. At moderate incidence angles, the sum of weights for Bragg and breaking mechanisms is close to 1. For the sake of simplicity, the Bragg phase velocity is used as a measure of the inherent scatterer velocity, Equation (4), which is evaluated using the angular distribution of Bragg wave spectrum reported in [54] and represented here as, The Bragg scatter velocity, Equation (4), is not valid at small θ where quasi-specular reflection dominates. However, since line-of-sight projection of Bragg velocity tends to zero at θ → 0 due to sin θ-factor, we assume that Equations (3), (4), and (12) determine the inherent scatterer velocity in the whole range of θ in our model.
Combining Equations (11) and (6), the unified expression for semi-empirical DC model (KaDOP) valid from small to moderate incidence angles, θ < 65 • , reads: or in terms of directional frequency spectrum, S(ω, φ), The frequency integral in Equation (14) is defined by the low-frequency part of the spectrum and thus not sensitive to the upper limit of integration. The empirical MTF developed in [40] is an essential component of this model. Notice that this MTF was originally developed for wave frequencies 0.2-0.8 Hz, and strictly speaking, does not apply to lower frequencies, <0.2 Hz. However, in the absence of a better alternative, we will use it for DC estimates even for wave frequencies out of the MTF validity range.
The third spectrum moment can be evaluated as: For practical applications, when several, N, wave systems are present (mixed sea), Equation (14) can be reduced to where N is the number of wave systems, v dr , φ dr are the sea surface current speed and direction, respectively, v sc is the inherent scatterer velocity determined by Equations (3), (4) and (12), φ waN is the radar-to-Nth-wave-system direction, U is the 10 m-wind speed, H sN and ω pN are the significant wave height (SWH) and peak frequency of the N-th wave system, G is the geometrical coefficient determined by Equation (7), M(θ, φ waN , U) is the empirical MTF [40] (see also Appendix A). The coefficient, β N , in Equation (16) depends on the shape of the frequency spectrum and thus differs for broad wind wave spectrum and narrow swell spectrum. As a first guess, it is estimated from the Pierson-Moskowitz spectrum [55], for which β ws ≈ 0.2. For very narrow (delta-function type) swell spectrum, β sw = 1/16.

Model Validation
To validate the KaDOP model, Equation (14), it was first applied to the data subset with known wind and wave conditions, but unknown surface currents (Figure 7, left three columns). As a proxy for 10-m current shear, 1.5%U was used that was half of the "classical" 3%U but close to 1.3%U dependence observed at the platform site [56]. In order to emphasize the importance of angular spreading of wave spectrum on DC, we performed simulations with both, 2D and 1D wave spectra obtained from wave gauge measurements. For the latter case, the directional spectrum S(ω, φ) was replaced by the omni-directional one, the MTF was specified for the mean wave-direction, and the integration over azimuth was omitted. Further, simulations with 2D/1D spectrum are referred to as KaDOP-2D/KaDOP-1D, respectively.
The KaDOP-2D worked better than the KaDOP-1D which systematically overestimated the measured DC. This is explained by unrealistically narrow azimuth distribution with all the wave energy coming in the same direction, φ wa . In the 2D case, the real azimuthal distribution was accounted for, to enable a better DC simulation.
At small incidence angles, the range DC as well as its errors are "amplified" by the 1/ sin θ factor that leads to larger data scatter. At moderate to large θ, the empirical model works better than at small θ except at θ >65-70 • where the model exhibits significant variability while the measurements do not. This scatter is more pronounced at HH polarization. This is explained by the breaking wave contribution that is indirectly accounted for through the linear empirical MTF. The increased model scatter at large θ suggests that the concept of linear MTF is applicable only for θ < 65-70 • . At larger θ with large NRCS modulation magnitudes, the linear MTF approximation is no longer accurate [41]. However, we leave this issue out of the scope of this study by limiting the applicability of the KaDOP to small and moderate incidence angles, θ < 65 • .
The subset with known drift currents (Figure 7, right three columns) demonstrates similar behavior but with a smaller root-mean-square error. Note that the model is sensitive to the angular width of the wave spectrum (compare Figure 7d,j,e,k light blue points). But, 2D wave spectra were evaluated indirectly from wave gauge records. Data shown above are based on the maximum likelihood method. Application of an alternative maximum entropy method yielded similar spectra but with a narrower angular width of spectrum peak. Thus, uncertainties in wave spectrum estimation also contributed to the uncertainties of the empirical model.
Overall, the correlation coefficient between the KaDOP and measurements (>0.9) was high at θ < 65-70 • that allows extrapolating the KaDOP on an arbitrary wave spectrum. The KaDOP evaluated using Pierson-Moskowitz spectral shape, Equation (16), are shown in Figure 7c,f,i,l. Comparing these simulations with KaDOP-1D simulations (Figure 7b,e,h,k) one may conclude that they are very similar. This fact indicates that the simplified DC model based on observed SWH and spectral peak frequency had similar skill as more complex model utilizing measured wave spectra. For further discussions, Equation (16) is used with wave input parameters expressed in terms of SWH and ω p .  (14); 1D-sea, Equation (14); equivalent Pierson-Moskowitz (PM) spectral shape, Equation (16). (right three columns) are the same, but for cases with known wind drift. Correlation coefficient, R 2 , and root-mean-square error (RMSE) are shown in each panel.

Discussion
In this section, we compare the KaDOP with available measurements in different bands: • X-band, Wavemill data collected in the Irish Sea [35] (their Table 1 To conduct calculations, the drift current is assumed to vary from 0 to 3% with the mean value of 1.5% of the wind speed. Bragg wave phase speed was set to its characteristic value in the Ka-band, ≈0.35 m/s at moderate incidence angles, which is close to values in the C/X/Ku-bands, ≈0.23-0.26 m/s. As a wave input to the KaDOP, the classical Pierson-Moskowitz relationships for a fully developed sea [55] are used, For simplicity, it is assumed that waves propagate in the wind direction, thus radar-to-wind and radar-to-wave azimuths are the same, φ wi = φ wa = φ.

Look Geometry Dependence
In general, the KaDOP, Equation (16), reproduces fairly well all independent observations, but the CDOP that is somewhat underestimated (Figure 8). The DC at HH polarization is remarkably higher than at VV polarization in line with the dual-polarized SAXON-FPN, IAP and CDOP data. The IAP data are in good agreement with KaDOP except for downwind HH polarization. Note that in the KaDOP, the polarization sensitivity is determined by the MTF polarization dependence only. The AirSWOT is somewhat overestimated, probably due to the fact that a full-developed Pierson-Moskowitz spectrum is not exactly the case for fetch limited AirSWOT waves. The same applies to the Wavemill which is based on complicated mixed sea conditions. The sign of crosswind DC discussed earlier in the raw data analysis (Figure 4) is also reproduced by the KaDOP. The crosswind DC zeros at θ ≈ 60 • in line with the DopplerScatt at θ = 56 • . In contrast with the KaDOP, the CDOP and Wavemill models predict small but positive crosswind DC at moderate θ = 20-40 • . This is also demonstrated by azimuth dependencies of KaDOP DC (Figure 9) that better agree with the DopplerScatt than with the CDOP and Wavemill.

Wind Speed Dependence
KaDOP wind speed dependence of the DC is generally consistent with the DopplerScatt and CDOP, except for the CDOP in the upwind direction ( Figure 10). The difference between DC wind speed dependencies at different incidence angles can be explained by the balance between orbital velocity and MTF magnitudes that differently respond to the wind speed. For an omni-directional elevation ω −5 -spectrum, the magnitude of orbital velocity, u = aω, is inversely proportional to wave frequency, u ∼ ω −1.5 . For a given wave age, U/c p = const, the magnitude of orbital velocity is approximately a linear function of wind speed, u ∼ U 1.5 . On the other hand, the MTF increases with wind speed decreasing. To the first order, the tilt-MTF is inversely proportional to the MSS, M T = ∂ ln σ/∂θ ∼ 1/ζ 2 , that in turn linearly increases with U. But, the contribution of hydrodynamic modulations makes the total MTF a more complex function of wind speed. Thus, the DC wind speed dependence is determined by relative changes of magnitudes of orbital velocity and MTF. In addition, it also depends on the wind drift velocity, which is also a linear function of U sin θ.

Upwind/Downwind Asymmetry
There exists a slight but important difference between DC magnitudes in the upwind and downwind directions, an upwind/downwind asymmetry (UDA), UDA = |V up |/|V down |. It is attributed to the hydrodynamics modulation that impacts the DC in opposite ways for upwind and downwind look directions. This effect is expected to decrease with wind speed following the related decrease in the magnitude of hydro-MTF. There are noticeable differences in UDA behavior for different empirical models (Figure 11). In the CDOP, the difference between upwind and downwind DC is always positive and the UDA is about 0.2-0.4 higher at VV than at HH polarization. In contrast, the KaDOP does not show such strong UDA polarization difference. At θ = 56 • , KaDOP UDA decreases with increasing wind speed from ≈ 1.3-1.4 at U = 5 m/s to ≈ 0.8-0.9 at U = 15 m/s in line with the DopplerScatt. The latter UDA < 1 values indicate that the magnitude of downwind DC exceeds the magnitude of upwind DC. Such behavior is explained by the impact of hydrodynamics modulation that changes its phase depending on wind speed (see e.g., Figure 11d in [40]). At moderate winds, U = 7 m/s, the peak of hydro-MTF locates on crests/forward slopes of long modulating waves but shifts on rear slopes at stronger winds, U > 10 m/s. Figure 11. Upwind-to-downwind Doppler centroid ratio versus incidence angle for various wind speeds, U.

Crosswind Doppler Centroid
One of the most challenging KaDOP features is its non-zero crosswind behavior. Crosswind KaDOP DC is remarkably negative for θ < 60 • . In contrast, CDOP and Wavemill crosswind DC is about zero in this range of θ. This discrepancy may be interpreted in the following way.
Crosswind DC (its sign and magnitude) is controlled by phase and magnitude of the hydro-MTF (impact of crosswind tilt-MTF is zero). The empirical MTF, which is the component of the suggested KaDOP model, predicts the shift of the backscattering modulation peak on the windward (rear) slope of modulating waves, that leads (as follows from Equation (16)) to negative crosswind DC. As discussed in [40,57], the observed shift of modulation peak on the windward wave slope results from wave-induced airflow acceleration over windward wave slopes. This effect occurs only for developing wind waves for which the shift of maximum wind velocity on the rear slope is directly related to the appearance of positive air pressure anomaly over the leeward wave slope, the pressure anomaly that provides energy transfer from wind to waves. Thus, the shift of small-scale surface roughness onto the windward wave slope is a very natural phenomenon, which is, nevertheless, present if (and only if) these dominant wind waves are rather young. This is typical for the Black Sea platform observations characterized by fetch-limited conditions. A similar Ka-band MTF phase behavior was found in the SAXON-FPN measurements at θ = 45 • [58] (see their Figure 4) conducted from an offshore platform in the North Sea.
When wind waves are fully developed, wave-induced wind undulations almost vanish over dominant waves (at least, because the wave phase velocity is close to wind speed). In this case, there is no mechanism which could shift the backscattering modulation peak onto the windward slope. As a consequence, the crosswind DC should be about zero. In case of a swell, swell-induced wind acceleration occurs either over wave crests or in wave troughs depending on the relative wind-swell propagation direction [59,60]. Thus, a swell does not produce a shift of aerodynamic modulations on either of swell wave slopes. In case of a swell, we should not expect the appearance of slope correlated component of hydro-MTF which could contribute to the crosswind DC.
Overall, the observed negative crosswind DC can be most plausibly attributed to the young wind sea effect. This type of waves is typical for coastal platform measurements and, therefore, is a specific feature of our empirical MTF, M ws (with "ws" standing for wind sea). This empirical MTF predicts a shift of hydro-MTF to the rear slope of the modulating wave. But, it should not probably be extrapolated to open ocean conditions where fully developed sea and swell are much more typical. Such specific limitation of the M ws can easily be fixed (to the first order) by its modification via forcing the crosswind MTF phase to zero (see Appendix A for details). We suggest using this modified MTF, M sw (with "sw" standing for swell), for applications of the KaDOP model in open ocean conditions.

Swell Impact
In the open ocean, swell waves may have extremely long wavelength and phase speed, C sw ≈ 30 m/s. Having only very small steepness, (AK) sw = 0.01-0.05, these waves, nevertheless, induce orbital motions with magnitudes, C sw (AK) sw = 0.3-1.5 m/s, which may be of the order of wind-wave-induced DC. Swell may travel thousands of kilometers from its generation source and thus have arbitrary direction relative to the local wind wave direction. The KaDOP, Equation (16) with M ws and M sw used for wind sea and swell contributions, allows estimating the possible influence of ocean swell on the DC for various swell-to-wind directions.
As a first guess for possible swell impacts, the mixed sea DC is simulated in Figure 12 using wind sea M ws at U = 6 m/s and the Pierson-Moskowitz spectrum. The impact of the swell is added using M sw for various swell-to-wind azimuths. Parameters of the swell are somewhat exaggerated: For co-aligned wind and swell, the impact of swell increases the magnitude of upwind and downwind DC by ≈ 100% while the crosswind DC remains unchanged, as expected (Figure 12a). The same holds for the case of upwind swell (Figure 12c), but wind sea and swell contributions compensate each other in this case. Notice, this is only a qualitative simulation, in part, because M sw phase is set to zero in the crosswind direction. In line with discussions in the previous section, it could be noted that more realistic MTF should also include the impact of wave-induced wind stress modulation and the corresponding crest/trough-correlated hydro-MTF which may impact upwave and downwave (relative to the swell) DC estimates.

Conclusions
This paper presents measurements of the Doppler spectrum centroid (DC) conducted from the Black Sea research platform in an attempt to systematize the multi-year in situ co-polarized data in a wide range of incidence angles, azimuths, and winds obtained in well controlled sea state (wave) conditions, and provide a model to describe the Doppler characteristics. In contrast with the NRCS, the observed DC is primary dependent on the characteristics of energy containing waves. The DC cannot be fully described in terms of just wind speed like it is usually done for radar cross-section geophysical model functions. It is preferable to use the observed wave spectrum, or its parameterization, rather than the wind-alone as an input parameter to the DC model.
In this work we use a simplified semi-empirical model to fit Ka-band DC data (KaDOP, source codes are available in the Supplementary Materials). It is based on the Doppler Radar Imaging Model (DopRIM, [11,29]) involving two-scale surface decomposition. Further simplifications include the usage of an empirical MTF [40], instead of the theoretical MTF. Empirical MTF naturally includes all mechanisms (tilt, hydro, Bragg, wave breaking) and thus can be used to describe wave-induced DC contribution at moderate incidence angles, θ. At small θ, the mean specular point velocity is determined by the same equation as the wave-induced contribution at moderate θ. Thus, these two terms are combined and the empirical MTF is solely used to model the DC for small to moderate incidence angles. The average wave breaking Doppler velocity, originally considered significant in the DopRIM, is shown to be close to Bragg wave phase velocity due to the effect of embedding (slowing down) of wave breaking scatterers into the water [40,46].
The KaDOP adequately reproduces the measured DC when in situ wave gauge spectra and surface currents are used as model inputs. The only exception is the large incidence angle case, θ = 70 • , for which HH polarization model DC is inadequate and clearly indicates missing effects of wave breaking and resulting MTF non-linearity. Thus the KaDOP applicability is limited to θ ≤ 65 • . This limitation allows to keep the issue of wave breaking impact on the DC at large incidence angles out of the scope of the present study.
Based on fairly good KaDOP performance versus the ground truth training data set, we check the effect of substitution of observed wave spectrum by empirical spectrum models. In this approximation, the KaDOP is a function of the peak wave frequency and SWH. With Pierson-Moskowitz parameters [55] as an input, the KaDOP agrees well with the Ka-band DopplerScatt [21] and AirSWOT [19] data, as well as with the Ku-band SAXON-FPN data [32], X-band IAP [33] (except for HH downwind). The C-band CDOP [14] and X-band Wavemill [35] data are systematically underestimated in the upwind direction.
The KaDOP suggests negative DC values in the crosswind look direction (equivalent motion away from an observer) at θ < 50 • . Qualitatively, this is explained by the hydrodynamics MTF that peaks on windward wave slopes. Speculatively, this effect reflects an acceleration of air flow above windward slopes of developing waves. Indirectly, this feature is confirmed by independent Ka-band MTF measurements from the SAXON-FPN experiment [58]. These latter data also reveal a negative MTF phase in near-to-crosswind look directions. However, the airflow acceleration is possible only for wind faster than wave (developing wind sea), while for swell-dominated sea, the hydro-MTF is not slope-correlated. As the first guess for the swell-MTF, unavailable from our data, we refit the original wind-sea MTF with crosswind phase forced to zero.
With both wind-sea and swell-MTF, the KaDOP shows a clear impact of swell on the DC measurements. Thus, the disagreement between KaDOP and mixed-sea observations (CDOP and Wavemill), besides differences in microwave bands and corresponding MTF, can be attributed to peculiarities of the real spectrum that can differ from the Pierson-Moskowitz model and can be impacted by a swell.
Overall, the main purpose of this study is to show that proposed simplified semi-empirical approach based on observed MTF and integral parameters of long wave spectrum, is promising. With empirical MTF and a simple Pierson-Moskowitz spectrum, the KaDOP is capable of reproducing most of the available independent data sets. The KaDOP is validated against our Black sea fetch-limited measurements (fetch < 400 km, wind speed < 15 m/s, dominant wavelength < 40 m, SWH < 1.5 m). Further improvements of DC predictions are tied with better wave spectrum description and MTF parameterization, no matter empirical or theoretical.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Empirical MTF Modification
The empirical MTF fit [40] used in this study is a polynomial function of θ, φ wa , and U, with magnitude (expressed as the natural logarithm), ln |M|, and phase, M/|M|, fitted independently: Originally, the MTF fit was not constrained at θ = 0 • that resulted in azimuth oscillations of MTF and wave-induced DC. Obviously, any variable measured by a "perfect radar" in the nadir direction should be azimuth independent. At θ = 0 • , the phase of MTF is close to 180 • indicating that peak of backscattering intensity locates in wave troughs which are smoother than wave crests and thus produce stronger specular backscattering. To adjust the MTF at θ = 0 • , it was simulated from Equation (A1) on a uniform {θ, φ wa , U}-grid. Then, the MTF phase was set to 180 • at θ = 0 • . After the above adjustment, it was fitted again using the same polynomial function.
The crosswind behavior of this MTF is explained by the effect of aerodynamic hydro-MTF attributed to developing wind sea conditions typical of our observations. Thus we recall it as wind-sea-MTF, M ws . New coefficients for M ws , differ from the original ones [40] only in a "cosmetic" way due to the nadir phase adjustment. They are presented in Table A1.
For a swell-dominated sea, the crosswind MTF phase should be set zero that cannot be directly obtained from our developing wind sea observations, but can be easily fixed by forcing the crosswind phase of M ws to zero in a way similar to that described above for the nadir correction. The wind-sea-MTF is simulated on a uniform {θ, U}-grid, but now only for upwind, crosswind, and downwind directions. Then, the phase of crosswind samples is set to zero, and the MTF is refitted. Coefficients for this alternative fit, swell-MTF, M sw , are presented in Table A2.