Fabrice Ardhuin
email: fabriceardhuin@ifremer.fr

Matias Alday

Maria Yurovskaya

Total Surface Current Vector and Shear from a Sequence of Satellite images: Effect of Waves in Opposite Directions

Phase shifts in lagged pairs of satellite images yield

Introduction

Surface current velocities play an important role in many ocean processes, including the flux of kinetic energy from the atmosphere to the ocean [START_REF] Wunsch | Ocean circulation kinetic energy:reservoirs, sources, and sinks[END_REF], air-sea fluxes [START_REF] Cronin | Air-sea fluxes with a focus on heat and momentum[END_REF], and the transport of buoyant material [START_REF] Maximenko | Toward the integrated marine debris observing system[END_REF]. Different observation systems have been proposed to measure the surface current in a wide range of conditions. [START_REF] Barrick | Extraction of wave parameters from measured HF radar seaecho Doppler spectra[END_REF] and many others have developed land-based HF radars that rely on the dispersion relation of surface gravity waves, while open ocean conditions are very sparsely monitored by a wide range of techniques that differ in their effective depth of measurement, as illustrated in Fig. 1. In situ moorings are typically limited to measurements at depths larger than 5 m, away from the layer where the Stokes drift of surface gravity waves is strong. In particular, Surface Velocity Program (SVP) drifters have been designed to have the least influence of wave motions in their measurements thanks to a drogue centered around 15 m depth [START_REF] Niiler | Wind-driven motions in the Northeast Pacific as measured by Lagrangian drifters[END_REF][START_REF] Lumpkin | Advances in the application of surface drifters[END_REF]. In the absence of that drogue, the drifter measures a not so clear combination of wind and surface cur-rent speeds [START_REF] Elipot | A global surface drifter data set at hourly resolution[END_REF]. The surface drifts of Argo floats have also been used [START_REF] Lebedev | Yomaha07: velocity data assessed from trajectories of argo floats at parking level and at the sea surface[END_REF], and, for the lack of a better alternative, satellite remote sensing can be used, combining scatterometer winds, sea level anomalies from altimeters, and a combination of drifters and satellite gravimeters for the Mean Dynamic Topography [START_REF] Rio | Beyond GOCE for the ocean circulation estimate: Synergetic use of altimetry, gravimetry, and in situ data provides new insight into geostrophic and Ekman currents[END_REF]. DopplerScatt [START_REF] Rodríguez | Estimating ocean vector winds and currents using a Ka-band pencil-beam doppler scatterometer[END_REF], CARTHE drifters [START_REF] Novelli | A biodegradable surface drifter for ocean sampling on a massive scale[END_REF], HF radars at 12

MHz (Stewart & Joy, 1974), near nadir Ka-band radars such as KaRADOC [START_REF] Marié | Measuring ocean surface velocities with the kuros airborne near-nadir doppler radar: a multi-scale analysis in preparation of the skim mission[END_REF].

