Joint retrieval of ocean surface wind and current vectors from satellite SAR data using a Bayesian inversion method

This paper presents a method for joint retrieval of the ocean surface wind and current vectors using the backscatter and the Doppler frequency shift measured by spaceborne single-beam single-polarization synthetic aperture radar (SAR). The retrieval method is based on the Bayesian approach with the a priori information provided by atmospheric and oceanic models for surface wind and currents, respectively. The backscatter and Doppler frequency shift are estimated from the along-track interferometric SAR system TanDEM-X data. The retrieval results are compared against in-situ measurements along the Swedish west coast. It is found that the wind retrieval reduces the atmospheric model bias compared to in-situ measurements by about 1 m/s for wind speed, while the bias reduction in the wind direction is minor as the wind direction provided by the model was accurate in the studied cases. The ocean model bias compared to in-situ measurements is reduced by about 0.04 m/s and 12∘ for current speed and direction, respectively. It is shown that blending SAR data with model data is particularly useful in complex situations such as atmospheric and oceanic fronts. This is demonstrated through two case studies in the Skagerrak Sea along the Swedish west coast. It is shown that the retrieval successfully introduces small scale circulation features detected by SAR that are unresolved by the models and preserves the large scale circulation imposed by the models.


Introduction
Ocean surface winds and currents play an important role in atmosphere-ocean dynamics. Wind stress is the primary force for ocean surface circulation, thus an important input to ocean numerical models. Ocean models are usually forced by atmospheric reanalyses, which contain errors (Belmonte Rivas and  and thus induce uncertainties in ocean models. In addition, coupled atmosphere-ocean modelling can be improved by better exploitation of satellite radar stress measurements (Trindade et al., 2020). Ocean currents transport heat, salinity and nutrients and thus affect the weather and the marine life. There is, for instance, an increased interest in the effect of surface currents on the estimation of the air-sea fluxes (Dawe and Thompson, 2006). Surface currents are also highly relevant for applications involving prediction of the drift of floating objects, e.g. plastic, (van Sebille et al., 2020). Given the largeness of the ocean, the observations obtained from moored buoys and ships are sparse and drogued drifters measure currents at 15 m depth. Remote sensing technologies and particularly spaceborne sensors represent complementary means for obtaining global measurements of surface winds and currents. While ocean surface wind estimation from spaceborne sensors, particularly scatterometers, can be considered a mature technology, ocean current estimation from space is still in its early phase. Simultaneous estimation of winds and currents are highly desirable for understanding the air-sea coupling .
There are several methods for estimating surface winds and currents from spaceborne sensors data. Here, we focus only on the active microwave sensors (radar). The widely used satellite ocean wind products are derived from the so-called virtual constellation of scatterometers , which consists of an international collaboration to build a global constellation of scatterometers and standard wind products. The second source of ocean wind products from satellites is altimetry (Witter and Chelton, 1991). Scatterometers have the advantage of large spatial coverage and azimuth angle diversity hence the ability to retrieve the wind direction. However, the scatterometerderived wind direction usually requires ambiguity removal (Hoffman et al., 2003;Vogelzang et al., 2009). Altimeters are furthermore commonly used to estimate geostrophic currents from sea level measurements (Isern-Fontanet et al., 2017). Recently, (Rodriguez et al., 2018) proposed a concept design based on a pencil-beam rotating Doppler scatterometer for estimating both ocean vector winds and currents. Another mission called Surface Water and Ocean Topography (SWOT) is proposed to address the shortcomings of conventional altimetry by providing a wide swath and high-resolution mapping of the water elevation (Rodriguez, 2015). The disadvantage of both scatterometers and altimeters is the coarse spatial resolution (10-50 km), this might be suitable for synoptic and mesoscale studies but barely resolve eddies and baroclinic coastal currents. These are particularly relevant in the shelf seas and coastal areas, where the topography, river runoffs, etc. contribute to the complexity of the circulation. In these areas, scatterometers and altimeters are limited by land contamination.
Synthetic aperture radar (SAR) is a high resolution imaging radar. SAR and scatterometers measure the power backscattered from the surface of the Earth. Over ocean, the backscattered power is a measure of the roughness of the sea surface. The backscattered power is sensitive to both wind speed and direction (Horstmann and Koch, 2008). In fact, the backscattered power relates to the wind stress, which in turn relates to the wind vector . Thus, the wind stress is usually expressed as the stress-equivalent wind (de Kloe et al., 2017). The received power is converted to the normalized radar cross section (NRCS), also called the backscatter, backscattering coefficient or sigma naught (σ 0 ). This conversion, ideally, eliminates the instrument dependent parameters (e.g. antenna gain) to obtain a quantity which depends on the dielectric and geometric properties of the surface, frequency, polarization and incidence angle. SAR also measures coherently the phase of the return signal from which the Doppler shift can be estimated. The Doppler shift can be estimated using a single SAR Doppler centroid analysis (DCA) technique (Hansen et al., 2011) or along-track interferometric SAR (ATI-SAR) (Romeiser et al., 2010a). The latter requires two SAR images of the same area acquired within a short time delay. This Doppler shift is directly related to the relative motion between the satellite and the ocean surface (Ardhuin et al., 2019). After removal of the non-geophysical contribution due to the relative velocity of the satellite and the solid Earth, the residual (geophysical) Doppler shift is directly related to the ocean surface velocity. The main contribution to this velocity is due to the orbital velocities of the sea surface waves, hence the relation of the Doppler shift to the surface wind  and due to surface currents (Rouault et al., 2010).
Due to the complexity of the relation between the sea surface parameters and the backscatter, to date, the common approach for wind inversion from spaceborne radar measurement (scatterometer, altimeter and SAR) is based on the so-called empirical geophysical model functions (GMF). These GMFs are usually built by matching the measured backscatter with surface winds provided by reanalysis or in-situ measurements. They are frequency, polarization and incidence angle dependent. The advantage of using GMFs is efficient computation and limited parameterization (incidence and azimuth angles, wind speed and direction). These GMFs are applicable to open slick-and-ice-free water, which assumes that presence of sea ice, oil spill, etc. is well detected and filtered out. However, other effects such as atmospheric stability and air mass density (de Kloe et al., 2017), rain (Xu et al., 2020) and temperature (Wang et al., 2017) influence the relationship between the measured backscatter and the wind. More relevant for SAR, effects such as bathymetry, swell, wave-current interaction, ship wakes and wind mills also influence the relationship between the backscatter and wind.
The wind retrieval from SAR is an inverse problem. The available GMF is a forward model, i.e. it simulates the backscatter from the state variable (wind vector). The inversion can be formulated in terms of retrieving the wind speed and direction or the zonal u and meridional v components of the wind vector. For each pixel in the image, a singlebeam single-polarization SAR provides one backscatter value. Thus, the wind vector retrieval problem from SAR, based on backscatter measurement only, is underdetermined (Portabella et al., 2002). The problem does not become determined even by using the Doppler shift in addition to backscatter  as shown later. Even multibeam sensors such as scatterometers produce ambiguous wind solutions (Hoffman et al., 2003;Vogelzang et al., 2009). Wind retrieval methods can be classified into non-Bayesian and Bayesian.
In non-Bayesian wind retrieval algorithms (called here direct method) the objective is to find an estimate of the wind vector closest to the radar measurements. In this method, the wind direction obtained from an external source is directly inserted into the forward model to extract SAR wind speed. This makes the problem determined, and can be solved for example by least squares method. The wind direction is usually provided by a numerical weather prediction (NWP) model or reanalysis. This method has two disadvantages. It implicitly assumes an error-free direction information. The second disadvantage is the coarse resolution of the NWP wind direction (≈5 − 10 km). It was shown (Fetterer et al., 1998;Horstmann and Koch, 2005) that the wind direction can be estimated from linear streaks in the SAR image and used in the retrieval of wind speed instead of obtaining the wind direction from an external source. This provides a relatively higher resolution wind direction (~1 − 2 km) but suffers from 180 o ambiguity that needs to be resolved. This ambiguity is usually removed using an atmospheric model. However, in many cases the absence of wind features in the SAR image hinders the wind direction extraction.
The Bayesian method has been first proposed and applied to scatterometer-derived winds (Hoffman et al., 2003;Cornford et al., 2004;Stoffelen and Portabella, 2006;Vogelzang et al., 2009). SAR wind retrieval using variational methods has been first proposed by (Portabella et al., 2002) and then investigated by a few other authors (Danielson et al., 2008;Choisnard and Laroche, 2008;Jiang et al., 2017). In this approach (Portabella et al., 2002), the backscatter measurements are used with a priori information about the wind speed and direction, typically obtained from NWP or reanalysis. In the following, the a priori information is referred to as background. It was found that the retrieval relies heavily on the background wind and particularly the wind direction (Portabella et al., 2002). A similar conclusion was drawn by (Choisnard and Laroche, 2008). The Bayesian inversion, applied to Envisat/ASAR data, was assessed in (Adamo et al., 2014). The author also concluded that the performance strongly depends on the quality of the NWP wind direction. In , it was suggested that the use of Doppler shift in addition to the backscatter and the background wind improves the retrieval results compared to both the direct and the Bayesian approach with backscatter only. The main improvement was found in the wind direction when the background was inconsistent with in-situ winds.
For the current retrieval, there is to our knowledge no reported study based on Bayesian approach for estimating ocean surface currents from a single-beam SAR. A few authors, e.g. (Rouault et al., 2010) investigated the direct method, i.e., retrieval of the radial component of the ocean surface current. For instance, a study of the Agulhas current using Envisat/ASAR was reported in (Rouault et al., 2010). It is worth noting that other authors suggested advanced technologies to estimate the total current vector from SAR. For instance, (Frasier and Camps, 2001) proposed a dual-beam interferometric SAR. Similar design was investigated and demonstrated through an airborne campaign in (Martin and Gommenginger, 2017). The authors studied, through numerical simulations, the performance of wind and currents retrieval  based on a hypothetical concept of a dual-beam dual-polarization ATI system.
In this paper, the joint retrieval of the wind and current vectors over the sea surface is demonstrated using the backscatter, Doppler shift and background wind and current fields. The Doppler shift used in this work is derived from the X-band ATI-SAR system (TanDEM-X) data. To our knowledge, this is the first attempt to simultaneously retrieve the total wind and current vectors from a single-beam single-polarization spaceborne ATI-SAR data.

Observables and forward models
The observation vector y and the state vector x are related by the forward model F as follows where ε o = ε mod + ε rep + ε inst is the observation error which is composed of the forward model error, representation error and instrument error respectively and p groups all the model parameters.
In the context of wind and current retrieval from SAR, x is the wind and current vector and y is the NRCS (σ 0 ) and Doppler shift (f D ). Empirical GMFs are commonly used as forward models. These GMFs are frequency dependent and parameterized by incidence angle, azimuth angle, radar polarization. A detailed description of the surface scattering and Doppler using theoretical and empirical GMFs can be found in, e.g. (Fois et al., 2015). The used GMFs are described in more details in the following.

Backscatter
In the Bragg scattering theory, the backscatter (NRCS or σ 0 ) is proportional to the spectral density of the Bragg waves (Wright, 1966). Longer waves also affect the NRCS through modulation of the Bragg waves spectrum (Wright, 1968). At, e.g. C-and X-band, the Bragg waves correspond to the gravity-capillary waves with wavelengths in the order of few centimeters. These waves are directly related to the instantaneous local surface wind stress. The backscatter is also sensitive to the relative angle between the radar antenna beam and the Bragg waves direction (wind direction). In addition, the backscatter depends on the polarization of the electromagnetic wave and the incidence angle. It was found that sea surface temperature (SST) effects on the backscatter are more important at higher frequencies (e.g. Ku-band) and almost negligible at low frequencies (e.g. C-band) (Wang et al., 2017). There is to our knowledge no study of this effect at X-band, but based on the abovementioned study, a moderate effect should be expected at X-band.
The surface currents play a minor, but non-negligible, role in the backscatter through several effects. The first effect is due to the relative motion between the air and sea surface. Imprints of surface currents in satellite-derived winds have been shown in (Kelly et al., 2001;Plagge et al., 2012). The authors suggested that satellite radars measure wind relative to the moving ocean surface. This effect depends on the relative direction of the wind and currents. The wind stress will decrease (increase) when the surface current is in the same (opposite) direction as the wind. In general, ocean surface wind speed is an order of magnitude larger than surface current speed, so this effect is usually small. It might however be important for low wind over a strong ocean current.
The second effect is due to wave-current interaction Johannessen et al., 2005). This interaction will increase or decrease the surface roughness, thus the backscatter. That is, for the same wind, a higher or lower backscatter will be measured due the presence of current gradient. This effect also depends on the relative direction of the wind and currents and on their directions relative to the radar beam. It is usually manifested in SAR images as bright or dark linear features (Johannessen et al., 1991(Johannessen et al., , 1996. Thus a direct inversion of the backscatter to wind speed, ignoring currents, might produce spurious (high or low) wind speed values.
The third effect is due to a change in the marine atmospheric boundary layer (MABL) stability by the surface temperature gradients which are usually associated with surface current systems (Chelton et al., 2004). This alteration of stability will cause a change in surface wind stress which, in turn, alters the surface roughness and consequently the radar backscatter. For instance, on the warmer water side of an oceanic front, the radar backscatter will be larger than on the colder side (Beal et al., 1997;Krug et al., 2018).
The fourth effect is due damping of the Bragg waves by surfactants accumulated in the current convergence zones (Beal et al., 1997). Usually, SAR observes a combination of these four effects and is often difficult to separate them. Depending on wind and current conditions, one might be more dominant than the other.

Backscatter forward model
The model maps the observation space (σ 0 ) into the state space (wind) as follows where u 10 is the wind speed at 10 m height and φ w is the angle between the radar look direction and the wind direction, θ is the incidence angle and pol is the polarization.
The forward model used here is an empirical GMF called XMOD2 (Li and Lehner, 2014). This GMF was built using the X-band (9.65 GHz) data from TerraSAR-X, which is identical to TanDEM-X. The variation of the XMOD2 GMF with wind speed and direction is depicted in Fig. 1. The forward model is, for a given polarization and incidence angle, a nonlinear function of wind speed and wind direction. Note that the sensitivity to wind speed is higher at up/downwind than at crosswind. Note also that the sensitivity to wind direction increases with wind speed.
As mentioned above, it is known that spaceborne radars such as SAR actually measure the relative motion between the air and the sea surface (Kelly et al., 2001;Plagge et al., 2012). Thus we, usually, assume that the current vector contributes to a difference between the SAR-retrieved wind vector and a wind vector measured by an anemometer. To take this effect into account, u 10 and φ w should be replaced by an ocean-relative wind.

Doppler frequency shift
The radar measures the Doppler shift due to the relative motion between the satellite and the Earth surface. The two common techniques used to estimate the Doppler shift from SAR data are Doppler centroid analysis (DCA) Johannessen et al., 2008;Hansen et al., 2011) and along-track interferometry (ATI) (Romeiser et al., 2010b;Elyouncha et al., 2019). The geophysical Doppler shift (also called Doppler anomaly) is calculated as the difference between the total measured Doppler shift and the calculated geometric Doppler shift. The DCA technique requires a single SAR image, while the ATI requires two SAR images of the same area.
The velocity derived from the measured Doppler shift represents the total relative (between the satellite and the ocean surface) velocity projected on the line of sight (LOS). This total velocity U D can be decomposed into three components The first term is the non-geophysical component (U geo ) which is due to the motion of the radar platform relative to the solid rotating Earth. The second component (U nwd ), which includes geostrophic currents (barotropic and baroclinic) and tidal currents, is not driven by the instantaneous local wind. For high resolution radars such as SAR, internal waves and swell also contribute with their orbital motions as local surface currents. All these currents are induced by forces independent of the local instantaneous wind.
The wind-driven component (U wd ) can be decomposed into three components where U E is the Ekman current induced by the local wind stress at the surface, U S is the Stokes drift and U wv is the integrated motion induced by wind-waves along the radar line-of-sight. The U wv term does not represent an actual water mass transport (current), but it induces a Doppler shift in the radar measurements. This means that even in the absence of currents, the radar would measure a Doppler shift. It includes the Bragg wave phase speed and a contribution induced by the correlation between the variation of the backscatter and the long waves orbital velocities (Thompson and Jensen, 1993). It was found by  that the term U wd contributes a significant amount (~30% of the wind speed) to the total geophysical Doppler shift (U D − U geo ) measured by SAR. Moreover, the U wv contribution is often larger than the surface current contribution (Elyouncha et al., 2019). Thus the observed Doppler shift is dominated by the wind induced motion. Finally, the reported values of the magnitude of the Stokes drift vary between 0.4% and 1.3% of the wind speed and directions between 0 and 45 ∘ (Ardhuin et al., 2009).

Doppler forward model
The Doppler model, here, refers to a model that predicts the contribution of the local wind and wind-waves to the SAR Doppler shift in absence of currents. This model has the following form Here, the model developed in Mouche et al., 2012), called CDOP, is used. The variation of the CDOP GMF with wind speed and direction is depicted in Fig. 2. The model is also nonlinear in wind speed and wind direction. Note that the sensitivity to wind speed is higher at up/downwind than at crosswind. The sensitivity to wind direction increases with wind speed. CDOP is an empirical model built from C-band (5.33 GHz) ENVISAT/ASAR data. Thus the predicted Doppler shift (f D mod ) must be scaled to be applicable to X-band. The CDOP model output is converted to X-band using a scaling factor as f mod where α = f r X /f r C is the scaling factor, and f r X and f r C are the radar carrier frequencies of TanDEM-X and Envisat/ASAR, respectively.

