Optimization of Sea Surface Current Retrieval Using a Maximum Cross-Correlation Technique on Modeled Sea Surface Temperature

Using sea surface temperature from satellite images to retrieve sea surface currents is not a new idea, but so far its operational near-real-time implementation has not been possible. Validation studies are too region specific or uncertain, sometimes because of the satellite images themselves. Moreover, the sensitivity of the most common retrieval method, the maximum cross correlation, to the parameters that have to be set is unknown. Using model outputs instead of satellite images, biases induced by this method are assessed here, for four different seas of western Europe, and the best of nine settings and eight temporal resolutions are determined. The regions with strong currents return the most accurate results when tracking a 20-km pattern between two images separated by 6–9 h. The regions with weak currents favor a smaller pattern and a shorter time interval, although their main problem is not inaccurate results but missing results: where the velocity is too low to be picked by the retrieval. The results are not impaired by the restrictions imposed by ocean surface current dynamics and available satellite technology, indicating that automated sea surface current retrieval from sea surface temperature images is feasible, for pollution confinement, search and rescue, and even for more energy-efficient and comfortable ship navigation.


Introduction
Knowledge about near-real-time ocean surface currents is essential to many commercial and societal sectors. These include, for example, life-saving search and rescue missions and oil spill confinement (Klemas 2012), but also navigation planning to reduce fuel consumption (Ronen 2011). Close to the coast, ocean surface currents can be measured by ground-based radars (Stewart and Joy 1974), but no such facilities exist in the open ocean. Furthermore, altimetry products are available only as daily 25-50-km grids (e.g., Pascual et al. 2007), which is not a high enough resolution for this type of applications.
Sea surface temperature (SST) in contrast is measured globally every 2-6 h by Advanced Very High Resolution Radiometer (AVHRR) instruments, at 1-km horizontal resolution, hence 30 years ago the idea of using SST to infer sea surface currents emerged. Emery et al. (1986) were the first to use the maximum cross-correlation (MCC) method for that purpose. Assuming that SST patterns are solely advected, this method tracks patterns that are most correlated between two consecutive satellite images and hence retrieves the current velocities by calculating the displacement of each pattern from one image to the next.
Since satellite images are affected by noise, clouds, geometric or geolocation mismatches, atmospheric distortion, etc., the validation of the MCC method is not straightforward and is better done with modeled sea surface current and temperature fields. Tokmakian et al. (1990) and Wahl and Simpson (1990) notably used models to assess the time interval that is needed between two images for the advection assumption to be valid. Wahl and Simpson (1990) found that horizontal diffusion was negligible for time intervals smaller than 24 h, and both studies concluded that a time interval smaller than 12 h is necessary for local heating or cooling effects (Wahl and Simpson 1990) and vertical mixing (Tokmakian et al. 1990) to be negligible compared to advection. These and subsequent studies (e.g., Emery et al. 1992;Simpson and Gobat 1994;Marcello et al. 2008) however used relatively simplistic ocean circulation models, so Doronzo et al. (2015) checked the validity of the MCC method using a high-resolution (1 km) 3-hourly realistic model, the Regional Oceanic Modeling System (ROMS; Shchepetkin and McWilliams 2005). But this and all the previous studies all reach the same conclusion: their results may be region specific, so they prove only that the MCC applied to SST is effective at retrieving sea surface current over the limited area that they studied.
In this paper, we use high-resolution ocean analysis products from the Copernicus Marine Environment Monitoring Service, available every hour at a 2-km horizontal resolution, to validate the use of the MCC method applied to sea surface temperature over four different regions and to show that the validation is in fact not region specific. To the best of our knowledge, we also provide the first sensitivity study of the settings of the MCC method. The model and methods that we use are briefly described in section 2. The results are presented and discussed in section 3, focusing first on the sensitivity study and then on the possible inaccuracies caused by the specificities of each sea. We conclude on the limitations and feasibility of an operational use of this method in section 4.