The depth of measurement for wave dispersion is taken as the depth where a linear velocity profile matches the contribution of the current to the phase speed, namely z = 0.08L where L is the considered wavelength (Stewart & Joy, 1974;[START_REF] Teague | A comparison of multifrequency HF radar and ADCP measurements of near-surface currents during COPE-3[END_REF].

These estimates of the near-surface current can have significant differences, in part due to the sampling of different depths as illustrated in Fig. 1. Each measurement system provides a horizontal current velocity that is a convolution of the vertical profile of the velocity. For simplicity, it is convenient to define a "measurement depth" that can be taken as the depth at which a linearly varying current takes the given value. We note that DopplerScatt involves an empirical Geophysical Model Function and thus the physics of the measurement are not completely understood but the backscatter dominated by short gravity waves suggests a measurement depth under 0.1 m, whereas near-nadir radar measurements, such as performed by the KaRADOC instrument [START_REF] Marié | Measuring ocean surface velocities with the kuros airborne near-nadir doppler radar: a multi-scale analysis in preparation of the skim mission[END_REF] give a velocity that is weighted by the surface slope spectrum and corresponds to a measurement depth that does not vary much around 1 m. It is thus desirable to measure the vertical shear of the current in order to be able to compare or combine these estimates from different observing systems. The shear is also an important indication of mixing or lack thereof, giving information on possible upper ocean stratification.

Shear estimates have used the wave dispersion modification due to the current vector, defined by the two components U x (z) and U y (z) of the horizontal current profile (Stewart & Joy, 1974). For completeness, a non-linear wave correction should also be included [START_REF] Broche | Mesure par radar décamétrique cohérent des courants superficiels engendrés par le vent[END_REF][START_REF] Ardhuin | Observation and estimation of Lagrangian, Stokes and Eulerian currents induced by wind and waves at the sea surface[END_REF], which is almost the same as replacing the Eulerian mean current by the Lagrangian mean current [START_REF] Andrews | On wave action and its relatives[END_REF]. We thus expect, for kD 1,

U (k, ϕ) U (k) cos(ϕ -ϕ U ) = 2k 0 -D U x (z) exp(2kz) cos ϕ + U y (z) exp(2kz) sin ϕdz. (1)
Obtaining current shear from a sequence of images has been done from many sensors including stereo-video imagery [START_REF] Fedele | Space-time measurements of oceanic sea states[END_REF], X-band radar [START_REF] Campana | The development of an inversion technique to extract vertical current profiles from X-band radar observations[END_REF] or polarimetric imagery [START_REF] Laxague | Observations of near-surface current shear help describe oceanic oil and plastic transport[END_REF]. In all cases it requires reliable estimates of U (k, ϕ), for different wavelengths, including the shortest components, and this is performed by identifying propagating waves in the three-dimensional (3D)

Fourier transform of the measured signals [START_REF] Young | A three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents[END_REF][START_REF] Peureux | Note on the directional properties of meter-scale gravity waves[END_REF].

A great opportunity is offered by satellite imagery with accurately co-registered views of the same ocean surface with short time lags. This is particularly the case of Sentinel 2 imagery has been used to estimate surface current [START_REF] Kudryavtsev | Sun glitter imagery of surface waves. part 2: Waves transformation on ocean currents[END_REF]. The Sentinel 2 Multispectral Instrument [START_REF] Drusch | Sentinel-2: Esa's optical high-resolution mission for gmes operational services[END_REF] has very strict co-registration requirements that make it possible to observe the signature of current velocities of the order of 1 m/s [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF]. Compared to methods that use a series of many images processed with a 3D Fourier transform, the analysis of only a few images is more difficult because of the very poor temporal resolution that does not give a full spectrum in the frequency domain. In particular the linear wave signal is not so easily separated from other contributions to the measurement.

The objective of the present paper is to discuss the influence of this limited time sampling on the accuracy of surface current estimates, in the presence of waves propagating in opposing directions, starting with the 2-image method used by [START_REF] Kudryavtsev | Sun glitter imagery of surface waves. part 2: Waves transformation on ocean currents[END_REF], as discussed in Section 2. In order to demonstrate the different processing steps and the influence of the image properties, we rely on the comparison of true data and simulated images generated using the simulator described in Appendix A. Due to the possible corruption of phase speeds by waves in opposing directions, we

propose a new method using sequences of 3 images, as described in Section 3 with details given in Appendix B. Discussions and conclusions follow in Sections 4 and 5. This paper does not address issues associated to systematic errors in the spatial registration on a global reference system with sub-pixel accuracy. These are partly discussed in [START_REF] Kääb | Glacier remote sensing using Sentinel-2. part I: Radiometric and geometric performance, and application to ice velocity[END_REF] and [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF] and will be the topic of future work.

2 Effect of waves in opposite directions with 2-image sun glint method

Short waves in opposing directions

Pictures of the sun glint reveal wave patterns that are caused by the tilting of the sea surface by waves with wavelength larger than the pixel, adding their long wave slope to the local slope probability density function, and thus changing the pixel brightness. This effect has been described in many papers including Kudryavtsev et al. (2017a), and the geometry of the measurement is defined in Fig. 2. A key concept is that the surface can be decomposed in facets with a size of the order of 1 mm by 1 mm, scale at which the sea surface is well approximated by a plane. There are thus a large number of such facets in a typical image pixel (10 m by 10 m for some of the bands of the MSI sensor on Sentinel 2) but the number of those that correspond to the specular direction can be relatively small, of the order of 100, while their brightness also varies, introducing random fluctuations in the image brightness.

As shown in Fig. 2.b for a spherical Earth, the satellite position S and observation point O correspond to a zenithal angle θ v , related to the off-nadir angle γ by the law of sines, sin

γ/R E = sin(π -θ v )/(R E + H). (2)
The time of acquisition of the different pixels is not available in the Level-1C Sentinel 2 product, but it can be retrieved from the provided view geometry. For example color band B01 is acquired at time t 1 when B02 is acquired at time t 2 , the time difference is given by the ratio of the angular distance α 1,2 between the two nadir points N 1 and N 2 , as depicted in Fig. 2.c, and the angular speed along the orbit Ω (in rad/s). The angular distance α 1,2 is obtained from the law of cosines on the sphere, cos α 1,2 = cos α 1 cos α 2 + sin α 1 sin α 2 cos(ϕ 2 -ϕ 1 ).

θ v θ sun C N O S R E R E θ v (a) (b) H γ (c) α N 1 N 2 O α 1 α 2 α 1,
(3) This typically gives distances and time lags within 1% of the expression given by eq.

(1) in [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF].

In order to illustrate the limitations of the 2-image method, we start from the same image example that was used in Kudryavtsev et al. (2017a), acquired off the California coast in the region of San Diego. The image processing method is illustrated in Fig. 3. In order to understand the processing results, we also have generated simulated images and applied the exact same processing to the them.

The image simulator is described in more detail in Appendix A, and corresponds to the forward model of Kudryavtsev et al. (2017a), combined with a noise model. For our first example, the model input parameters are the Sentinel 2 viewing geometry, an estimate of the surface wind vector given by satellite scatterometer data, and a directional wave spectrum that is estimated from an in situ buoy. The buoy is station number 220 of the Coastal Data Information Program (CDIP) located at 32.752N

117.501W, also identified by the World Meteorological Organization with the number 46258.

In order to obtain a more robust estimation of the current speed, we used a phase estimated from the coherent sum of the complex amplitudes obtained from individual image tiles that are 500 m wide. We first sum the 16 2 =256 tiles, and then add 15 2 tiles that are shifted by 250 m in each direction in order to use the signal that is otherwise much reduced by the 2-dimensional Hann window, but is these shifted windows are not independent samples. This gives 256 independent estimates of the phase and associated current for each spectral component. The main information that we shall use is the phase difference ψ i,j (k, ϕ) between the bands number i and j, for the spectral component with wavenumber k propagating in azimuth ϕ. By definition the phase changes sign for opposite directions ψ i,j (k, ϕ) = -ψ i,j (k, ϕ+ π), and this phase difference can be interpreted as the result of the propagation of a single wave component at speed

C(k, ϕ) = ψ i,j (k, ϕ)/[k(t i -t j )], (4)
where t i is the acquisition time for band i. Most of our analysis of Sentinel 2 data is based on the bands B04 and B02 giving a phase ψ 4,2 (k, ϕ).

The shortest waves that propagate along the x or the y axis in the image have a 20 m wavelength. Their phase speed, for zero current, is expected to be 5.6 m/s and thus the wave patterns should be displaced by 5.6 m between the red and the blue channels that are separated by 1.0 s, and only 2.8 m between the red and green. This distance is shorter than the 10 m pixel size, but comparable to the requirement for co-registration of the MSI sensor set to 3 m for 3 standard deviations [START_REF] Drusch | Sentinel-2: Esa's optical high-resolution mission for gmes operational services[END_REF] In other words, Sentinel 2 can detect wave motions, but is the instrument and processing accurate enough to detect the generally smaller variations in phase speed associated to currents? In fact, Fig. 3 shows that the phase speeds down to 25 m wavelength are consistent with linear wave theory but large fluctuations of the order of 1 m/s are found between 25 m and 20 m wavelengths, and these vary strongly with the choice of azimuth ϕ. Such fluctuations are not included in the surface current estimates made by [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF], because these authors exclude spectral components with a coherence under 0.8. This coherence, denoted "coh" in the following equations and figures, is also called magnitude-squared coherence. We note that this threshold is equivalent to a standard deviation of the co-spectrum phase ψ of 40 • , because for small values of the phase ψ in radians, std(ψ) 2 √ 1 -coh.

If the vertical shear in the top few meters is to be measured, we have to use these shorter wave components. Presumably we could use spectral components with a lower coherence, hence a larger uncertainty, and use the averaging over a larger number of spectral components to mitigate this larger uncertainty. For the shorter components, with k 40 cpkm, the coherence is under 0.35 for all directions, and highest for 110 < ϕ < 120 • , with a corresponding fluctuation of the phase std(ψ) 70 • . Interestingly, the same low coherence and high level of phase fluctuations are also present in the simulated data, even when the noise level is reduced to zero. We found that this pattern was not associated to the amplitude nor to the additive or multiplicative nature of the noise in eq. ( A1), nor even the non-linear modulation transfer function. The low coherence persists as long as some energy remains for waves in opposing directions. These fluctuations in the phase speed for the shortest wave components disappear in the simulation when the input spectrum is "chopped" to remove waves propagating from the east (with k x < 0, see Fig. 4). Clearly, the spurious large values of phase speeds for wavelengths 20 m < L < 25 m are associated to a significant level of energy in opposing directions.

Any spectral component (k, ϕ) contains information that propagates in both directions ϕ and ϕ + π. By interpreting the phase difference ψ 1,2 as the phase of a single travelling wave, in direction ϕ if the phase speed is positive, we are assuming that we can neglect the waves in the opposite direction. In fact, the data is in general the sum of two wave trains travelling in these opposite directions, each giving a different contribution to the phase difference ψ 1,2 , one of these two can be neglected if its energy is much weaker (typically with a difference of 20 dB or so), which is not the case in our example for L < 25 m. The magnitude of wave energy in opposing directions for wavelengths under 20 m has been particularly studied for the retrieval of surface wind direction using HF radar (e.g. Kirincich, 2016), but few studies have been performed

for longer wavelengths [START_REF] Tyler | Wave directional spectra from synthetic aperture observations of radio scatter[END_REF].

In order to quantify the magnitude of waves in opposing directions, we define an "opposition spectrum",

H(k, ϕ) = 4E(k, ϕ)E(k, ϕ + π) [E(k, ϕ) + E(k, ϕ + π)] 2 .
(5)

As defined, H ranges form 0 for waves propagating only in direction ϕ, to 1 for equal amplitudes in opposing directions ϕ and ϕ + π. This is the directionally-distributed counterpart of the "overlap integral" I(k) defined by [START_REF] Farrell | What do deep sea pressure fluctuations tell about short surface waves?[END_REF].

I(k)
was first used by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF] and [START_REF] Brekhovskikh | The radiation of infrasound into the atmosphere by surface waves in the ocean[END_REF] for the theory of generation of secondary microseisms and microbaroms (see [START_REF] Ardhuin | How ocean waves rock the earth: two mechanisms explain seismic noise with periods 3 to 300 s[END_REF]De Carlo et al., 2020, for recent reviews).