Retrieval algorithm: non-Bayesian
In the non-Bayesian approach (direct method), the wind and currents are retrieved separately. In this approach, only the wind speed and the radial component of the current can be retrieved using the backscatter and Doppler shift, respectively. The wind direction is required in the wind speed retrieval. Thus the wind speed retrieval is affected by the wind direction and the backscatter GMF errors. In the retrieval of radial current, both wind speed and direction are required to remove the wind induced contribution from the total measured Doppler shift. Similarly, the radial current retrieval is affected by errors in the wind speed, wind direction and the Doppler GMF.

Wind speed retrieval
For wind speed retrieval, the direct method is based on the maximum likelihood estimation (MLE). This is equivalent to minimizing the weighted least squares (WLS), assuming Gaussian error distribution. The cost function to minimize is expressed as in (Stoffelen and Anderson, 1997) where in this case φ wa is the a priori (known) wind direction and Δσ 0 is the standard deviation of the observed backscatter. The main  assumption in the direct method is that the GMFs and the external wind input to the GMFs are error free. If a constant standard deviation is assumed, this reduces to least squares (LS). Using this approach only one component (either wind speed or direction) can be retrieved if the other component is known (Portabella et al., 2002). Typically, the wind direction is provided by NWP model, reanalysis or in-situ measurements and the cost function is minimized for wind speed. Moreover, often for faster computation, a look-up table is constructed and the wind speed corresponding to the closest measured σ 0 is searched and interpolated.
Note that the direct method neglects the wave-current interaction, i.e. current effects are assumed to be small. Taking into account this effect would require a model that computes the spatial modulation of the wave spectrum by the current gradients, by solving the wave energy conservation equation Johannessen et al., 2005). This is very time consuming. One can however take into account the relative motion of the surface if the current information is available.