Data and methods
We use four weeks of hourly 0.038 horizontal resolution (2 km) output from the Atlantic-Iberian Biscay Irish-Ocean Physics Analysis and Forecast system provided by the EU Copernicus Marine Environment Monitoring Service. This model is based on an eddyresolving 0.038 application of the Nucleus for European Modelling of the Ocean (NEMO; Madec 2008), forced every 3 h by atmospheric fields from the European Centre for Medium-Range Weather Forecasts (ECMWF). For more information about this product and in particular its accuracy when compared with satellite imagery, the reader is referred to Sotillo et al. (2015). The outputs used here are the sea surface velocities u and y, and the SST, in January, April, July, and October 2016, from 0000 UTC on the 21st of the month to 2300 UTC on the 27th of the month. These were last accessed on 6 April 2017. We divided the output into the four seas that it spans (from north to south): To retrieve the sea surface velocities from the sea surface temperature we use the MCC method (Emery et al. 1986). It consists of tracking features in SST between two images (here two time steps). The user defines the size of the pattern from the first image to track, as well as the largest area to look around in the second image. The displacement of the water corresponds to the distance between the location of the original pattern and the location of the pattern in the second image, which is most correlated to the original one. Three settings hence have to be chosen to apply the MCC method, whose sensitivity is seldom studied: d the time interval between the two images (here considering only the images every 1 to 6, 9, and 12 h); d the size of the template in the first image (here 5 3 5, 10 3 10, or 20 3 20 km); d the size of the search area radius in the second image (here a radius of 10, 20, or 30 km around the template).
Hence in practice, the resolution of the retrieved velocity depends on the time interval and varies from 0.55 m s 21 for 1 h to 0.05 m s 21 for 12 h. To avoid any coastline effect, the algorithm as applied here does not consider the points within a search area of the coast and edges. The comparison between the reference velocities, directly provided by the model, and the retrieved velocities, obtained by applying the MCC method, is made using two main methods. We first assess which percentage of the retrieved field that is returned is within 0.2 m s 21 (as done by, e.g., Bowen et al. 2002;Warren et al. 2016) and 458 in direction of the original field. In a real-life situation, this would give an indication of the size of the work area and ensure that it is in the correct sector.
We also use Taylor diagrams (Taylor 2001), a common method when comparing climate model simulations. These diagrams provide a visual representation of which method is the closest to the reference (here the original model velocity fields), as explained in Fig. 1. The best setting is the one with the highest temporal correlation to the reference field, the smallest RMS difference to it, and the closest standard deviation to that of the reference. We plot one Taylor diagram per time interval setting and per sea, hence showing on the same diagram the nine settings corresponding to the different pattern sizes and search areas.

Results
In this study, we assess the performances of the MCC method in reproducing a known velocity field. A first requirement for validating this method is that a velocity is indeed returned. Table 1 shows the area, averaged over all the settings, where no value is returned, as well as the average surface current speed in the model. In fact, the four seas can be grouped in two categories: d the high-average-speed seas (North Sea and the Channel), which also have low areas where no value is returned (less than 10%); d and the low-average-speed seas (Bay of Biscay and western Mediterranean Sea), which have large areas where no value is returned by the MCC method (often over 20%).
For the four seas, there is no large difference between the four periods, although in both the Bay of Biscay and the Mediterranean Sea the MCC method returns the fewest values in January.
We shall now assess for the two categories how the accuracy of the MCC method varies with different settings, for the four periods.

a. The most accurate setting for each sea
Tables 2-5 present for each sea, from north to south, and for each period, the mean percentage of the study area for which a velocity value is returned by the MCC, where d the corresponding speed is within 0.2 m s 21 from the reference d the corresponding direction is within 458 from the reference.
Such criteria, in a real-life situation, would indicate in which sector of the water to search for a rescue operation even with no prior knowledge of the velocity field. For the North Sea (Table 2) and the Channel (Table 3), the individual settings corresponding to the maximum area with accurate velocities vary with the period and with the sea, yet common features are found. In particular the performances of the MCC method tend to TABLE 1. For the four regions and for the four periods considered, the average speed of the reference model's hourly output over the studied week (top) and the percentage of the region where no value is returned by the MCC method, averaged over all the time steps and the different settings (bottom).