If H(k, ϕ) is independent of ϕ then I(k) = H(k, ϕ).
Starting from the same wave spectrum as in Fig. 4, we have simulated images deviation of the value of U , it is clear that a larger opposition spectrum gives a larger error. Part of the larger error comes from a larger number of spectral bins for which the standard deviation of the cross-spectral phase is larger than 60 • and are thus not included in the average. That effect also explains why no value is shown for N t = 0.2 and H > 0.41: all spectral bins in that case had a std(ψ) > 60 • . Even for N t = 0, the presence of opposing waves leads to very large biases on U that cannot be detected by inspecting only std(ψ). In other words, it may not be feasible to flag errors caused by the the presence of waves in opposing directions when using the co-spectrum phase to estimate the surface current. For the case shown in Fig. 5.a, the error can be up to 0.25 m/s for H = 0.2, which would be the value given right off the coast by a 5% coastal reflection, which is a typical value for steep beaches or rocky shores [START_REF] Ardhuin | Coastal wave reflection, directional spreading, and seismo-acoustic noise sources[END_REF]. These results are very robust and do not change qualitatively when changing the shape of the wave spectrum. We suspect that the general larger errors for smaller wavenumbers are associated to the smaller phase shift of the longer waves, corresponding to a lower signal to noise ratio.

Coastal reflections and longer wave components

We may look for further evidence for the effect of waves in opposing directions by is that the phase difference between the B04 and B02 images is biased low because some of the estimates are dominated by noise, even though our coherent sum of the co-spectra was weighted by the spectral energy. In this case the strongest spectral component in (Fig. 7.e) with k x = 2.5 cpkm and k y = -2.5 cpkm is the one with the largest velocity magnitude (17 m/s) in (Fig. 7.f), but the random distribution of phases gives an uncertainty of ±3 m/s, and the coherence is 0.94.

Looking at the full spectrum, we find that all the velocities are also probably biased by an error in the relative position (co-registration error) of the two bands B04 and B02, and/or an error in the time lag. Whereas the given geometry of the measurements gives a time lag of 1.00 s, it would take a roughly 1.0 m/s current in azimuths 20 • < ϕ < 80 • to explain the measured phase speeds for k in the range 30 to 40 cpkm, and the estimated current should vary like cos(ϕ -ϕ U ). Instead, the observed wave dispersion is more consistent with a time lag of 0.87 s and a much weaker current. That time lag difference of 0.13 s, with a phase speed of 7 m/s is also equivalent to a bias of 0.9 m in the location of the pixels, that could be caused by a bias of 1 microradians in the knowledge of the relative pointing of the different bands for the same detector. We also note that the 1 m/s order of magnitude of the possible error on the current velocity is consistent with the spurious stripes appearing in maps of surface current estimated by [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF] and O(1 m) co-registration errors found by [START_REF] Kääb | Glacier remote sensing using Sentinel-2. part I: Radiometric and geometric performance, and application to ice velocity[END_REF].

More interestingly for the purpose of the present paper, the variation of phase speed as a function of wavenumber has a O(50 cm/s) anomaly in box 2 for k around 20 cpkm. Could that be the signature of a current maximum at a depth around 10 m?

In the absence of verifying measurements we cannot explore this with any certainty.

3 Least squares method applied to a sequence of 3 images Going back to the problem of estimating phase speeds for the shorter wave components that often have relatively large values of H, we propose to try to separate the waves in opposing directions, and for this, use more than two images. This problem is very similar to the problem of separating waves in opposing directions in wave laboratory experiment, which is necessary for implementing absorbing boundary conditions at paddle wave makers. The founding paper in this line of work was the method of [START_REF] Mansard | The measurement of incident and reflected spectra using a least squares method[END_REF] for computing wave reflection using a series of 3 wave gauges with a least square method. It was later improved on by [START_REF] Zelt | Estimating incident and reflected wave fields using an arbitrary number of wave gauges[END_REF]. We are not aware of an adaptation to image processing and generalization of the method to estimate currents at the same time as the amplitudes of the two opposing wave trains. As detailed in the Appendix B, this generalization is straightforward.

We first test the method for simulated monochromatic waves of 50 m wavelength propagating in one dimension and resolved at dx=10 m resolution with time lags of 0.5 and 1.0 s similar to the red, green and blue bands (B04, B03, B02) of Sentinel-2.

We found that adding one extra measurement at a 0.8 s lag, similar to band B08 on Sentinel-2, had a limited impact on the results. The method is illustrated in Fig. 8 with spatial series, with or without noise. With a small value of the opposition spectrum, here A = 1, B = 0.1 and thus H 0.04, (except for B = 0.2 in Fig. 8.a), the wave field looks like a single propagating wave with a modulated amplitude, changing from 1.1 to 0.9 over half a wave period, here 2.8 s, due to the partial standing wave.