Radial current retrieval
For the radial current retrieval, the assumption is that the measured Doppler shift can be decomposed into two independent contributions, the wind-correlated and wind-uncorrelated, e.g. tidal current, contributions to f D . Similar to wind speed retrieval, the wave-current interaction is neglected. Taking into account this effect would require a model that computes the modulation of the wave spectrum by the current gradients  and then calculate the U S + U wv term which depend on the wave spectrum. The wind-correlated Doppler shift is simulated using the CDOP GMF . This GMF was built by matching global ECMWF 10 m height winds and Envisat/ASAR derived Doppler shift. Thus, the Doppler GMF most likely accounts for U S + U wv which, as mentioned earlier, are correlated to the local wind. It is, however, not clear how much of the surface Ekman current the GMF accounts for. Finally, the simulated f D mod is subtracted from the estimated geophysical Doppler shift (after removing the geometric term from the measured f D ). The difference is converted to horizontal radial velocity as where k e is the electromagnetic wavenumber. U r is assumed to be the projection of the surface current on the ground range direction.

Retrieval algorithm: Bayesian
The inversion method is based on the maximum a posteriori estimation (MAP). It is also called variational method (VAR) in the data assimilation (DA) framework. In short, the principle of the Bayesian approach is combining a priori information with the measurements, referred to as background and observations in the DA framework, assuming that all sources of information contain errors and these errors are well characterized. The objective is to find the state vector x, given the observations y, that maximizes the posterior probability p(x/y) given that the a priori P(x) and the likelihood P(y/x) are known.

Cost function
Obtaining the MAP is equivalent to minimizing a cost function (Lorenc, 1988;Rodgers, 2000). In the context of wind and current retrieval, the state vector is composed of the wind vector u w (u w , v w ) and the current vector u c (u c , v c ). The observation vector is composed of σ 0 and f D . Note that for SAR the wind and currents are spatial 2D fields, thus u w , u c , σ 0 and f D refer to the flattened 2D matrix. Thus the MAP is defined as P(u w , u c | σ 0 , f D ). The cost function to be minimized is defined, assuming Gaussian distributions, as where J o and J b are the observation and the background terms defined as where the operator H is defined as where α = 9.65/5.33, u e is the ocean-relative wind vector which has components (u w − u c , v w − v c ) and u r = u c ⋅ e l is the projection of the current vector on the radar look direction (LOS). F σ 0 and F fD are the forward models described above, I is the interpolation operator from the background grid to the retrieval grid, u w and v w are the zonal and meridional components of the wind vectors, respectively, u wb and v wb are their a priori values. u c and v c are the zonal and meridional components of the current vectors, respectively, u cb and v cb are their a priori values. S σ 0 is the sum of the covariance matrices of the backscatter observation and backscatter GMF errors. S fD is the sum of the covariance matrices of the Doppler observation and Doppler GMF errors. S wb and S cb are the covariance matrices of the background wind and current fields, respectively. Additional weak constraints can be added to the total cost function (Lorenc, 1988;Laroche and Zawadzki, 1994), which becomes J s is a smoothness term defined by the Laplacian of the wind and current components and is implemented using the finite differences approximation (Laroche and Zawadzki, 1994) where L is the Laplacian operator without the denominator because it is not relevant for the minimization. It was found that this term filters the spurious solutions and also helps the solver to converge more rapidly to the minimum. J c is a constraint that is only applied to the current components u c over land. It penalizes the deviation of u c from zero over land In practice, the two observable vectors and covariance matrices are grouped in the same vector and matrix, respectively. The same applies for background vectors and matrices. The weights w 1 and w 2 are determined heuristically.
In some studies, e.g.  and , the spatial correlation in the covariance matrices is neglected and the error is assumed to be constant. The advantage of this simplification is that, the grid points can be processed in parallel which reduces the processing time. In this study, we assume independent observations (σ 0 and f D ) but with spatially varying error. We also assume that the errors in u and v components are independent. However, we take into account the spatial correlation in the background wind and current fields (see section 2.3.3). The disadvantage of this approach is a higher computational cost but on the other hand it is more likely to retrieve physically consistent fields.
For a SAR scene of size N x × N y , the size of observation and background covariance matrices (S o and S b ) are 2N × 2N and 4N × 4N, respectively, where N = N x × N y . Finally, the minimization is done using the Python limited-memory quasi Newton routine (L-BFGS-B). The initial state for both wind and current is set to the background values. The Jacobian is computed numerically by finite differences.