North Sea
The  (Taylor  2001). The reference (here, the velocity from the model) is on the x axis, and the best setting is the one whose symbol lies closest to the reference. The distance to the reference depends mostly on the correlation between the reference field and the retrieved velocities (blue lines) and the standard deviation of both fields (the black quarter circle; the value on the x axis for the reference and y axis for each setting).
improve as the time interval between the SST ''images'' increases. For all periods, more than 50% of the area is accurate in the North Sea for a time interval larger than 2 h, and larger than 5 h in the Channel. More than 66% of the area is accurate over 4 h in the North Sea, and over 9 h in the Channel. In both seas, the largest values are encountered for a time resolution between 6 and 12 h, and for a template size of 20 km. For all four periods, in these strong current seas, the maximum area is large and exceeds 80% (values in bold font; Tables 2 and 3).
In the Bay of Biscay (Table 4) and western Mediterranean Sea (Table 5) in contrast, the area where MCC-calculated currents and reference currents coincide is rather low. In the Bay of Biscay, although the maximum in January and April is encountered for a template size of 20 km and a time interval of 2 h (Table 4), for all periods the accuracy is increased and larger than 50% for a time interval larger than 3 h and for a template size smaller than 20 km. Only January, which is also the month when the fewest values are returned (Table 1), has an accuracy larger than 66%, for a time interval between 3 and 5 h. The results are even poorer in the Mediterranean Sea where only in October, for a template size of 5 km and a time interval between 2 and 9 h do we have an accurate area larger than 50%. For both seas, as long as the time interval is larger than 1 h, there is no large seasonal difference. In summary, strong currents (in the North Sea and the Channel) are more accurately retrieved for a large template size (20 km) and time interval (9 h; Tables 2  and 3). This is consistent with Matthews and Emery (2009), who used a template size of 24 km to retrieve the mean strong California Coastal Current, and Bowen et al. (2002), who used 30 km for the east Australian coastal current. The large time interval probably allows for a wider range of velocities to be accurately retrieved. For the weak currents in contrast (in the Bay of Biscay and the western Mediterranean Sea), a lower template size (5 or 10 km) and time interval (3-6 h; Tables 4 and 5) return the most accurate results. This is consistent with Warren et al. (2016), who used a 10-km template size to retrieve the low velocities (0.12-0.25 m s 21 ) of the Korea Strait. Heuzé et al. (2017) in particular found that a small time interval was necessary in the Mediterranean Sea because of the presence of many eddies, which destroy the temperature patterns that the MCC method is designed to track. In all regions, the size of the search window did not seem to significantly impact the results, nor did the different periods.
This assessment was based on criteria that take into account both components of the velocity at the same time, but maybe, as a result of the geometry of each sea, one component is easier to retrieve than the other? This is the topic of the next section.

b. Zonal versus meridional component accuracy
Since the accuracy of the MCC method does not seem to be season dependent in this study, we now concentrate on only one of the periods. We choose 21-27 October 2016, which is not an extreme in velocity for any of the regions (Table 1). Figures 2 and 3 show the performance of each setting in the four seas, for the two velocity components u (gray color scale) and v (warm color scale) separately using Taylor diagrams, for four selected time intervals (2, 4, 6, and 9 h). The 12-h setting was discarded, since it is not performing better than 9 h in the previous section, and it has the lowest number of time steps (14). The focus in this section is not to determine which setting is the best, but whether one setting performs significantly better than the others, and whether there is a strong difference between u and y. With this type of representation, the most accurate setting is the one whose dot is closest to its reference (Taylor 2001).
For all seas, the range of standard deviation decreases as the time interval increases (Figs. 2 and 3, compared from top to bottom). This is probably because brief, extreme values are smoothed by the long-time averaging. The correlations are relatively similar for all time intervals of each sea though, suggesting that the model's surface velocity (the reference) and sea surface temperature have the same dynamics, at least at the temporal scales considered in this study. This would not be true for times larger than 12 h, that is, when horizontal diffusion, local cooling or heating, and vertical mixing are no longer negligible compared to advection (Tokmakian et al. 1990).
In the North Sea at low time intervals (Fig. 2a), the zonal component (u, gray) has too large a standard deviation and a lower correlation than the meridional component (y, warm) when compared to the reference. The standard deviation of u decreases as the time interval increases, whereas that of y does not change significantly. Similarly, in the Channel (Fig. 2, right), for low time intervals, the standard deviation of u is larger than twice that of the reference and hence its points are off the chart (Figs. 2e and 2f). Besides, for y all settings perform similarly, but for u, the three settings with a search area of 10 km (circle symbols) outperform the other settings. As in section 3a, we find that these two seas with strong currents have the most accurate retrievals for large time intervals.
Moreover in the Channel, a small search area leads to an increased accuracy of the retrieved zonal velocity. This is probably because the zonal SST gradient is very weak in the Channel over the period 21-27 October 2016 [0.028C (8longitude) 21 ], 5 times weaker than the meridional SST gradient, and an order of magnitude smaller than in the North Sea. As a result, the MCC algorithm can fail to distinguish SST patterns and return spurious high velocities, corresponding to unrealistic correlations between the east and west parts of the basin. A small search radius then constrains the velocities and limits their maximum value.
In the Bay of Biscay (Fig. 3, left), in agreement with Table 4, the accuracy of both u and v is best for the time intervals 4 and 6 h. No setting seems significantly better than the others, and the majority of settings have a correlation lower than in the North Sea and Channel. Similarly in the Mediterranean Sea (Fig. 3, right), the correlation rarely exceeds 0.7. As in Table 5, the 4-h time interval and 5-km template size (small black and magenta symbols) setting returns significantly better results than the others for u and y in the Mediterranean Sea.
In summary, apart from the Channel, which has a weak zonal SST gradient, there is no clear difference between the accuracy of the zonal and meridional components of the velocity. For the North Sea and the Channel (Fig. 2), no setting is particularly inaccurate. The accuracy is less in the Bay of Biscay and the Mediterranean Sea (Fig. 3). In fact, only in the Mediterranean Sea did one setting clearly perform better than the others, and this setting is not the same as was found using Table 5. So what caused the differences in Tables 2-5? We now check to see whether that could be the strength of the current itself.