When multiplicative noise is added, the distribution of current estimates from phase differences is Gaussian for waves propagating in only one direction (Fig. 8.c). In the case of the least square method, the distribution has heavier tails than a Gaussian distribution and thus requires a very broad range of velocities to be properly characterized. In our example, the uncertainty is 20 times larger when using the least square method (see Fig. 8.b compared to 8.c). In contrast to A and B, we do not have an analytical expression for U and it is difficult to predict its distribution. In practice, we find that the median of the distribution of U is apparently not biased (Fig. 8.b,d,f), and the mean of M estimates of this median apparently converges following the central limit theorem with a standard deviation reduced by 1/ √ M . One may thus hope to retrieve the current with this method, even for noisy data.

We note that the uncertainty on the current U , as measured by the standard deviation of the distribution of U , is not affected by the presence of waves in both directions when using the least square method: changing the value of B from 0 (Fig. the least square method the current uncertainty generally decreases when reducing the noise (Fig. 8.d) or increasing the signal, for example when increasing the time differences between the measurements (Fig. 8.f). In contrast, when waves in opposing directions are present, the error reduction in the phase method is very limited when the noise is reduced (Fig. 8.e). Both methods are improved when the time lag is increased by a factor 10, as shown in Fig. 8.f,g.

We now apply the least-square method to actual Sentinel 2 imagery, with 3 bands (B02, B03, B04). We first note that the image amplitude and standard deviations are different for the different bands, so that the shift from one band to another is not just a propagation but also includes a change in mean value and amplitude. In order to mitigate that effect we have shifted and rescaled the pixel values so that each image has a zero mean and unit standard deviation before computing Fourier transforms.

We first take up our example off California, with results shown in Fig. 9. Noisy parts of the spectrum generally correspond to a low coherence in image pair (Fig. 9.a)

and a high uncertainty for the co-spectrum phase (already shown in Fig. 3) and hence current velocity. When using a least square fit, an obvious candidate for quantifying the noise is the residuals that we have normalized by the sum of the spectral densities of the images.

For each spectral component (k, ϕ) and each m-index sub-image of 100 by 100 pixels we fit the amplitudes Z A,m (k, ϕ) of a wave train travelling in direction ϕ and Z B,m (k, ϕ) travelling in direction ϕ + π and the current velocity U m (k, ϕ) that min-imizes the sum of square residuals that is the difference between the image spectral density B n,m (k, ϕ) and our model of two counter-propagating components (See appendix B). We only keep values of U that fall in the range from -5 to 5 m/s, for which there is a number M (k, ϕ) of estimates. For each spectral tile we have a normalized residual,

ε m (k, ϕ) = n |ε n (k, ϕ)| 2 / n |B n (k, ϕ)| 2 . ( 6
)
From these "successful fits", their number is M (k, ϕ), we take the current to be the median of the U m (k, ϕ) values and we define a root mean square residual, Inspecting (Fig. 9.b) we propose that a first not-too-conservative but reasonable threshold for acceptable results is ε r < 0.4, giving the current values shown in Fig. 9.d. This choice was motivated by the desire to include the spectral components for which we found that waves in opposite direction were a significant source of error for the phase method, but this also keeps spectral components with very low signals (with azimuth directions between 0 • and 30 • , and between 135 • and 180 • ). In addition to the current, we also estimate the opposition spectrum as

ε r (k, ϕ) = m ε 2 m /M (k, ϕ). ( 7
H(k, ϕ) = 4|Z A | 2 (k, ϕ)|Z B | 2 (k, ϕ) (|Z A | 2 (k, ϕ) + |Z B | 2 (k, ϕ)) 2 . ( 8
)
Although we have no direct measurement of the true directional spectrum E(k, ϕ) and thus of H(k, ϕ) using eq. ( 5), we may compare estimates H i to the values shown in Fig. 10.a and given by the wave spectrum estimate E(k, ϕ) using the Maximum Entropy Method, which is used in the image simulations. It is not clear at all if the MEM derived spectra give values of H(k, ϕ) that should have the same order of magnitude as those given true directional spectrum. This question could be investigated with stereo-video data (e.g., [START_REF] Guimarães | A data set of sea surface stereo images to resolve space-time wave fields[END_REF].

When using the least square estimate of Z A and Z B that correspond to the least square estimate of U for each image tile, the values recovered from the S2 image are typically much higher than those estimated from the buoy as shown in Fig. 10.a) : in the range 0.2 to 0.4 for the part of the spectrum that has a coherence squared higher than 0.64 which is already lower than the 0.8 threshold in [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF].

This value of H is probably at least 10 times too high around the spectral peak. These would typically give a background level of microseism sources that is too large by a factor 10 or more. What happens is that the fitting procedure puts some of the noise in the amplitude of the opposing waves. As a result, a lower threshold than ε r < 0.4 is necessary to give accurate estimates of the opposition spectrum H. But we can also force the current to the value estimated from the phase method and only fit Z A and Z B , in that case the values of H are more realistic, as shown in Fig. 10.b.

For that estimate we have also modified the equations in Appendix B to allow for a different current at times t 2 and t 3 in order to absorb the biases in the image position (δX, δY ) i,j = (U ij , V ij ) × (t j -t i ). Indeed the phase difference ψ 2,3 gives a velocity vector close to (-1.8,0.) while ψ 2,4 gives (-1.,0) corresponding to a 1 m eastward erroneous shift of the B02 image relative to B03 and B04. This inconsistency in the data is not included in the fitted model proposed in Appendix B and thus contributes to higher errors in the estimate of U . One possibility may be to recompute the least squares with different velocities over the different time lags, or to use the phase difference method on all image pairs to estimate deviations from a constant speed and shift the image before applying the least square method.

For our test image, it is thus dubious that the least-square method, as implemented here, has provided any additional reliable information for short waves compared to the phase method. Using a more conservative threshold ε r < 0.2 it is possibly able to slightly extend the part of the spectral plane from which a velocity can be derived to directions that are further away from the mean wave direction.

Looking beyond the particular case of the bands B02, B03 and B04 of the Sentinel 2 sensor, it is interesting to know how well this method may work. We have thus simulated the image and its processing, and reduced the noise level from N t = 0.15 (which looks similar to the true Sentinel 2 image) to no noise at all with N t = 0.

Without any noise, the least square fit is very good with ε r < 0.1 for the full spectral domain, except around the blind azimuth. As a result the input current vector U = (-1, 0) m/s is very well recovered. This would not be the case for the shortest components using the phase method except in the mean direction, giving only one component of the current vector.

The precision on the retrieval of the surface current is further illustrated in Figure 11, focusing on a narrow range of azimuths, between 110 and 120 • . The error bars give an estimate of the precision of the mean within each spectral bin that are all completely independent. For the phase-difference method, the smooth variation of the estimates across the spectra (within the error bar) confirm that the O(15 cm/s) precision for each spectra estimate is realistic. This does not say anything about the accuracy of the estimate that is dominated by an O(1 m/s) error due to relative pixel co-registration errors of the different bands.

For the least-square methods, the error bars are more difficult to define given the heavy tails of the U distribution and the sample size (256 independent spectra giving 256 estimates of U ). It might be possible to use the distribution of residuals ε m obtained for the M spectra as given by eq. ( 6), because they are correlated with errors on U , but we have not found a satisfactory parameterization that would work for both the academic 1D case of Figure 8 biases in the least square method may be dominated by artefacts of the processing method. We have not yet identified the source of these errors. We also note that the phase method, in contrast, has no trend in the simulated data for which the standard deviation of the phase is under 60 • .

We generally expect that errors can be reduced by increasing the time separation of the images so that the mean phase difference is much larger, making random phase differences comparatively smaller. Fig. 11.c,d shows that realistically noisy images with a doubled time lag are preferable to a noise-free image with the same time lag. This is easy to understand in the case of the phase difference method: the larger phase difference makes the random-phase noise a relatively smaller term in the phase difference. The uncertainty on U is inversely proportional to the time difference.

A first verification of this advantage of larger time lags is provided by using the B12 and B11 band, that are acquired 1.1 s and 0.5 s before B04, which is here 1 s before B02 (this ordering correspond to the even detectors on S2, it is reversed for the odd detectors). Hence combing B12 with B11 and B02, giving a maximum time lag of 2.1 s. Because the spatial resolution of B12 is only 20 m, we have averaged B02 over 2 by 2 pixel boxes to provide images at the same resolution, including a 1 m westward shift of B02 to corrected for the error noted above. These results are illustrated in Fig.

12.

We first note that the shape of the spectrum, here resolved at higher spectral resolution, shows a 3-lobe structure with minima of the image PSD and coherence for the azimuths 100 wavelength between 50 and 100 m. We note that the least square method gives rather puzzling results that we do not understand, with a variation of the estimated current as a function of wavenumber that is large and not random.

In the case of the waves shorter than 40 m wavelength, that are only resolved in the 10 m images such as given with bands B02 and B04 with 1 s time lag, the uncertainty of U from the phase difference method for wavenumbers from 30 to 40 cpkm is larger at 4.8 cm/s due to the opposing effects of a lower coherence and a larger number of spectral estimates.

4 Discussion: consequences for surface current velocity and shear re- From the consistency of the velocity estimates for all spectral components, and in the particular case of the image analysed in Fig. 2 and 11, we find that Sentinel 2 imagery is capable of providing a velocity precision of the order of 5 cm/s for spectral ranges of 10 cycles per kilometer. These uncertainties are of the order of the differences the first term of the matrix (A CA) -1 where A(m, 1) = k x,m , A(m, 2) = k y,m and C

is the diagonal matrix such that C(m, m) = 1/std(ψ m ) 2 [START_REF] Hogg | Data analysis recipes: Fitting a model to data[END_REF].

In our case, we take as a reference σ U,l = 0.018 m/s the uncertainty on retrieved U for the case of a linear MTF and linear surface. This uncertainty is not zero due to all the processing steps: use sub-images, windows, etc. Presumably the effect of non-linearities or waves in opposing directions are independent from these processing effects and we estimate the excess uncertainty due to effect Z, as

σ U,Z = σ 2 U -σ 2 U,l . (11)
The retrieved current and excess uncertainty are reported in Table 1 for 4 different simulations. We find that the effects of non-linearities introduced by the Choppy Wave Model for wavelengths larger than 20 m and in the case considered here, are comparable to the effect of the non-linear MTF and 9 times smaller than the effect of waves in opposite directions. Possibly nonlinearity becomes more important when waves much shorter than L = 20 m are taken into account.

Also, waves are not homogeneous in space, with gradient driven by the horizontal shear of small scale currents [START_REF] Ardhuin | Measuring ocean waves in sea ice using SAR imagery: A quasi-deterministic approach evaluated with Sentinel-1 and in situ data[END_REF][START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF][START_REF] Villas Bôas | Wave-current interactions at meso and submesoscales: Insights fromidealized numerical simulations[END_REF]. If the shorter waves correlate with currents in a way different from the longer waves, which can be the case at the smallest scales [START_REF] Suzuki | On the physical mechanisms of the two-way coupling between a surface wave field and a circulation consisting of a roll and streak[END_REF], what appears like a vertical shear in the difference of phase speed could be the effect of the horizontal shear. Detailed simulations of these effects will be needed to find the order of magnitude of horizontal shear contributions to the mean phase speed difference.

In general, the vertical shear of the current is a priori not sensitive to image co-registration errors because all wavelengths are affected by these errors in the same way, and the shear is associated by a difference in phase speed of the different wave components. We find that a 10 cm/s difference in phase speed between 50 m and 25 m wavelengths (k=20 cpkm and k=40 cpkm) should be detectable with Sentinel-2 using data from a 8 km by 8 km region of the ocean. Such a difference correspond to a fairly large current shear in the top 10 m of the ocean.

Beyond Sentinel 2: a possible STREAM-O instrument

Resolving weaker and more typical shears would require more sensitive measurements such as provided with larger time lags and higher spatial resolution. Fig. 13.b

shows that extending the spectrum to 100 cpkm (10 m wavelengths) would double the difference in velocity that can be detected. Such capability will be available with the next generation Sentinel 2. Using these shorter components will probably require methods that are less sensitive to the presence of waves in opposite directions, such as the least square method proposed here.

Based on the results shown in Section 3, increasing time lags is probably the most important feature needed to obtain more accurate currents using both the phase difference or the least square methods, with or without waves in opposing directions.

Due to data downlink constraints and cloud cover, optical monitoring of currents using wave dispersion cannot compete with the possible coverage given by radar systems [START_REF] Ardhuin | Satellite doppler observations for the motions of the oceans[END_REF]). Yet, optical imagery are uniquely capable of providing estimates of vertical current shear, and may provide very useful data for reducing pointing errors.

With these two goals in mind, the "ocean Surface Transport kinetic Energy, Air-sea fluxes and Mixing" (STREAM) concept combines a Ka-band Doppler scatterometer radar STREAM-R with a 1000 km wide radar swath, with a optical system STREAM-O that takes sequences of 5 ocean images at 5 m resolution with time lags of 1, 4, 9 and 10 s, overlapping over a 10 km wide strip along the orbit. This concept was proposed for the 11th Earth Explorer of the European Space Agency.

In addition to the increased time lag, The STREAM-O design uses a push-frame acquisition system, instead of the push-broom of the MSI sensor on Sentinel 2. Namely, the raw data for one band consists of 2-D images with large overlaps instead of a single line. Such a feature makes on-board averaging possible, increasing the integration time to a value larger than several times the life time of specular points, i.e. 10 milliseconds or more. This is clearly not feasible for a pushbroom system in which the duration of acquisition of each pixel is less than the pixel size (10 m) divided by the ground velocity (7 km/s), i.e. 1.4 ms. Future work is needed to understand the parameters that control the phase difference noise in optical imagery. It is possible that other parameters are relevant besides the number and life time of specular points. Also, the present paper considered images in the sun glint, but several useful applications (in particular for water depth retrieval) have been made with different observation geometries, with different noise characteristics.

Conclusions

In order to retrieve a surface current vector and current shear from observed wave dispersion it is necessary to obtain separate and robust estimates of the phase speed of different components of the wave spectrum, with different directions to obtain a current vector, and with different wavelengths that have different sensitivities to different depths.

Although the present work did not demonstrate a full solution method, we have highlighted difficulties associated to the retrieval of phase speed from a small number of ocean surface images using either a phase difference method or a least square fitting of the current velocity and the amplitude of waves in opposing directions. Both methods have complementary advantages and should probably be combined and modified for further improvements. We particularly highlighted how the presence of waves in opposite directions cause an error in the phase difference method. In one specific case analyzed here, this is particularly a problem for retrieving phase speeds from waves with wavelengths shorter than 4 times the dominant wind sea. The least square method using 3 or more images is not sensitive to waves in opposing directions, but it provides relatively noisy estimates of the current velocity when applied to Sentinel 2, due to the short time lags (about 1 s). As a result, the least square method may not provide much more useful additional information on the current velocity than the phase difference method. We also note that anomalously low coherence in image pairs may be an indication of the presence of waves in opposite directions, which may have application to the identification of strong microseism or microbarom sources.

Our simulations show that when applied to other sensors with lower image noise and/or larger time lags, the least square method may allow to use the shortest wave components that are more likely to be associated to high levels of energy propagating in opposing directions. We find that a 2 s time separation and the same pixel noise as Sentinel 2 should be sufficient to retrieve reliable phase speeds of shorter waves, all the way to the Nyquist wavelength. In that case it should be viable to reliably estimate the magnitude of waves in opposing directions as quantified by the opposition spectrum introduced in Section 2. Future work will be needed to refine and verify the error model for the two methods and their possible combination.

Figure 1 .

 1 Figure 1. Left: typical day and night velocity profiles of the total current in the Atlantic at 26N , 36W (adapted from Sutherland et al. 2016). Center: sensitivity kernels for surface gravity wave phase speeds. Right: depth of measurement of different instruments. From top to bottom:

Figure 2 .

 2 Figure 2. (a) Definition of viewing angles corresponding to a given sun and satellite sensor positions. The image brightness of a pixel is defined by the area of sub-pixel facets (in green) that gives a specular reflection and thus must have a given surface slope vector (sx, sy). That area is proportional to the probability density function within that pixel for the slope (sx, sy). This slope corresponds to the zenith angle β and azimuth ϕa. The perpendicular azimuths ϕ b = ϕa ± π/2 are "blind azimuths" in which the waves contribute a second order change to the pixel brightness and cannot be observed. (b) Position of satellite (S), observation point (O) and center of the Earth (C) in a vertical plane. (c) Triangle on the sphere joining the observation point O and the nadir positions N1 and N2 at observation times 1 and 2.

Figure 3 .

 3 Figure 3. Example of processing from Level-1C images to phase speeds, using 500 x 500 m tiles over a 8 by 8 km area. Top: data from Copernicus Sentinel 2 on 29 April 2016 off California (See Figs. 3-9 in Kudryavtsev et al. 2017), with β = 9 • , U10 = 6 m/s. Bottom: simulated Sentinel 2 data based on in situ wave spectrum determined from directional moments using the Maximum Entropy Method, and with random phases. The multiplicative noise amplitude is set to Nt = 0.15. The present paper was motivated by the phase speed anomalies, highlighted with the dashed magenta circle near the Nyquist wavelength L = 20 m.

Figure 4 .

 4 Figure 4. (a) Phase speeds for the simulated image in which either the full spectrum F (kx, ky) is taken or the right half of the spectrum F (kx ≥ 0, ky) is set to zero to have zero opposing wave energy, for directions 130 • ≤ ϕ ≤ 140 • . (b) Wave spectrum estimated from buoy data and used in the simulation, energy is represented in the direction from which waves are coming, i.e. corresponding to negative frequencies

Figure 5 .

 5 Figure 5. Errors in current retrieval for directions 110 • < ϕ < 120 • as a function of the opposition spectrum H using the same simulated spectrum as in Figure 4. (a) Mean current value (and thus error since the input current is zero) (b) standard deviation of the current.

Figure 6 .

 6 Figure 6. Sentinel 2 image off the North Shore of Oahu. (a) True color image and location of analysis boxes 1 and 2. (b) Wind speed and direction from ASCAT. (c) and (d) are the channel B04 values for Box 1 and Box 2 respectively.

Figure 7 .

 7 Figure 7. Analysis of 2 pieces of the S2 image shown in Fig. 6. For Box 1, (a) shows the PSD of the image intensity obtained with 2 km by 2 km tiles to give a better spectral resolution for wavelengths around 300 m, with dots marking the "low noise spectral components" that give a standard deviation of the co-spectrum phase under 60 • and (b) the phase speed of long wave components with directions between 130 • and 150 • , with blue bars corresponding to those low noise spectral components. (c) and (d) were obtained with 500 m by 500 m tiles focusing on shorter waves with directions between 45 • and 55 • . (e) to (h) show the same quantities for Box 2. The dashed lines in panels (a),(c),(e),(g) indicate the blind azimuth (see Fig. 2 for its definition) and the blue dots in (a) and (e) indicate the spectral components for which the std(phase) is less than 60 • , with velocity estimates shown with blue symbols in (b) and (f).

8

  .b) or 0.05 (Fig. 8.d) to 0.1 has no influence on the distribution (not shown). With