Cost function analysis
In this section we illustrate the contribution of each term in the cost function, decomposed into three terms, namely backscatter, Doppler and background components Fig. 3 (upper row, panels a-e) depicts the cost function for different combinations of the three terms. The dark blue colour indicates the area where the cost function has its minima and the red areas indicate the maxima. It can be observed that using the NRCS alone J σ 0 (a) or the Doppler shift alone J fD (b), it is impossible to retrieve the full wind vector. By adding the NRCS and the Doppler terms (d), the space of the solutions is reduced to two ambiguous and symmetric minima. It can be readily noticed that without an additional information, on the wind direction, it is impossible to resolve the two ambiguities. The two cost function minima have similar values hence the two solutions have equal likelihood. By adding the background term to the NRCS (c) and to the NRCS+Doppler (e), the ambiguity is resolved. The background error makes the solution slightly deviate from the true value. Note that in the NRCS+f D +background (e) case the solution is closer to the true value than in the NRCS+background (c). Thus even in the presence of the background, there is an improvement by using the Doppler shift. The improvement is minor in this first case (upper row) because the background wind direction is very close to the true value. Fig. 3 (lower row, panels f-j) illustrates a second case where the background wind direction is biased, i.e. it deviates by 90 ∘ from the true wind direction. In this case, the role of the Doppler term is important in finding the correct solution for the wind vector. The retrieved wind direction using NRCS+background only (h) will be the wind direction of the background, thus biased. The NRCS+f D +background (j) yields the wind vector solution closer to the true wind vector. This shows that the combination of the NRCS+f D +background may yield the optimum solution, given a well calibrated NRCS and Doppler shift.

Covariance matrices
Appropriate characterization of the error statistics for both observational and background data is essential for an optimal retrieval. The covariance matrix can be decomposed as (Daley, 1991;Bannister, 2008), S = Σ ∧ Σ, where ∧ is the non-diagonal correlation matrix and Σ is a diagonal matrix with the standard deviation as the diagonal elements. We neglect the correlation between variables, because these correlations are poorly known.
Background covariance matrix S b : Estimation of the background covariance matrix has been a topic of interest for several authors in the field of atmospheric remote sensing for many years and different estimation methods have been proposed, e.g. (Bannister, 2008). It can be estimated empirically from differences between observations and background (O-B) assuming the observation errors are uncorrelated (Vogelzang and Stoffelen, 2012). This requires a large number of densely distributed observations. Alternatively, it is modeled using a structure function. The most common parameterization of the structure functions is Gaussian and exponential (Daley, 1991). Gaussian structure function was used for instance in (Danielson et al., 2008). Both functions assume homogeneous and isotropic error distribution.
We assume that errors in the wind components u w and v w are independent and normally distributed (Stoffelen, 1998). We also assume that errors in the current components u c and v c are independent and normally distributed. Furthermore, the wind vector u w and the current vector u c are assumed to be independent. Thus, the covariance matrix S b is a block diagonal sparse matrix. We take into account the spatial correlation of each variable. The spatial correlation is modeled here by the exponential structure function (Rodgers, 2000) where i and j are the grid point indices, δr is the distance between the grid points and L s is the horizontal correlation length (length scale). The correlation matrix has the so-called Toeplitz-block-Toeplitz pattern as shown in Fig. 4. The advantage of using this type of analytical structure function is that it generates a symmetric positive definite matrix. This matrix is factorized using Cholesky decomposition which makes the inversion efficient and numerically stable. The disadvantage, however, is that the value at a given grid point is affected by all the other points. This yields an ill-conditioned covariance matrix. The consequence is an increased computational cost and slow convergence.
The key factors in the covariance matrix are the standard deviation (elements of Σ) and the correlation length scale L s . The horizontal length scales of the correlation functions are based on the internal Rossby radius of deformation L R , L s = βL R , where β is between 0.5 and 2. The horizontal length scales in the atmosphere and the ocean are very different being O(100-1000 km) and O(1-100 km) in the atmosphere and ocean, respectively. In the ocean, the Rossby radius depends on latitude, stratification and thickness of the upper buoyant layer (Chelton et al., 1998). It decreases with increasing latitude and deceases with decreasing thickness (Chelton et al., 1998). In the Baltic Sea (our study area) the Rossby radius varies between 1.3 km and 7 km (Fennel et al., 1991). The correlation length for u c and v c is heuristically set to 5 km. In the atmosphere, the length scale also depends on latitude and height. Thus it is smaller in the high latitudes and near the surface. The length scale for the meridional and zonal wind components has been estimated using ASCAT scatterometer and an ECMWF model (Vogelzang and Stoffelen, 2012). However, the estimated length scale ~300 km in the above-mentioned study represent the open ocean, while our study area is coastal hence smaller length scale is expected due to topography, e.g. land shadow. Therefore, L s is set to 100 km for u w and v w .
For the wind components standard deviation, we use the values reported in literature (Portabella et al., 2002) f D is neglected. The correlation between the two variables is expected due the fact that both depend on the surface wind, which can be noticed in Fig. 6 (third and fourth columns). However, the correlation between their errors is more complicated to determine.
The observation error include instrument error, forward model error, and representation error. Instrument error, which is mainly due to thermal noise, is usually uncorrelated but the forward model and representation errors might lead to correlation. This correlation is however neglected since it is poorly known. Moreover, no evidence reported in literature showing the benefit of accounting for spatial correlation in SAR data. Thus, the covariance matrix S o used here is diagonal.
The standard deviation of backscatter measurements has been studied in literature (Portabella et al., 2002;Stoffelen and Portabella, 2006;Danielson et al., 2008;Choisnard and Laroche, 2008) and it is usually described as a fraction of the backscatter Δσ 0 = K p σ 0 , where K p is a dimensionless quantity (normalized standard deviation of σ 0 ) with typical value 0.05 < K p < 0.3 (Choisnard and Laroche, 2008). Here we follow (Portabella et al., 2002) and we set K p to 0.078.
For the Doppler shift, values of standard deviation estimated from Envisat/ASAR data are also reported in literature (~5-7 Hz (Hansen et al., 2011;Mouche et al., 2012)). In this study, the standard deviation is estimated from the TanDEM-X data over land pixels with high coherence (>0.8). The estimated standard deviations are 3.65, 3.77, 4.70 and 4.80 Hz for four satellite acquisitions. These are comparable (slightly smaller) than the Envisat/ASAR reported values. Finally, in order to account for probable forward model errors, the standard deviation is slightly inflated, thus Δ fD is set to 7 Hz.

Study area
The study area is shown in Fig. 5. It is located along the Swedish west coast in the transition area between the Kattegat and the Skagerrak Seas. The large scale circulation in this area is dominated by the Baltic current, a surface outflow (of the Baltic Sea) of low salinity water along the Swedish coast, the Jutland current (along the Danish west and north coast), the Atlantic current which is an inflow of deep saline water and the Norwegian coastal current which is an outflow along the Norwegian coast (Gustafsson, 1997;Christensen et al., 2018). The small scale circulation is more complex due to local wind forcing variation, interaction with topography, rivers run-off, tides, etc.