c. Better if faster, stronger
Studies that aim at retrieving surface currents from SST usually concentrate on a small region and conclude that their result is probably region or current dependent (e.g., Emery et al. 1986;Tokmakian et al. 1990;Bowen et al. 2002;Doronzo et al. 2015). We verify this assumption here by looking at four different regions. Figure 4 shows for each region, of all the time steps of their respective setting with the largest number of accurate points returned (Tables 1-5), the one that is most accurate when comparing the original model velocities (left) to the ones retrieved (right).
In fact, the most striking feature is the amount of points where a velocity is not retrieved (white areas in Fig. 4). As was already shown by Table 1, in the North Sea and in the Channel (Figs. 4a and 4b, respectively), this number is relatively low and is mostly limited to the coastal regions that we purposely excluded from our algorithm. In the Bay of Biscay however (Fig. 4c)  even more in the western Mediterranean (Fig. 4d), there are many nonreturned points in areas with low velocities (dark blue). This is consistent with many studies, from Emery et al. (1986) to Heuzé et al. (2017), which showed that the MCC method cannot work if the temperature patterns do not move much during the time interval.
For all seas, both the velocity and the direction are correctly reproduced, even in the seemingly more complex areas. In the North Sea, for example, the current accurately turns westward north of 538N to follow the English coast (Fig. 4a). In the Channel, both the reference and retrieved currents are strongest around 28W and dip southward as they flow westward (Fig. 4b).
Even in the complex swirling Mediterranean Sea, when a direction is returned it is relatively accurate (Fig. 4d).
In fact, in our study, the accuracy is not location specific. For all the seas, the velocities returned are accurate. It is the number of these velocities that are returned that depends on the location. If the region has FIG. 4. A comparison between the (left) reference and (right) retrieved velocities for the (a)-(d) four regions of this study for their settings with the largest amount of accurate points returned. The same scale is used for all the panels. The number of wind arrows is decreased by four in latitude and longitude for clarity.
OCTOBER 2017 low velocities, as is the case in this study for the Bay of Biscay and the Mediterranean Sea (Figs. 4c and 4d), then a very large time interval between two images becomes necessary to retrieve more velocity value. But over 12 h, the assumption that the sea surface temperature is purely advected is not valid anymore (Tokmakian et al. 1990) and the MCC method can no longer be used.

Conclusions
In this paper, we compared the sea surface current velocities from a model with the current velocities retrieved by applying the maximum cross-correlation method to the sea surface temperature of the same model. The aims of this study were as follows: d to validate this cross-correlation method using a known velocity field from a model as a reference in order to see which biases are introduced by the algorithm; d to perform a sensitivity study in order to determine the best settings when applying this method, using the model output as a reference; d to compare the results over four different periods and four different seas in order to assess whether our results are region specific.
The biases induced by applying this algorithm-that is, the differences between the known model field and the retrieved field-are lower than 0.2 m s 21 in speed and 458 in direction for retrievals separated by only a few hours. We have thus shown that this algorithm is accurate enough and could theoretically be used operationally on actual spaceborne sea surface temperature data (e.g., to determine a search area for search and rescue operations; Davidson et al. 2009). The results in this study are not region specific, but in regions with low velocities the algorithm can less often retrieve values: the values retrieved are accurate, but there are fewer of them. We found that the highest accuracy was obtained for all regions for a time interval between two images of 4-9 h. This is compatible with already-existing satellite technology. In fact, Advanced Very High Resolution Radiometer images have an even higher horizontal resolution than the model used in this study (1 km), and since this instrument is aboard several satellites, images can be obtained every 2-6 h. Our study suggests that automated near-real-time retrieval of sea surface currents from AVHRR sea surface temperature is possible to implement.
One limitation though is that the presence of clouds can reduce the number of available images. The MCC method also assumes that the current is purely advected and would probably struggle in frontal zones and eddyrich areas. Yet, for the confinement of an oil spill, search and rescue missions, open ocean navigation fuel efficiency, or even the historical variability of mean currents, retrieving the surface currents from sea surface temperature seems to be an accurate and effective solution.