Figure 8 .

 8 Figure 8. (a) Example of 3 spatial series at 3 different times in the presence of rightward propagating and leftward propagating waves of wavelength L = 50 m, amplitudes A = 1 and B = 0.2 and a U = 0.2 m/s current. The time lags of 0.5 and 1 s are typical of Sentinel 2 data for the 10 m resolution R, G, and B channels. (b-g) Monte-Carlo simulations of the estimationof the current velocity U from 50000 monochromatic spatial series with white multiplicative noise of amplitude Nt, using least squares in (b,d,f), and using the phase difference between 2 images in (c,e,g). For reference the distribution of estimated currents is also shown in (b,c) when the amplitude of leftward propagating waves is zero.

Figure 9 .

 9 Figure 9. All the panels show quantities in spectral space (0 < kx < 50 cpkm, -50 < ky < 50 cpkm). Top panels (a,b,e,f,g): Error metric (coherence for phase method, standard deviation of normalized residual for least-squares method) Middle panel (c,d,h,i,j): Estimate of velocity in ϕ direction. The dashed circle corresponds to k = 50 cycles per km. The simulated images includes a non-zero current vector U = (-1, 0) m/s, so that the current component in direction ϕ should be -sin(ϕ), very close to what is retrieved in panel (i).

Figure 10 .

 10 Figure 10. Estimates of the opposition spectrum H(k, ϕ), for the spectral range (0 < kx < 50 cpkm, -50 < ky < 50 cpkm). Panel (a) shows the H spectrum based on the MEM-estimated spectrum from the in situ buoy, (b) and (c) show two estimates from the Sentinel 2 image using either the current estimated from phase differences or the current also given by the least squares method. (d) and (e) are the results for simulated images with different levels of noise.

  Fig. 11.a. Although some of these errors could be caused by instrument errors (such as errors in the retrieved observation angles that could change the estimate time lags and distort the dispersion relation), it is striking that the simulated data shown in Fig. 11.b gives similar errors, but slightly weaker, which leads us to think that the