Background winds and currents
The wind vectors provided by a reanalysis project called Uncertainties in Ensembles of Regional Reanalyses (UERRA) are used as background. The horizontal grid spacing of UERRA is 11 km and the temporal resolution is one hour. Note that the effective resolution of atmospheric models is usually several times larger than the grid spacing (Skamarock, 2004), which prevents the accurate representation of coastal features.
The background wind fields and the observations have different grids and spatial resolutions. Unlike in data assimilation, the observation grid is used here as the retrieval (analysis) grid. Since the retrieval grid is finer than the background grid, the background winds are interpolated to the retrieval grid using bi-linear interpolation. Moreover, the background data are linearly interpolated to the satellite acquisition time. All satellite acquisitions fall between 5 AM and 6 AM (see Table 1).
The near surface currents (at − 0.5 m depth) are provided by the Nucleus for European Modelling of the Ocean (NEMO) (Gurvan et al.,Fig. 4. A subset of the correlation matrix ∧ generated using the exponential structure function.

SAR data
In this study, the backscatter and the Doppler frequency shift are derived from the ATI-SAR TanDEM-X mission data. This mission consists of two identical satellites (TerraSAR-X and TanDEM-X) flying in close formation and carrying identical SAR instruments operating at X-band. The data used in this study, called Co-registered Single-look Slant-range Complex Products (CoSSC), are provided by the German aerospace center (DLR). These images are all acquired in descending pass, VV polarization, stripmap and single receive antenna (SRA) mode with the antenna looking to the right of the flight direction. The raw spatial resolution is approximately 3 m and 3.3 m in ground range (acrosstrack) and azimuth (along-track), respectively. In order to reduce the backscatter and Doppler shift noise and to reduce the computational cost during the inversion, the resolution is degraded to 200 m in both range and azimuth. The relevant parameters of the used satellite dataset are given in Table 1.
The processed and calibrated σ 0 and f D derived from the 16 TanDEM-X acquisitions (Table 1) are depicted in Fig. 6. Every satellite scene in the figure is a mosaic of four SAR images (see Table 1) and every scene extends over an area of ~202 km in azimuth and ~34.6 km in range. The figure shows also the model (used as background) wind vectors (first column), current vectors (second column) and the SAR-estimated Doppler shift and NRCS in the third and fourth columns, respectively. The model current field (second column) shows a complex small scale circulation. The magnitude of the surface current is small (below 0.5 m/ s), except few local flows which reach 0.7 m/s. It can be observed that the backscatter spatial variation is mainly wind speed induced, i.e. high/ low backscatter correspond to strong/weak wind. A few other signatures probably due to rain cells, ships and ships wakes can be observed. The Doppler shift contains mixed signatures of wind and currents. The low order signature is wind driven. Note that the Doppler shift is generally positive (yellow) / negative (blue) when the wind is blowing toward (first acquisition) and away (other acquisitions) from the radar, respectively. In addition to the wind signature, the current adds higher order signatures which might be positive or negative. The current signature is hidden by the wind because the Doppler shift induced by the latter is much stronger (dominant) than the current induced Doppler shift.

In-situ data
Wind measurements were available from two weather stations, along the Swedish west coast, called Nordkoster A and Väderöarna A. The data are provided by the Swedish Meteorological and Hydrological Institute (SMHI). Current measurements were available from two SMHI's buoys called Nordkoster and Väderöarna. The current buoys provide the current speed and direction at five different depth levels (4, 16, 28, 40 and 60 m). The locations of the in-situ stations are shown in Fig. 5. Both weather and ocean in-situ stations provide measurements every hour. The measurements are linearly interpolated to the satellite acquisition time. The interpolated values are reported in Table 2 for every satellite acquisition. Note for instance on 2014-08-25, the wind speed increases southwards and the direction changes from north-east to north-west going from Nordkoster A down to Väderöarna A. This corresponds to the negative-to-positive Doppler shift and low-to-high backscatter, in the upper part of (c) and (d) in Fig. 6. On the other acquisitions, the wind direction is almost constant along the coast. The current direction is very different from the wind direction, indicating a low correlation, and sometimes varies between acquisitions by more than 90 o . This can not be seen directly in the Doppler shift since it is dominated by the wind signature.
Figs. 7 and 8 show the time series of wind and current, speed and direction, from August 20 to September 202014, recorded by Väderöarna A (wind) and Väderöarna (current). The wind is measured at roughly 10 m height. The plotted current is measured at -4 m depth. Both wind and current signals are highly variable in time. The correlation between the instantaneous wind and current is low. The correlation coefficient is 0.27 and 0.087 for speed and direction, respectively. The correlation coefficient increases to ~0.6, for both speed and direction, when a moving average filter of 6 days is applied. This is in agreement with (Christensen et al., 2018), who estimated the correlation time between the local wind forcing and current to be 4 and 7 days for the Jutland and Norwegain coastal currents, respectively. This also agrees with the global study in (Trindade et al., 2020) which reports NWP wind biases compared to ASCAT scatteromter with a temporal scale of ~5 days, which are attributed to ocean circulation.

Auxiliary data
These data are only used for analysis and discussion but not for the retrieval. The in-situ measured air temperature (T a ), sea surface temperature (SST) and salinity (SSS) are reported in Table 3. The air temperature is measured, at 2 m height, by the wind station (Väderöarna A). The SST and SSS are measured, at 1 m depth, by the current buoy (Väderöarna). The air-sea temperature differences indicate unstable MABL for the four acquisitions. Wave parameters significant wave height (SWH), peak period and direction are also reported in Table 3. These parameters are measured by the wave buoy (Väderöarna WR), which is located ~13 km offshore from the wind station (Väderöarna A). It can be noted that the SWH increases with wind speed but the wave direction deviates from the wind direction which might be due to refraction and diffraction by topography.
It has been shown by (Clarke and Van Gorder, 2018) that the surface Stokes drift is highly correlated with the local wind stress, hence essentially driven by the local wind and is in the local wind direction. (Clarke and Van Gorder, 2018) also showed that swell does not substantially contribute to the Stokes drift. Thus, Stokes drift can be estimated using the wind stress (Clarke and Van Gorder, 2018) or measured spectral moments (Webb and Fox-Kemper, 2011) (SWH and peak wave period T p ). (Webb and Fox-Kemper, 2011) showed that the induced error using the monochromatic approximation is smaller than the uncertainty between different Stokes drift products. In (Webb and Fox-Kemper, 2015), wave spreading loss factor (0.8) was introduced to compensate for not using a full 2D wave spectrum. Thus, the Stokes drift can be estimated as where a, ω and k are the peak wave amplitude, angular frequency and wavenumber, respectively, and g is the acceleration due to gravity. z is the water depth (z=0 at the surface and positive upward). Using the wave parameters in Table 3, the calculated surface Stokes drift is reported in the last row of the table.

Results and discussion
The retrieval of u w and u c is performed for the eight satellite acquisitions, covering the in-situ stations, using the Bayesian retrieval algorithm described above. For comparison, the wind speed and the radial Table 2 In situ measurements: Wind speed and direction at 10 m height and current speed and direction at 4 m depth. The locations of the weather stations are displayed in Fig. 5. The convention for the direction is coming from and moving toward wrt north for wind and current, respectively.  current are also retrieved using the direct method. The results are assessed quantitatively by comparing the background and the retrieval to the in-situ measurements.

Wind retrieval comparison
As mentioned earlier, SAR-measured NRCS responds to the surface roughness that is driven by the stress (τ). For practical reasons, the NRCS is related to the atmospheric real wind (u 10 ), while in fact the measured NRCS is affected by other geophysical effects than u 10 , e.g. atmospheric stability and air density. Thus, the SAR-retrieved wind using present GMFs is more related to τ, that includes all these effects, than to u 10 and is called stress-equivalent wind (u 10s ) (de Kloe et al., 2017). If the used GMF is trained on neutral wind (u n ), the retrieved wind is also neutral. Note that, the used GMF is tuned to u 10 (Li and Lehner, 2014) and the background data and the validation data both provide u 10 . To convert from u 10 to u 10s , air-sea parameters such as SST, air temperature, specific humidity and air density are required, which are not available. In this paper, the retrieved wind is assumed to be u 10 and the u 10 -to-u 10s conversion is left for future work.
For the wind validation, the SMHI stations measure the wind at 10 m height averaged over 10 min, so the comparison does not require height adjustment. However, the spatial representativeness is different between the different datasets, i.e. point against spatial average. Fig. 9 shows the retrieval, the background and the in-situ wind speed and direction for the eight satellite acquisitions for two weather stations (Nordkoster and Väderöarna). The first four acquisitions (1,5,9,13) are compared to Nordkoster and the last four acquisitions (2, 6, 10, 14) are compared to Väderöarna. The results of the comparison for the two stations are given in Table 4. It can be noticed that the wind speed retrieved using the Bayesian method (blue) is systematically closer to the in-situ measurement than both the background (green) and direct method (black), i.e. the Bayesian retrieval gives the smallest wind speed bias. An exception is the first acquisition, where the background is closer to the in-situ. During this acquisition the backscatter which is the main parameter for wind speed retrieval is very low (~ -30 dB) around the wind station. The wind speed at this acquisition is 1.36 m/s (see Table 2), which is below the validity limit of the GMFs (2 m/s).
For wind direction, the first acquisition corresponds to a crosswind situation and low wind speed in the area close to the Väderöarna station (see Figs. 5 and 6). The other acquisitions correspond to a downwind situation around the same station. As shown later in section 4.3, the Doppler effect on the wind direction bias correction is smaller in downwind situations and the wind direction error is higher at low wind speed. It can be observed from Fig. 9 (right) that the retrieved wind direction is very close to the background. This suggests that the retrieved wind direction using the Bayesian method is biased toward the background in agreement with (Portabella et al., 2002). Since the background wind direction is accurate around the wind stations in most cases, the retrieval improvement is minor. This is in agreement with the cost function analysis in section 2.3.2. Note that the background winds (UERRA) are obtained from a data assimilation system, which may explain the accurate wind direction. The retrieval method should be tested in cases where the background wind direction is inaccurate. Finally, the bias values averaged over the eight acquisitions are reported in Table 4. These values indicate an error reduction in wind speed of ~1.14 m/s averaged over all cases. The error reduction in wind direction is insignificant.

Current retrieval comparison
Comparison of satellite data, model data and in-situ measurements, for oceanic parameters, is a delicate task. First, because these data represent different physical components of the sea current. For instance, the SAR-derived current includes, in principle, all surface current components (see section 2.1) while the ocean model does not include the wave-induced Stokes drift. Second, because they represent the current at different depths. The SAR senses only the surface, i.e. the upper few millimeters. The ocean model, provides the depth averaged current in the upper meter. The buoy measures the instantaneous current at a specific depth (at − 4 m). Third, the different data represent different horizontal resolutions. SAR and the ocean model report the average over the pixel which are, here, ~200 m×200 m and ~1850 m×1850 m, respectively. The in-situ data report a point measurement.
Theoretically, the vertical profile of the magnitude of the Ekman current and the Stokes drift follow an exponential law, and can be approximated as (Clarke and Van Gorder, 2018) where z is the water depth, i stands for Ekman or Stokes, L S and L E are the Stokes and Ekman e-folding depth. The Stokes"e-folding"depth is about 1 to 2 m, e.g. (Axell, 2002). Thus, the Stokes drift should be negligible at − 4 m depth, even if the buoy was affected by the waves orbital motion. The ocean model (NEMO-Nordic), that generated the background currents, does not account for the Stokes drift. Thus, in principle only SAR reports the Stokes drift but the Doppler GMF removes it, since this GMF is fitted using the local wind . The e-folding scale of the Ekman is much larger, L E is between 10 and 100 m. It increases with wind speed and decreases with latitude. For instance, at a latitude 55 N and U 10 =10 m/s, the Ekman depth is about 10 m. Thus, in theory, the The tidal current decreases with depth, but with a lower rate following a power law (Lewis et al., 2017) where β and γ are parameters of the model, z is the height above sea bed, h is the total water depth and U is the depth averaged velocity. For typical values of γ and β (0.32 and 7), the difference between the surface and 4 m depth is about 1%, hence negligible. Thus, SAR, the ocean model and the buoy will be affected by the tidal current with slightly decreasing magnitude. The current shear in the last meter below the surface was investigated in (Laxague et al., 2018). The author shows that currents in the upper few centimeters of the ocean may have drastically different magnitudes and directions than the average over the upper meter. Thus, usually in-situ observations are extrapolated to the surface. However, with our in-situ data, we don't have access to the last centimeters to investigate the shear. Fig. 10 shows the vertical profile of the current speed and direction, respectively, before (at 05 AM) and after (at 06 AM) the satellite acquisition. A high temporal variability of the profile can be noticed. The slope of the current speed near the surface varies within the hour. Note for instance, the change in 2014-08-30, at -4 m, by about 0.12 m/s between 05 AM and 06 AM. The current direction profile also varies within an hour around the acquisition time. Note, for instance, the current direction measured by the buoy on 2014-09-05 at -4 m depth has changed from 135 deg. at 05 AM to 80 deg. at 06 AM, i.e. 55 deg. difference. Furthermore, the current profiles do sometimes not follow the Ekman spiral, e.g. lower speed at -4 m than at -16 m or direction turning to the left. Given this variability of the current profiles, it is difficult to extrapolate the in-situ current speed and direction to the surface. Therefore, the current retrieval is validated against the in-situ measurement at -4 m as it is, keeping in mind all these differences discussed above. That is, according to the discussion above, it is expected that the surface currents provided by SAR and the model should be in general larger than the in-situ measured current. Fig. 11 shows the retrieval, the background and the in-situ current speed (CS) and direction (CD) for the eight satellite acquisitions. First, note the generally low in-situ current speed (≲0.2 m/s) during most acquisitions. This is just above the limit of the retrieval uncertainty (see section 4.3). For this comparison, the retrieval is averaged within 1 km radius, which is about half the model grid spacing. It can be observed from Fig. 11 that the retrieval gives lower bias than the model in both current speed and direction at most acquisitions to different extents, i.e. the bias difference is sometimes significant and sometimes negligible. Similar to the wind direction retrieval, the CD is clearly imposed by the model. Encouragingly however, the retrieval has a non negligible impact on the current direction, i.e. it is never worse than the model. The bias values are given in Table 5. The values provided in the table indicate that, on average, the model bias is reduced by ~0.04 m/s and ~12 o in CS and CD, respectively. Moreover, the comparison between the direct method and the Bayesian method is illustrated in Fig. 12, by comparing the radial current (U r ) to the Bayesian retrieved current projected on the SAR look direction. It can be observed that the Bayesian retrieval gives the lower bias than both the background and the direct method.

Case studies
The comparison against in-situ data is local which does not illustrate the spatial variability of the wind and current fields. In this section, we analyze the retrieval of atmospheric and oceanic features through two case studies to show the benefit of blending models with SAR data on the spatial representation. In order to make the direct method, background and Bayesian retrieval comparable, only the wind speed and the radial current are compared. Thus, the background and the retrieval current vectors are projected on the SAR look direction. Note that the satellite is flying from north to south and looking to the right. Thus in the following figures, negative (blue) values of U r indicate a flow from east to west and vice-versa. Finally, the figures are depicted in the satellite geometry range (across-track) and azimuth (along-track). Fig. 13 (upper row, panels a-c) shows the wind speed retrieved from SAR using the direct method (panel a), provided by the background (panel b) and retrieved using the Bayesian method (panel c). The direct Table 4 Comparison of wind speed and direction against in-situ measurements. The values represent the bias (background, direct and Bayesian retrieval minus in-situ) for the eight satellite acquisitions.  method shows a strong and sharp front with wind speed difference of ~6 m/s between the two sides of the front. This is due to the very high sensitivity of the backscatter to wind speed. On the other hand, the front is not well resolved by the background. It shows, however, a shallow wind speed gradient in the north-east south-west direction. The Bayesian retrieval (panel c) slightly smoothes the wind speed variation but preserves well the location and magnitude of the front. This is an example where the Bayesian method has successfully combined both SAR and background wind speed information.

Case study 1
For the current (lower row, panels d-f), the direct method and the background radial current show to some extent similar features with different magnitudes. The front is still visible in the direct method radial current for two reasons. The first reason is that the background wind speed and direction are used to remove the wind-induced Doppler shift (see section 2.2). These winds do not carry a strong front signature as shown in the first row. The second reason is that an atmospheric front might induce an oceanic surface current which has a component in the radial direction. The Bayesian retrieval of radial current (panel f) lies between the direct method and the background radial current. In this case the background slightly mitigates the front signature in the current. The SST and SSS obtained from the Nemo-Nordic model are depicted in Fig. 14. This figure does not show a clear oceanic front, which suggests that the wind speed gradient is likely to be due to an atmospheric front.
For comparison, wind speed is extracted from the Metop-B ASCAT Level 2 coastal product. The spatial grid spacing of this product is 12.5 km. The wind speed is interpolated into the SAR grid and extrapolated toward the coastline for the acquisition 2014-08-30. The interpolated wind speed is plotted for the whole SAR frame in Fig. 15. The time of ASCAT acquisition is 08:27 that is approximately 3 h after the TanDEM-X acquisition. The time difference is too large for a quantitative (pixelby-pixel) comparison. Qualitatively, despite the coarse resolution of ASCAT compared to SAR, it can be observed from Fig. 15 that the front is persistent. Note that the front is tilted and displaced toward the north and the wind speed has increased on both sides of the front compared to SAR acquisition. This is due to the temporal variation of the wind during the 3 h time lag in agreement with Fig. 7, which shows an increase of wind speed from ~8 m/s (at SAR acquisition time) to ~12 m/s at ASCAT acquisition time.
Finally, the effect of the retrieval on the wind speed gradient across the front is illustrated in Fig. 16. For simplicity, the azimuth profile is averaged over the range direction. It can be observed that the direct Fig. 11. Comparison of the current speed (left panel) and current direction (right panel) against in-situ measurements for the eight satellite acquisitions. Insitu (red), background (green) and the Bayesian method (blue). The current direction convention is trigonometric, i.e. 0: toward east, 90: toward north, − 90: toward south. For the correspondence between acquisition number and acquisition date, see Table 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 5 Comparison of current speed and direction against in-situ measurements. The values represent the bias (background, direct and Bayesian retrieval minus in-situ) for the eight satellite acquisitions.  method shows the sharpest gradient with change in wind speed of ~6 m/s across the front. The background shows a smaller change (~3 m/s) and the retrieval lies between the SAR and the background. For comparison, the wind speed profile extracted from ASCAT data is also plotted. The latter shows a displaced front with a wind speed change of ~ 5 m/s, which is between the direct method and background values.

Case study 2
In this section, we show an example illustrating a fragment of the Baltic Sea fresh water outflow along the Swedish west coast. This outflow from the Kattegat into the Skagerrak Sea which joins the Norwegian current is an important feature in the Baltic Sea surface circulation (Gustafsson, 1997). Fig. 17 depicts the retrieval results, for this case study, of the wind speed (upper row, panels a-c) and radial current (lower row, panels d-f). The wind speed of both direct method (panel a) and the background (panel b) are similar on average, varying within the same range (between 0 and 7 m/s). The spatial variation is however very different. The background wind speed (panel b) is smooth with a shallow southwest-northeast gradient and wind speed decreasing toward the coastline. The wind speed derived using the direct method shows a frontal signature, which is probably due to the oceanic front associated with the coastal current (see also Fig. 18). Note for instance the sharp decrease of the wind speed at the current boundary. This is likely due to the current-induced SST front effect on the SAR backscatter and consequently on the retrieved wind speed, discussed earlier (see section 2.1 and compare Fig. 17 with Fig. 18). The retrieval preserves the large scale gradient and mitigates the current signature introduced by SAR. Residual current signatures can still be observed in the retrieval. This suggests that the wind speed is strongly imposed by the SAR observables.
The retrieval result for the radial current is depicted in the lower row of Fig. 17. The figure shows the radial current derived from SAR using the direct method (panel d), from the background ocean model (panel e) and retrieved using the Bayesian method (panel f). There are some  A. Elyouncha et al. similarities and differences between the direct method and the background radial current. The main similarity is the strong coastal outflow (negative U r , dark blue) with a magnitude about ~0.6 m/s in the core of the flow. It can be observed that the current derived with the direct method shows small scale circulation features such as the meandering along the front (also visible in wind speed), which is typical of baroclinic surface currents. The retrieved radial current is a combination of the features present in both the background and observation. It preserves the large scale circulation (outflow core) provided by the model and introduces small scale features provided by SAR and not resolved by the model. The acquisition 2014-09-05 coincides with a favorable wind direction for the Baltic outflow as discussed in (Gustafsson, 1997), i.e. the wind has been blowing from the south (150 o -200 o ) for few days preceding the satellite acquisition.
Note from Fig. 17 (lower row) that the current buoy Väderöarna (red circle) is, according to the direct method (panel d), located just outside the outflow (in white-red area). This agrees with the in-situ recorded low current speed (0.065 m/s) roughly in the opposite direction (northwestto-southeast) of the outflow (see Table 2), depicted by the red arrow. On the other hand, according to the background, the buoy is inside the outflow, hence the high current speed and southeast-to-northwest direction (see also Fig. 11). This is an example of how the high resolution information introduced by SAR is determinant for the current retrieval in complex coastal circulation.
For comparison, the SST is extracted from the Metop-A/AVHRR product. The SST is interpolated into the SAR grid and extrapolated toward the coastline for the acquisition covering the outflow (2014-09-05). The interpolated SST is plotted in Fig. 18 (panel c). Note that the time of AVHRR acquisition is 00:00, that is approximately 5.40 h before the TanDEM-X acquisition, thus some mismatch is expected. The SST values provided by AVHRR are very close (slightly higher) to the model SST and to the values recorded by the in-situ buoy on 2014-09-05 (see Table 3). Despite the coarse resolution and lag time, it can be observed, from Fig. 18, that SST spatial pattern is roughly similar to the current pattern present in the retrieved (compare Fig. 18 with lower row of Fig. 17). Note that the cold side of the front corresponds to lower SAR wind speed and vice versa (compare Fig. 18 with panel a of Fig. 17). The contours do not match exactly due to the time lag and due to the fact that correlation between SST and wind is not immediate, e.g. (O'Neill et al., 2003).
Finally, Fig. 19 depicts the wind speed variation and the SST gradient across the front. For simplicity, the azimuth profile is averaged over the range direction. Note that there is an SST change of about ~0.6 o Celsius between the core of the outflow and the surrounding area. The SAR wind speed (solid line) shows a change of about 3 m/s across the front. The gradient of the SST and the SAR wind speed are shifted, probably due to the time lag. The model (dashed line) shows a shallower decrease of wind speed toward the coastline. The SST gradient effect is removed in the Bayesian retrieved wind speed (dash-dotted). Thus, in this case the SAR observes the SST effect on ocean stress, where effects of SST variation on this scale are not taken into account in the background model winds at 10 m.

Estimation and analysis of retrieval uncertainty
In the linear problem, i.e. using a linear forward model, the retrieval error can be calculated analytically using the MAP covariance expression (Rodgers, 2000). In a non linear problem, the practical approach is to characterize the retrieval uncertainties numerically, by Monte Carlo simulation. Moreover, the errors in the retrieved wind and current usually depend on the algorithm used for the retrieval, e.g. constraints, and thus can only be partially characterized by simulation.
The simulation is conducted as follows, a large number of "true" wind vectors and current vectors u t are generated. The background fields are formed by adding a random noise ε b according to S b to the true vectors u b = u t + ε b . The measurements are simulated using the forward models y = F(u t ), and a random noise, according to S o , is added to the measurements simulated from the non-perturbed "true" winds y n = F (u t ) + ε o . The noisy measurements together with the background fields are inverted using the same retrieval algorithm applied to the data, û = map ( y n , u a ) , where map refers to the inversion operator. The difference and the standard deviation of the difference (û − u t ) represent the bias and the rms error, respectively. The bias in both the background and measurement is assumed to be zero.
The nonlinearity of the forward models, yields a dependence of the sensitivity of the observables σ 0 and f D on the wind speed (WS) and wind direction (WD). This consequently yields a modulation of the error (bias and rms) as a function of WS and WD (Stoffelen and Portabella, 2006). The sensitivity of the GMFs as a function of WD is depicted in Fig. 20. Thus, the error and its variation as a function of WS and WD are analysed separately in the following sections.

Effect of the Doppler term on the wind bias
It is known (Portabella et al., 2002;Choisnard and Laroche, 2008) that the Bayesian retrieval of wind speed from SAR has WD-dependent error. This is due to the variation of sensitivity of the GMF (F σ 0) with wind direction (see Fig. 20). The sensitivity of the GMF is lowest at up/  down and cross wind. Moreover, the sensitivity to the wind speed is lowest in the crosswind direction (see Fig. 1). This yields the worst wind speed retrieval performance in the crosswind. Fig. 21 (left panel) shows the simulated wind speed bias as a function of the true wind direction, with and without the Doppler term. The systematic mean bias about − 0.1 m/s is due to nonlinear relation between the wind components and the wind speed, i.e. zero-mean random errors in (u,v) domain induce a bias in the (WS,WD) domain (Stoffelen, 1998). Note that the bias is higher, in absolute value, in the crosswind direction which is due to the GMF nonlinearity effect mentioned above.
It can be observed that inclusion of the Doppler term (with σ fD =5 Hz) in the cost function, flattens the bias curve. This is due to the fact that the Doppler GMF (F fD ) compensates for the lack of sensitivity around the crosswind direction as shown in Fig. 20. When we increase/decrease the background and Doppler uncertainties, to 3 m/s and 3 Hz, respectively, the bias around the crosswind continues to decrease but the bias at up/ down wind slightly deteriorates.
It is also known (Portabella et al., 2002;Choisnard and Laroche, 2008), that the Bayesian retrieval of wind direction is compromised by the background, i.e. the wind direction is imposed by the background. Fig. 21 (right panel) shows the simulated wind direction bias as a function of the true wind direction, with and without Doppler term. In this simulation the background wind direction is biased by 20 o . It can be observed that, without the Doppler term the retrieved wind direction follows the background. When the Doppler term is introduced in the cost function, the retrieved wind direction is corrected toward the true wind around the crosswind where the sensitivity of Doppler GMF is highest (see Fig. 20). Closer to the up/downwind directions, where both GMFs have low sensitivity, the retrieved wind direction is still biased toward the background. If we increase/decrease the background and Doppler uncertainties, to 3 m/s and 3 Hz, respectively, the area around the crosswind broadens and draws the retrieval toward the true wind. To conclude, Doppler GMF (F fD ) compensates for the lack of sensitivity around the crosswind direction but it does not help much at up/down

Error wind speed dependence
In this section, we estimate the error of the wind and current speed and direction as a function of the true wind speed. Fig. 22 shows the estimated error (rmse) as a function of the true wind speed. The wind speed error (left panel) increases slightly with wind speed but remains below 1.25 m/s up to 20 m/s, which reflects the direct dependence of the WS error on the backscatter error. The wind direction error (right panel) is higher at low wind speed but it drops rapidly at WS~ 3 m/s and remains below 20 deg. up to 20 m/s. This is because the GMFs direction sensitivity increases with wind speed (see Figs. 1 and 2).
The current speed (CS) and current direction (CD) errors are only weekly dependent on wind speed, being below ~0.15 m/s and ~20 o , respectively. They are slightly worse at the lowest and highest wind speeds. Note that the current speed error is mainly sensitive to the Doppler error, it is only weekly dependent on the backscatter through the effective wind (see Eq. 9). We have found that the current retrieval improves with increasing incidence angle (not shown). This is due to larger projection on the radial direction which increases sensitivity and also due to lower contribution of the wind-induced Doppler shift. Fig. 23 shows the error (rmse) of the wind and current speed and direction retrievals as a function of the true wind direction. The wind speed error (left panel) varies around 0.5 m/s with a small modulation as a function of wind direction. This reflects the GMF modulation discussed above. The current speed error is almost constant, with no strong modulation as a function of wind direction. It varies between 0.08 and 0.14 m/s with slight bumps around crosswind. The WS and CS errors are mainly imposed by the backscatter/background and Doppler/background uncertainties, respectively. The backscatter plays a major and a minor role in the WS and current retrieval, respectively.

Error wind direction dependence
The wind direction error (right panel) varies around 10 o with minima at crosswind and maxima at up/down wind directions, reflecting the GMF nonlinearity effect discussed above. Finally, for the current direction, the error varies around 15 o with no apparent modulation.

Conclusion
We have proposed and demonstrated a joint retrieval of ocean surface wind and current vector fields using single-beam single-polarization SAR data and background wind and current vector fields. The SAR observables used in the retrieval are the backscatter and the Doppler frequency shift. The backscatter is calibrated (including noise subtraction)   using the TanDEM-X in-product provided calibration factors. The Doppler shift is estimated from the TanDEM-X ATI-SAR interferogram and calibrated using land as a reference, after removing the topographic phase from the interferogram. The retrieval algorithm is based on the Bayesian approach (MAP estimation). The background wind and current fields are provided by a reanalysis product (UERRA) and an ocean circulation model (NEMO-Nordic), respectively.
The retrieved wind and current speed and direction were compared to in-situ measurements. It was found that the Bayesian retrieval method gives the lowest bias compared to both the direct method (SAR only) and the background (model only). On average, over eight satellite acquisitions, the bias was reduced by 0.04 m/s and 12 ∘ in current speed and direction, respectively. The wind bias reduction varies slightly between buoys, i.e. it depends on wind, but on average over the eight acquisitions it is about 1.14 m/s for wind speed. The average bias reduction in wind direction is minor in most cases since the wind direction provided by UERRA is quite accurate. Moreover, it was shown that the retrieval improves the representation of the spatial variation of the wind and current fields by introducing small scale features unresolved by the atmospheric and oceanic models. This was illustrated by examining two case studies.
The first case study is a cyclonic front over the transition area between the Kattegat and the Skagerrak Sea. This front was characterized by a sharp change in wind speed and direction. In this case, the atmospheric model alone was unable to locate the front precisely, it only shows a shallow gradient of wind speed decreasing toward the coastline.
On the other hand SAR shows a very clear and sharp front in the backscatter, the Doppler shift and the retrieved wind speed using the direct method. This case study illustrates the benefit of blending SAR with atmospheric models. For comparison, the wind speed from Metop-B/ ASCAT acquired about 3 h later is analysed. Though the absolute values of the wind speed and the location of the front have changed between the two acquisitions, a sharp spatial variation of the same order (~5 m/s) as recorded by SAR is observed.
The second case study represents an oceanic front induced by a fragment of the Baltic Sea surface outflow along the Swedish west coast. The estimated wind speed using the direct method shows a clear current signature. This signature is absent in the background, which is due to the fact that the atmospheric model is not affected by the SST changes. The backscatter modulation by the current gradient via wave-current interaction is usually manifested as bright/dark linear features along the front (Johannessen et al., 1991(Johannessen et al., , 1996. This is not strongly visible in our case suggesting a minor effect of the current gradient. This is probably due to the fact that the current and the wind have roughly the same direction. The relative motion effect is ideally taken into account by introducing the ocean-relative wind in the cost function. Thus the  remaining explanation is the effect of the SST gradient on the MABL. This was investigated by examining the SST acquired by Metop-A/ AVHRR over the same area. Investigation of the SST suggests indeed a correlation between the SST and wind speed gradient. It is shown that the Bayesian retrieval reduces the modulation of the wind speed induced by the SST variation. For the current, the background and the direct method agree roughly on the large scale circulation. Both indicate a strong, compared to the surrounding, surface outflow of ~0.6 m/s but differ on the small scale variation. SAR, for instance, shows a meandering current with variable width which is a typical feature of baroclinic currents. While in the background, the current width and magnitude are almost constant over larger distance. The Bayesian retrieval preserves the large scale circulation imposed by the model and introduces a smoothed version of the small scale variation, e.g. meandering.
An analysis of the retrieval uncertainty was conducted. First, the effect of the Doppler term on the wind bias as function of wind direction was analysed. It was found that this term is particularly important at crosswind direction, where the backscatter GMF sensitivity is low, and in cases where the background wind direction is biased. Error analysis shows that an uncertainty below 1.5 m/s and 0.15 m/s can be achieved for the wind and current speed, respectively. The current speed is slightly worse at low wind speed, mainly due the degraded wind direction performance. The achievable wind and current direction errors are below 20 ∘ . The wind direction performance improves with increasing wind speed. The current direction performance is slightly worse at the lowest and highest wind speeds.
In this paper, the used GMF was tuned to real winds (u 10 ), the used background and the validation data both provide u 10 , hence the retrieved wind is assumed to be u 10 . This can be improved by converting the buoys and background winds to stress-equivalent wind (u 10s ) following (de Kloe et al., 2017) and using a GMF that is trained on u 10n or u 10s (if available). This corrects for the atmospheric stability and air density effects before the retrieval. The wind retrieved in this way should be the stress-equivalent wind, which can be converted to real wind or stress for meteorological and oceanographic applications, respectively.
The validation study has been restricted to eight satellite acquisitions, due to the limited number of acquisitions collocated with in-situ measurements. Therefore, further validation of the wind and current fields retrieval using extended data is necessary. It is necessary to explore, in more details, the impact of oceanic and atmospheric parameters other than wind, e.g. SST, SSS, atmospheric stability and air density, on the retrievals. It is also necessary to investigate the dependence of the method on the background accuracy by e.g. comparing retrievals based on different backgrounds. The proposed retrieval method is directly applicable to other SAR data, e.g. Sentinel-1, and also applicable to future missions such as a Doppler scatterometer. Finally, the accuracy and high spatial resolution of SAR data suggest a potential benefit of assimilating these data in air-sea coupled models. Note however, that the poor temporal resolution of SAR limits its capability for tracking fast evolving atmospheric features. On the other hand, SAR is efficient in tracking ocean features which require less temporal coverage. The temporal coverage can be improved by launching more SAR missions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.