Figure 11 .

 11 Figure 11. Comparison of different current estimates for waves in azimuths 110 • to 120 • for (a) Sentinel 2 data using bands B04, B03 and B02, and simulated data with the (b) same time lag and similar noise level, or (c) no noise, or (d) a doubled time lag. For the phase difference method (red and blue symbols) the error bars shows the mean value obtained for each spectral component plus or minus one standard deviation divided by the square root of the number of estimates. We have also tested (in green) using a sub-sample of the least-squares, keeping only those with small values of the residual εm.

  • and 125 • , these are probably due to artefacts of the Level 1-C processing. For the waves in the direction of highest coherence, 108 • < ϕ < 118 • , the uncertainty on U obtained in the range of wavenumbers 10 to 20 cpkm is as low as 0.1 m/s in spite of the average of only 64 independent tiles (compared to 256 for Fig.11.a). Combining all the 25 spectral components available from 10 to 20 cpkm gives an uncertainty of 3.4 cm/s, which we estimated from the mean of the uncertainties divided by the square root of the number of spectral components. Performing the same analysis on 20 m box averages of B03 and B02 gives a 5.8 cm/s uncertainty. It is therefore beneficial to use the largest time lags for estimating the current speed from

Figure 12 .

 12 Figure12. Example of results with a larger time lag of 2.1 s but coarser (dx=20m) using B12 and B02 bands. In order to better resolve the longer waves, the spectral analysis was done here with 1 km by 1 km tiles.

Table 1 .

 1 Effect of nonlinearities and opposing waves on retrieved current

	Simulation	U value σ u,Z
	Linear Model	-0.974	0.000
	Choppy Wave Model	-0.976	0.011
	nonlinear MTF	-0.975	0.004
	Opposing waves	-0.899	0.107

Acknowledgments

We acknowledge the use of Copernicus Sentinel 2 data, obtained from the Copernicus Science Hub https://scihub.copernicus.eu, and buoy data were furnished by the Coastal Data Information Program (CDIP), Integrative Oceanography Division, operated by the Scripps Institution of Oceanography, under the sponsorship of the U.S. Army Corps of Engineers and the California Department of Parks and Recreation https://doi.org/10.18437/C7WC72. F.A. and M.A. were supported by CNES as part of the SKIM preparation program and ANR grants for ISblue (ANR-17-EURE-0015) and MIMOSA (ANR-14-CE01-0012), with additional support from the European Space Agency, through the IASCO contract 4000129945/19/NL/FF/gp. M.Yu. was supported by Ministry of Science and Education of the Russian Federation under State Assignment No. 0555-2021-0004 and Russian Science Foundation through the Project No. 21-47-00038. We thank Bertrand Chapron for fruitful discussions and two anonymous reviewers that have helped make the manuscript much more readable. 10.3390/rs807059 Kirincich, A. (2016). Remote sensing of the surface wind field over the coastal ocean via direct calibration of HF radar backscatter power. J. Atmos. Ocean Technol., 33 (7), 1377-1392. doi: 10.1175//JTECH-D-15-0242.1 Kudryavtsev, V., Yurovskaya, M., Chapron, B., Collard, F., & Donlon, C. (2017a). Sun glitter imagery of surface waves. part 1: Directional spectrum retrieval and validation. J. Geophys. Res., 122 . doi: 10.1002/2016JC012425 -24-manuscript submitted to JGR: Oceans

in the advection speed of the different spectral components due to a typical vertical current shear in the top 20 m of the ocean. Along the equator with differences of the order of 50 cm/s between can be found between 1 m and 15 m depth and these should be detectable by Sentinel-2. A probably more typical shear, from subtropical Atlantic measurements, is shown in Fig. 13 and measuring it requires detecting 3 cm/s differences between k = 20 cpkm and k = 40 cpkm. Resolving this demands a reduction of the uncertainty by at least a factor 3, possibly obtained by averaging over at least 24 by 24 km. At that level of accuracy, other phenomena may be involved, such as non-linear wave effects, that include both an advection of short waves by long waves [START_REF] Weber | On the nonlinear theory for gravity waves on the ocean's surface. Part I: Derivations[END_REF][START_REF] Broche | Mesure par radar décamétrique cohérent des courants superficiels engendrés par le vent[END_REF] which corresponds to a fraction of the Stokes drift [START_REF] Ardhuin | Observation and estimation of Lagrangian, Stokes and Eulerian currents induced by wind and waves at the sea surface[END_REF], and the presence of harmonics (i.e. [START_REF] Leckler | Analysis and interpretation of frequency-wavenumber spectra of young wind waves[END_REF]. In order to test that latter effect we have made some simulations by replacing the linear wave model with the Choppy Wave Model [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF], which is a first order Lagragian theory. We use a current (U x = -1 m/s, U y = 0.) in our simulations.

Following [START_REF] Yurovskaya | Ocean surface current retrieval from space: The sentinel-2 multispectral capabilities[END_REF], we estimate U x and U y by a least square fit to the measured phase differences after removing the linear wave dispersion,

which is related to the current vector (U x , U y ) by

where ε is a small residual. We only keep the M spectral components with index m corresponding to one pair (k x , k y ) such that |k m | > 10 cpkm, |k m | < 40 cpkm and with std(ψ) < 60 • . This gives a length M vector Y = (ψ a,m ) 1≤m≤M . We invert eq.

(10) for all m to find the vector X = (U x , U y ) that gives the smallest sum of square residuals weighted by 1/std(ψ) 2 . The resulting uncertainty σ U is the square root of

Appendix A Image simulator

The first 5 directional moments are converted to a 5-degree resolution directional frequency spectrum using the Maximum Entropy Method [START_REF] Lygre | Maximum entropy estimation of the directional distribution in ocean wave spectra[END_REF].

This spectrum is then interpolated onto a regular grid in (k x , k y ) space to obtain power spectral densities of wave-induced surface elevation with a spectral resolution of 1/16000 cycles per meter, i.e. with a largest wavelength of 16 km, twice as large as the region analyzed. Drawing random phases for each spectral component, the wave power spectral density is used to define complex amplitudes that are inverse-Fourier transformed to generate 8 km square grids of the surface elevation and long wave slopes, (s x (x, y, t i ), s y (x, y, t i )), with x and y regularly discretized at 10 m resolution, and t i the discrete time sampling corresponding to the time of image acquisition.

The input to our image simulator are thus

• the wave spectrum F (f, θ) is transformed to F (k x , k y ) using linear dispersion with a maximum frequency f 0.7 Hz corresponding to a cut of wavelength of 3.2 m. This spectrum is then interpolated on a regular grid in (k x , k y ) space, covering the range [-k N , k N ] for each dimension where k N is the Nyquist frequency (i.e. corresponding to a 20 m wavelength when simulating 10 m resolution images). This avoids any aliasing.

• the direction of the dominant slopes ϕ mss (which is generally close to the wind direction)

• the mean square slope in that direction mss u and the mean square slope in the perpendicular cross-direction mss c .

• the bistatic view angles β and ϕ a for each image (see Figure 2 of specular points that contribute to the signal in each pixel [START_REF] Longuet-Higgins | Reflection and refraction at a random moving surface. II. number of specular points in a gaussian surface[END_REF].

The image pixel value is then taken as the nearest integer of a mean intensity I times (1 + n t )B/B 0 where n t is a random white noise of a amplitude N t that parameterizes the "twinkle" of the sea surface.

The noise of the detector is treated as an additive noise n d , represented as a Gaussian noise of standard deviation N d . For each channel j which corresponds to a time t j we have the pixel value

where the value E(x) is the largest integer value that is less or equal to x. The quantization effect of rounding to an integer pixel value is not very relevant in the present paper with examples that have a relatively bright sea surface. In contrast, the twinkle noise has a very important influence on the estimation of the surface current, as discussed in Sections 2 and 3.

Appendix B Adaptation of 3-probe least squares method to an unknown current

Let us have A and B the complex amplitudes of the waves propagating in the ϕ direction and the opposite direction ϕ + π, and N is the number of measurements at times t n . In the case N = 3, the system of equations for the 3 measured complex amplitudes F 1 , F 2 , F 3 at times t 1 = 0, t 2 , t 3 is, for each spectral component (k, ϕ),

Ae -i(σt2-kU t2) + Be +i(σt2+kU t2) -F 2 = ε 2 , (B2)

We look for the solution that minimizes the sum of the modulus of ε n squared,

Ae -i(σtn-kU tn) + Be +i(σtn+kU tn) -F n Ae iσ(tn-kU tn) + Be -i(σtn+kU tn) -F n (B5) where the overbar corresponds to the complex conjugate. Taking derivatives with respect to the real and imaginary parts of A and B and taking derivative with respect to U gives, respectively, n e -i(σtn-kU tn) Ae -i(σtn-kU tn) + Be +i(σtn+kU tn) -F n = 0 (B6) n e i(σtn+kU tn) Ae -i(σtn-kU tn) + Be +iσtn+kU tn -F n = 0 (B7) n t n Im Ae -i(σtn-kU tn) + Be i(σtn+kU tn) Ae -i(σtn-kU tn) + Be +i(σtn+kU tn) -F n = 0, (B8)

where Im(X) is the imaginary part of X.

Using t 1 = 0, this can be re-arranged as

where we have defined Finding the solution for f = 0 gives an estimate of the value of U . This operation can be repeated for each Fourier transform (each tile) and each spectral component.

Different averaging procedures are discussed in Section 3. In particular we find that the square root of the sum of |ε n | 2 is linearly correlated to the error on U , in particular when the phase differences are large. Finally, this approach is easily extended to more than 3 images.


Coline Monchanin

Mathieu Lihoreau

A toxic error

A toxic error

Heavy metal pollution limits shouldn't just keep humans safe, but other animals too, say Coline Monchanin and Mathieu Lihoreau ACROSS the world, insects are in decline. Intensive pesticide use, new diseases, habitat destruction and climate change are all contributing. Sadly, we are increasingly discovering that there is another impact of our everyday activities that is just as important: heavy metal pollution.

These substances are all around us. They are naturally present in Earth's crust and are released at low levels through weathering of rock and volcanic activity. But this gets a significant boost from human activity. We release these metals in various ways, ranging from dust that comes from vehicle brakes to the burning of fossil fuels for power and transport. All of this raises concentrations above natural levels. And once metallic dusts are out there, they stay for millennia.

While some of the compounds of these metals are essential for living organisms, most of them are highly toxic even at low concentrations. There are international guidelines designed to protect us from such pollution, but it turns out they aren't strong enough to do the same for insects.

In our recent work, we surveyed the scientific literature from the past 45 years that looked at the most monitored metals: arsenic, cadmium, mercury and lead. From this, we were able to identify the concentrations of these that are harmful to terrestrial invertebrates, the majority of which are insects. We then compared them with thresholds for the metals recommended by international regulatory bodies, for food, water and soil.

Though amounts of these in the environment should be below "human-safe" limits, in almost half of the studies the levels in natural conditions exceeded these figures. Not so surprisingly, at these concentrations the metals almost always killed invertebrates. More alarmingly, 90 per cent of the studies investigating metal levels within "humansafe" limits reported harmful effects on insects.

What's more, while it is clear that contamination harms these animals, we only have a partial picture as studies focus on sentinel species such as moths and bees, so-called because we use them to give an idea of possible harm to people, and on organisms favoured in many lab studies such as Drosophila. While we expect these to reflect a broad pattern, they only represent a tiny fraction of insect biodiversity, which is by far the most abundant of the terrestrial animals on the planet. A deeper understanding is critical to find the best possible solutions to mitigate these effects.

We need to start now, so we are calling for urgent action. While many local initiatives have worked to reduce emissions of some specific metal pollutants and emission sources over the past 40 years, environmental metallic pollution globally is still high. Take lead, for example, which started to be banned from petrol in the 1980s, but remains high in the environment because of the processing of ore and metals, and use of leaded aviation gasoline.

We could effectively reduce metal emissions now by dramatically cutting the use of fossil fuels, gasoline vehicles and metal-based pesticides, which are all sensible ways to protect our environments.

Politicians and scientists also have to rethink guidelines about what a safe level of metal pollution is, to take account of non-human species. This kind of shift in thinking was achieved a few years ago for broad scale use of harmful pesticides in intensive agriculture in many parts of the world.

Metal pollution is still an underappreciated threat. Its potential contribution to the huge insect biodiversity decline, called "insectageddon" by some, hasn't been recognised in recent conservation plans, such as the proposal for a more sustainable European Union Common Agricultural Policy. It is time for that to change.