Ji Liu
email: ji.liu@inria.fr

Esther Pacitti
email: esther.pacitti@lirmm.fr

Patrick Valduriez
email: patrick.valduriez@inria.fr

Daniel De Oliveira
email: danielcmo@ic.uff.br

Marta Mattoso

Scientific Workflow Execution with Multiple Objectives in Multisite Clouds

INTRODUCTION

Large-scale in silico scientific experiments typically take advantage of Scientific Workflows (SWfs) to model data operations. A SWf is the assembly of scientific data processing activities with data dependencies among them. An activity is the description of a piece of work that forms a logical step within a SWf representation. A SWf Management System (SWfMS) is the tool to manage SWfs [START_REF] Liu | A survey of data-intensive scientific workflow management[END_REF]. In order to execute SWfs efficiently, SWfMSs typically exploit High Performance Computing (HPC) resources in a cluster, grid or cloud environment. Clouds have become an interesting solution for SWf execution because of various advantages. A cloud is typically made of several sites (or data centers), each with its own resources and data. Thus, in order to use more resources at different sites, SWfs could also be executed in a distributed manner at different sites. Nowadays, the computing resources or data of a cloud provider such as Amazon or Microsoft are distributed at different sites and should be used during the execution of SWfs. As a result, a multisite cloud is an appealing solution for large scale SWf execution. As defined in [START_REF] Liu | Scientific workflow partitioning in multi-site clouds[END_REF], a multisite cloud is a cloud with multiple sites, each at a different location (possibly in a different region) and being explicitly accessible to cloud users, typically in the data center close to them for performance reasons. In addition, there is a stored data constraint, i.e. some data cannot be moved to other sites because of big amounts or proprietary issues, in a multisite cloud.

To enable SWf execution in a multisite cloud, the execution of each activity should be scheduled to a corresponding site. Then, the scheduling problem is to decide where to execute the activities. In general, to map the execution of activities to distributed computing resources is an NP-hard problem. The scheduling problem may have multiple objectives, i.e. multi-objective, e.g. reducing execution time or monetary cost etc. Thus, the multisite scheduling problem must take into account the impact of resources distributed at different sites, e.g. different bandwidths, data distribution and costs to use Virtual Machines (VMs) at different sites. To the best of authors' knowledge, there is no solution to execute SWfs in a multisite cloud while considering both multiple objectives and dynamic VM provisioning. The related work either focuses on static VM provisioning, single objective or single site execution. Static VM provisioning refers to the use of the existing VMs (before execution) for SWf execution without modifying VMs during execution. In addition, existing cost models are not suitable for the SWfs that have a big part of sequential workload.

In this short paper (see [START_REF] Liu | Multi-objective scheduling of scientific workflows in multisite clouds[END_REF] for the extended version), we propose a general solution based on multi-objective scheduling to execute SWfs in a multisite cloud with the following main contributions: the design of a multi-objective cost model, SSVP VM provisioning approach, ActGreedy scheduling algorithm and an extensive experimental evaluation.

FRAGMENT SCHEDULING

In this section, we present the multisite SWfMS architecture, cost model, SSVP and ActGreedy.

The architecture of a multisite SWfMS is composed of four modules: workflow partitioner, multisite scheduler, single site initialization, and single site execution. The workflow partitioner partitions a SWf into fragments. After SWf partitioning, the fragments are scheduled to sites by the multisite scheduler. After scheduling, in order to avoid restarting VMs for the execution of continuous activities, all the activities scheduled at the same site are grouped as a fragment to be executed. Then, the single site initialization module prepares the execution environment for the fragment, using two components, i.e. VM provisioning and multisite data transfer. At each site, the VM provisioning component deploys and initializes VMs for the execution of SWfs. Finally, the single site execution module starts the execution of the fragments at each site. This can be realized by an existing single site SWfMS, e.g. Chiron [START_REF] Ogasawara | Chiron: a parallel engine for algebraic scientific workflows[END_REF].

It is difficult to estimate the execution time and monetary costs for the whole SWf even with a scheduling plan according to existing cost models [START_REF] Oliveira | A provenance-based adaptive scheduling heuristic for parallel scientific workflows in clouds[END_REF] since it is hard to generate a VM provisioning plan for each site with global desired execution time and monetary costs. As shown in Formula 2.0.1, we decompose the cost model as the sum of the cost of executing each fragment.

Cost(Sch(SW f, S)) = Schedule(wf i ,s j )=1 wf i ∈SW f Cost(wfi, sj) (2.0.1)
In the cost model, we also take the cost to provision VMs and the sequential workload into consideration in order to estimate the cost more precisely. In addition, we take advantage of Amdahl's law to estimate the execution time.

We propose a single site VM provisioning algorithm, called SSVP, to generate VM provisioning plans, which minimizes the cost to execute a fragment at Site s. First, SSVP calculates a near-optimal number of virtual CPUs to instantiate based on the cost model. Then, it optimizes the provisioning plan to reduce the cost to execute a fragment at Site s. Afterward, SSVP calculates the cost to execute the fragment at the site based on the cost model and improves the provisioning plan by inserting a new VM, modifying an existing VM or removing an existing VM. If the cost to execute the fragment at Site s can be reduced by improving the provisioning plan, the provisioning will be updated, and the improvement of provisioning plan continues.

Finally, we propose our multisite scheduling algorithms, i.e. ActGreedy, which schedules the most suitable site to each fragment. In ActGreedy, all the fragments of a SWf are not scheduled at beginning. During the scheduling process, all the fragments are scheduled at a corresponding site, which takes the minimum cost with the consideration of the stored data constraint while the cost is the minimum compared to other sites. As a result, the cost of executing a SWf in a multisite cloud is minimized. ActGreedy schedules fragments of multiple activities. ActGreedy can schedule a pipeline of activities to reduce data transfer between different fragments, i.e. the possible data transfer between different sites. A pipeline is a group of activities with a oneto-one, sequential relationship between them. In addition, ActGreedy makes a trade-off between time and monetary costs by using SSVP. ActGreedy schedules the available fragments, while choosing the best site for an available fragment rather than choosing the best fragment for an available site.

EXPERIMENTAL EVALUATION

In this section, we present an experimental evaluation of ActGreedy algorithm. All experiments are based on the execution of a SciEvol SWf in Microsoft Azure multisite cloud with three sites. We compare ActGreedy with LocBased [START_REF] Liu | Scientific workflow partitioning in multi-site clouds[END_REF] and SGreedy [START_REF] Oliveira | A provenance-based adaptive scheduling heuristic for parallel scientific workflows in clouds[END_REF], as well as with two general algorithms, i.e. Genetic and Brute-force. In the experiments, workflow partitioner, multisite scheduler and single site initialization are simulated, but the execution of fragments is performed in a real environment by Chiron.

The results show that LocBased corresponds to up to 21.75% higher cost than ActGreedy and SGreedy takes up to 74.51% higher cost than ActGreedy based on our proposed cost model. Based on an existing cost model, the advantage of ActGreedy can reduce the total cost up to 10.7% compared with LocBased and 17.2% compared with SGreedy. In addition, when the weight of reducing execution is low, ActGreedy may correspond to bigger execution time but when the weight is high, ActGreedy corresponds to less execution time. The execution with ActGreedy always takes less monetary cost compared with LocBased (up to 14.12%) and SGreedy (up to 17.28%). In addition, the experiments also show that the scheduling time of Genetic and Brute-force is much longer than ActGreedy (up to 577 times and 128 times). For instance, with more than 22 activities or 10 sites, the scheduling time of both Genetic and Brute-force exceeds the execution while the scheduling time of ActGreedy remains small.

CONCLUSION

In this paper, we proposed a general solution based on multi-objective scheduling to execute SWfs in a multisite cloud (from the same provider). We first proposed a novel multi-objective cost model, based on which, we proposed a dynamic VM provisioning approach, namely SSVP, to generate VM provisioning plans for fragment execution. The cost model aims at minimizing two costs: execution time and monetary costs. Our proposed fragment scheduling approach that is ActGreedy, allows for considering stored data constraints while reducing the cost based on our multiobjective cost model to execute a SWf in a multisite cloud. We evaluated our approaches by executing SciEvol in Microsoft Azure cloud. The results show that since it makes a good trade-off between execution time and monetary costs based on SSVP, ActGreedy leads to the least total normalized cost, which is calculated based on our multi-objective cost model, than LocBased (up to 21.75%) and SGreedy (up to 74.51%) approaches. Additionally, ActGreedy scales very well, i.e. it takes a very small time to generate the optimal or near optimal scheduling plans when the number of activities or sites increases, compared with general approaches, e.g. Genetic and Brute-force.

(c) 2016 ,

 2016 Copyright is with the authors. Published in the Proceedings of the BDA 2016 Conference (15-18 November, 2016, Poitiers, France). Distribution of this paper is permitted under the terms of the Creative Commons license CC-by-nc-nd 4.0. (c) 2016, Droits restant aux auteurs. Publié dans les actes de la conférence BDA 2016 (15 au 18 Novembre 2016, Poitiers, France). Redistribution de cet article autorisée selon les termes de la licence Creative Commons CCby-nc-nd 4.0. BDA 2016, 15 au 18 Novembre 2016, Poitiers, France.


Carlyna Bondiombouy

Boyan Kolev

Patrick Valduriez

Oleksandra Levchenko

Extending CloudMdsQL with MFR for Big Data Integration *

INTRODUCTION

Multistore systems [START_REF] Bondiombouy | Query processing in multistore systems: an overview[END_REF] have been recently proposed to provide integrated access to multiple, heterogeneous data stores through a single query engine. Compared to multidatabase systems [START_REF] Özsu | Principles of Distributed Database Systems[END_REF], multistore systems typically trade source autonomy for efficiency, using a tightly-coupled approach. In particular, much attention is being paid on the integration of unstructured big data (e.g. produced by web applications) typically stored in HDFS with relational data, e.g. in a data warehouse.

Existing solutions to integrate such unstructured and structured data do not directly apply to solve our problem, as they rely on having a relational view of the unstructured data, and hence require complex transformations. SQL engines, such as Hive, on top of distributed data processing frameworks are not always capable of querying unstructured HDFS data, thereby forcing the user to query the data by defining map/reduce functions.

Our approach is different as we propose a query language that can directly express subqueries that can take full advantage of the functionality of the underlying data processing frameworks. Furthermore, the language should allow for query optimization, so that the query operator execution sequence specified by the user may be reordered by taking into account the properties of map/filter/reduce operators together with the properties of relational operators. We respect the autonomy of the data stores and pay special attention to this throughout our experimental evaluation.

In this short paper (see [START_REF] Bondiombouy | Multistore big data integration with cloudmdsql[END_REF] for the long version), we propose a functional SQL-like query language (based on Cloud-MdsQL) and query engine to retrieve data from two different kinds of data stores -an RDBMS and a distributed data processing framework such as Apache Spark or Hadoop MapReduce on top of HDFS -and combine them by applying data integration operators (mostly joins). However, users need to be aware of how data are organized across * Work partially funded by the European Commission under the Integrated Project CoherentPaaS the data stores, so that they write valid queries. The query therefore contains embedded invocations to the underlying data stores, expressed as subqueries. As our query language is functional, it introduces a tight coupling between data and functions.

QUERY LANGUAGE

The query language is based on a more general common query language, called CloudMdsQL [START_REF] Kolev | The cloudmdsql multistore system[END_REF], designed in the context of the CoherentPaaS project [1] to solve the problem of querying multiple heterogeneous databases (e.g. relational and NoSQL) within a single query while preserving the expressivity of their local query mechanisms. The common language itself is SQL-based with the extended capabilities for embedding subqueries expressed in terms of each data store's native query interface. The common data model respectively is table-based, with support of rich datatypes that can capture a wide range of the underlying data stores' datatypes, such as MongoDB arrays and JSON objects, in order to handle non-flat and nested data, with basic operators over such composite datatypes.

In this section, we introduce a formal notation to define Map/Filter/Reduce (MFR) subqueries in CloudMdsQL that request data processing in an underlying big data processing framework (DPF). Then we give an overview of how MFR statements are combined with SQL statements to express integration queries against a relational database and a DPF.

MFR Notation

An MFR statement represents a sequence of MFR operations on datasets. A dataset is considered simply as an abstraction for a set of tuples, where a tuple is a list of values, each of which can be a scalar value or another tuple. Although tuples can generally have any number of elements, mostly datasets that consist of key-value tuples are being processed by MFR operations. In terms of Apache Spark, a dataset corresponds to an RDD (Resilient Distributed Dataset -the basic programming unit of Spark). Each of the three major MFR operations (MAP, FILTER and REDUCE) takes as input a dataset and produces another dataset by performing the corresponding transformation. Therefore, for each operation there should be specified the transformation that needs to be applied on tuples from the input dataset to produce the output tuples. Normally, a transformation is expressed with an SQL-like expression that involves special variables; however, more specific transformations may be defined through the use of lambda functions.

Queries that integrate data from both a relational data store and a DPF usually consist of two subqueries (one expressed in SQL that addresses the relational database and another expressed in MFR that addresses the DPF) and an integration SELECT statement. The syntax follows the CloudMdsQL grammar introduced in [START_REF] Kolev | Cloudmdsql: querying heterogeneous cloud data stores with a common language[END_REF]. A subquery is defined as a named table expression, i.e. an expression that returns a table and has a name and signature. The signature defines the names and types of the columns of the returned relation. Thus, each query, although agnostic to the underlying data store' schemas, is executed in the context of an ad-hoc schema, formed by all named table expressions within the query. A named table expression can be defined by means of either an SQL SELECT statement (that the query compiler is able to analyze and possibly rewrite) or a native expression (that the query compiler considers as a black box and passes to the wrapper as is, thus delegating it the processing of the subquery).

EXPERIMENTAL VALIDATION

The goal of our experimental validation is to evaluate the impact of query rewriting and optimization on execution time. More specifically, we explore the performance benefit of using bind join under different conditions.

Datasets

We generated data to populate the PostgreSQL table scientists, the MongoDB document collection publications, and text files with unstructured log data stored in HDFS.

Experimental Results

To evaluate the impact of optimization on query execution, we use a cluster of the GRID5000 platform (www.grid 5000.fr), with one node for PostgreSQL and MongoDB and 4 to 16 nodes for the HDFS cluster. The Spark cluster, used as both the DPF and the query processor, is collocated with the HDFS cluster. Each node in the cluster runs on 16 CPU cores at 2.4GHz, 64GB main memory, and the network bandwidth is 10Gbps. To demonstrate in detail the optimization techniques and their impact on the query execution, we prepared 1 query with different selectivity factors of the bind join condition. We execute the query in three different HDFS cluster setups -with 4, 8, and 16 nodes. Then we compare the execution times without and with bind join to the MFR subquery, which are illustrated in each query's corresponding graphical chart.

Query 1 involves all the data stores and aims at finding experts for publications of authors with a certain affiliation, it uses the MFR subquery experts_alt , which uses more sophisticated map functions, but makes only one shuffle, where the key is a keyword. For comparison, the MFR expression experts makes two shuffles, of which the first one uses a bigger key, composed of a keyword-author pair. Therefore, the corresponding Spark computation of Query 1 involves much smaller size of data to be shuffled, which explains its better overall efficiency and higher relative benefit of using bind join.

- This experimental evaluation illustrates the query engine's ability to perform optimization and choose the most efficient execution plan. The results show the significant benefit of performing bind join in our experimental scenario, despite the overhead it produces (see [START_REF] Bondiombouy | Multistore big data integration with cloudmdsql[END_REF]).

CONCLUSION

In this short paper, we proposed a functional SQL-like query language and query engine to integrate data from relational, NoSQL, and big data stores (such as HDFS). Our validation demonstrates that the proposed query language achieves the following requirements. First, it provides high expressivity by allowing the ad-hoc usage of specific map/filter/reduce operators through the MFR notation, as it was demonstrated with the hdfs subqueries. Second, it is optimizable as was demonstrated through performing bind join by rewriting the MFR subquery after retrieving the dataset from the MongoDB database. Our performance evaluation illustrates the query engine's ability to optimize a query and choose the most efficient execution strategy.

  [1].

(c) 2016 ,

 2016 Copyright is with the authors. Published in the Proceedings of the BDA 2016 Conference (15-18 November, 2016, Poitiers, France). Distribution of this paper is permitted under the terms of the Creative Commons license CC-by-nc-nd 4.0. (c) 2016, Droits restant aux auteurs. Publié dans les actes de la conférence BDA 2016 (15 au 18 Novembre 2016, Poitiers, France). Redistribution de cet article autorisée selon les termes de la licence Creative Commons CCby-nc-nd 4.0. BDA 2016, 15 au 18 Novembre, Poitiers, France.

  AND p.author = s.name AND e.kw IN p.keywords -Query 1.2: selectivity factor 30% SELECT p.author, p.title, e.kw, e.expert FROM scientists s, publications p, experts_alt e WHERE s.affiliation IN ('affiliation1', 'affiliation2', 'affiliation3') AND p.author = s.name AND e.kw IN p.keywords

	Query 1.1: selectivity factor 10%
	SELECT p.author, p.title, e.kw, e.expert
	FROM scientists s, publications p, experts_alt e
	WHERE s.affiliation = 'affiliation1'


, where several inferred sources fell on the Antarctic continent.

Using spherical ray tracing we investigate the hypothesis that the refraction of waves by surface currents produces the mirages. With reconstructions of surface currents inferred from satellite altimetry, we show that mesoscale vorticity significantly deflects swell away from great-circle propagation so that the source and receiver are connected by a bundle of many rays, none of which precisely follow a great circle. The ±10 • directional fluctuations at the receiver result from the arrival of wave packets that have travelled along the different rays within this multipath. The occasional failure of great-circle backtracking, and the associated mirages, probably results from partial topographic obstruction of the multipath, which biases the directional average at the receiver.

Introduction

Following the disruption of the 1942 Anglo-American landings on the Atlantic beaches of North Africa by six-foot surf [START_REF] Atkinson | An Army at Dawn -the War in North Africa[END_REF], the forecasting of surface gravity waves became a wartime priority. These first surface-wave forecasts were based on weather maps, and the relatively predictable propagation of swell, and were used to determine optimal conditions for amphibious assault [START_REF] Von Storch | Seventy Years of Exploration in Oceanography: a prolonged weekend discussion with Walter Munk[END_REF][START_REF] Bates | Utilization of wave forecasting in the invasions of Normandy, Burma, and Japan[END_REF]. Wartime work showed that long surface waves, generated in stormy regions of the globe, can travel for many thousands of kilometers before breaking on distant shores [START_REF] Barber | The generation and propagation of ocean waves and swell. I. wave periods and velocities[END_REF]. In the following decades, Munk, Snodgrass and collaborators observed that swell generated in the ocean surrounding Antarctica travels half way around the Earth [START_REF] Munk | Measurement of southern swell at Guadalupe Island[END_REF][START_REF] Munk | Directional recording of swell from distant storms[END_REF][START_REF] Snodgrass | Propagation of ocean swell across the pacific[END_REF]. Thus, after a transit of 5 to 15 days, the waves created by winter storms in the Southern Ocean produce summer surf in California. [START_REF] Barber | The generation and propagation of ocean waves and swell. I. wave periods and velocities[END_REF] used linear wave theory to relate the range of a distant storm to the rate of change of peak wave frequency at an observation point. But the formula of Barber & Ursell -see ( 2) below -provides no information about the direction of the source. The first attempt at measuring the direction of the incoming swell was made by [START_REF] Munk | Directional recording of swell from distant storms[END_REF] using an array of three pressure-transducers on the sea bottom offshore of San Clemente Island. The method is analogous to astronomical interferometry and was used to infer the direction of incoming wave-trains from June to October of 1959. Assuming that swell travels on great-circle routes, these observers combined Barber & Ursell's estimate of the range with interferometric direction to locate storms in the Southern Ocean. Although these inferred sources could be related to storms on weather maps, the location errors for some events were as much as 10 • of arc along the surface of the Earth, or 1000km. Moreover, there is a systematic error in the observations of [START_REF] Munk | Directional recording of swell from distant storms[END_REF]: the inferred source is most often to the south of the actual storm (Munk, 2013b). Because of this southwards shift, three of the thirty inferred storms even fell on the Antarctic continent, and several others on sea-ice. The unavoidable conclusion is that surface waves do not travel precisely on great-circle routes. Effects such as planetary rotation [START_REF] Backus | The effect of the earth's rotation on the propagation of ocean waves over long distances[END_REF], or the slightly spheroidal figure of the Earth, were shown to be far too small to explain the observed departure from great-circle propagation. We now believe that these mirages were caused by a combination of topographic obstruction [START_REF] Munk | Corrigendum to "Directional recording of swell from distant storms[END_REF] and refraction by surface currents. [START_REF] Kenyon | Wave refraction in ocean currents[END_REF] advanced the hypothesis that the refraction of surface gravity waves by major ocean currents, particularly the Antarctic Circumpolar Current, might explain the discrepancies reported by [START_REF] Munk | Directional recording of swell from distant storms[END_REF]. A packet of surface waves propagating through these currents is in a medium that varies slowly on the scale of a wavelength. The currents induce a doppler-shift in the surface wave frequency,

Ω(x, k) = gk + u(x) • k , (1)
with k the wavevector, k = |k|, u the surface velocity, and g the acceleration of gravity. Over many wavelengths refraction by the currents changes the direction of propagation of a wave packet so that it departs from a great circle. Given the surface current u(x, t), the ray equations describe the coupled evolution of the local wavevector k(t) and the position of the wave-packet x(t), just as the propagation of light in a slowly varying medium fol-lows the rules of geometrical optics [START_REF] Whitham | A note on group velocity[END_REF][START_REF] Bühler | Waves and Mean Flows[END_REF][START_REF] Landau | of a Course of Theoretical Physics[END_REF]. Using an idealized model of the Antarctic Circumpolar Current as a parallel shear flow, [START_REF] Kenyon | Wave refraction in ocean currents[END_REF] solved the ray-tracing equations and found appreciable deflections for rays making a grazing incidence with the current.

Here we re-visit Kenyon's hypothesis with the benefit of modern reconstructions of ocean surface currents and wave climate. In section 2 we analyze modern swell measurements from pitch-and-roll buoys [START_REF] Kuik | A method for the routine analysis of pitch-and-roll buoy wave data[END_REF]. We observe that the direction of incident swell measured by a deep-water buoy close to San Clemente Island varies by ±10 • on a time scale of hours. After averaging over these directional fluctuations to determine the mean direction, we use great-circle backtracking and find that for most events there is excellent agreement between the location of the inferred source and wave maxima on weather maps. But, in a few cases, great-circle backtracking places the wave source on land, as much as 10 • from the nearest storm on weather maps. Section 3 develops the theory of spherical ray tracing, including the effects of surface currents. In section 4 we compute ray paths through realistic surface currents and show that because of deflection by currents the source and the receiver are connected by a bundle of rays (a "multipath"). None of the rays in a multipath follow the great circle. In section 5 we show that refraction by surface currents can explain both the magnitude of the directional fluctuations at the receiver and their quantitative dependence on wave frequency. We conclude in section 6 by proposing a mechanism based on the interplay between surface currents and topography to explain the systematic shift towards the South and the frequent inference of sources on land in the observations by [START_REF] Munk | Directional recording of swell from distant storms[END_REF].

Modern data

We have reproduced the observations of [START_REF] Barber | The generation and propagation of ocean waves and swell. I. wave periods and velocities[END_REF] and [START_REF] Munk | Directional recording of swell from distant storms[END_REF] using modern measurements provided by recordings of waves from pitch-and-roll buoys deployed by NOAA's National Data Buoy Center (NDBC). Using three-dimensional accelerometers, these buoys measure the wave height, period and direction [START_REF] Kuik | A method for the routine analysis of pitch-and-roll buoy wave data[END_REF]. The intensity of the swell can be characterized by the Significant Wave Height (SWH), defined as the average height of the highest one-third of the waves, and recorded by the station every hour during a 20-minute sampling period. The SWH is usually of the order of four times the root-mean-square surface elevation.

NOAA Station 46086 is located in the San Clemente basin at latitude 32.5 • N and longitude 118.0 • W, close to where [START_REF] Munk | Directional recording of swell from distant storms[END_REF] used bottom-pressure transducers to measure swell direction. The accelerometer buoy at station 46086 has the great advantage of floating in 2000 meters of water, so that the influence of bathymetry is negligible. By contrast, the station of [START_REF] Munk | Directional recording of swell from distant storms[END_REF] was in approximately 100 meters of water, and refraction by local bathymetry impacted the measured swell direction (Munk, 2013a).

The spectral intensity is provided in frequency bins spaced by 5mHz in the lowfrequency range, together with the mean direction of the signal for each one of these bins. These directional spectra are based on an observation length of twenty minutes and are delivered every hour. The arrival of swell from a distinct source at buoy 46086 is indicated by a strong spectral peak that shifts towards higher frequencies as time goes on, with a timescale of several days: see Fig. 1. The progressive shift of the spectral peak to higher frequencies is caused by the dispersive propagation of surface waves: longer waves travel faster, so that fast low-frequency swell reaches the receiver first [START_REF] Barber | The generation and propagation of ocean waves and swell. I. wave periods and velocities[END_REF].

To isolate swell events, we track the peak frequency of the low-frequency part of the spectrum. To smooth the discretization from frequency bins, we first locate the lowfrequency maximum in spectral intensity, before performing a weighted average of this frequency with the frequencies of the two neighbouring bins, using weights equal to the spectral intensity in each of the three bins.

From each directional spectrum, we extract the incident angle of the swell, which we define as the mean direction for the frequency bin with maximum spectral intensity: this incident angle is computed every hour using 20-minute long samples.

a. A case study: the event of November 2007

The left-hand side panels of Fig. 2 focus on the particularly striking event in November 2007: there are four consecutive days of swell with SWH at around 1.2 meters in panel (e). Fig. 2(a) shows the wave propagation diagram. As first observed by Barber & Ursell, the peak frequency increases linearly with time and, using ω = √ gk, one can infer the distance L between the source (the storm) and the receiver (buoy 46086) from the slope of the frequency versus time on the wave propagation diagram:

df dt = g 4πL . (2)
Above, f = ω/(2π) is the frequency in Hertz. Using a linear fit to the early part of the swell event reported in figure 2, we find a distance of 86.7 • of arc between the storm and buoy 46086. The intersection of the fitting line with the f = 0 axis gives the date of birth of the storm which is October 29th 2007, at around 23:00 GMT.

To infer the direction of the source we turn to the incident-angle signal in Fig. 2(c). There are ±10 • fluctuations in the direction of the incident waves at the buoy. We remove these fluctuations by averaging, and so find that the wave signal at station 46086 comes from 204 • , measured clockwise from North. We hypothesize that the ±10 • directional fluctuations in Fig. 2(c) are too large to be simply instrumental noise, and that there might be physical information in the directional measurements. We return to further discussion of this point in section 5 e.g., see Fig. 11 and the supporting discussion.

The inferred range and direction locate the wave source on weather maps that are made available by the European Center for Medium-Range Weather Forecasts (ECMWF). The ECMWF interim reanalysis [START_REF] Dee | The ERA-Interim reanalysis: Configuration and performance of the data assimilation system[END_REF] provides sea-ice cover, 10-meter wind speed, and also SWH computed using the WAM wave model and assimilation of altimeter data [START_REF] Hasselmann | The WAM model -a third generation ocean wave prediction model[END_REF][START_REF] Komen | Dynamics and Modelling of Ocean Waves[END_REF]. We identify the southern sources of swell, corresponding to southern-ocean storms, as strong local maxima in 10-meter wind speed, and as large SWH. In Fig. 3 we compare the relevant ECMWF SWH with the source inferred from the accelerometer buoy recordings (the red dot). The buoy data analyzed above predict a storm at a range of 86.7 • on a great-circle route making an angle of 204 • going clockwise from North at the buoy, on October 29th 2007 at 11 pm. The inferred location, shown as a red dot in Fig. 3, is in excellent correspondence with the SWH maximum. In the example of Fig. 3, great-circle backtracking works very well: refraction by surface currents and topographic effects do not spoil the inference of the

b. Another case study: the event of April 2007

The agreement shown in Fig. 3 is the most frequent situation we found in our analysis of the data from accelerometer buoys: for 14 of the 18 swell signals we analyzed, the inferred source corresponds to a maximum in ECMWF surface wave height within 5 • of arc. But we also found a few examples for which the swell signal was very clean e.g., Fig. 1(b), yet great-circle backtracking resulted in a bad estimate of source location. In the right-hand panels of Fig. 2, and in Fig. 4 we focus on this event in April 2007. The SWH at station 46086 is of the order of 1 meter and the frequency of the swell increases linearly with time. Yet in Fig. 4 the direction of the inferred source is approximately 10 • from the maximum in ECMWF surface wave height. Moreover, the 10 • error puts the inferred source of this April 2007 event on New Zealand. This mirage is reminiscent of the Antarctic sources inferred by [START_REF] Munk | Directional recording of swell from distant storms[END_REF].

Mirages seem to occur preferentially when there is shallow topography, or even land, close to the great-circle route between the storm and the receiver. In the case of Fig. 4, a dense part of the Tuamotu Archipelago between 14 • to 18 • South and 148 • to 140 • West -see Fig. 50 of [START_REF] Munk | Directional recording of swell from distant storms[END_REF] -blocks the wave packets propagating on the great-circle route between the storm and the receiver. Only rays that are deflected strongly enough by surface currents to go around this Tuamotu blockage can reach the receiver. In anticipation of results from section 4, a possibility is that because of the distribution of surface currents, some wave-packets were deflected west of the Tuamotu blockage, so that the inferred source appears on New Zealand: the mirage in Fig. 4 results from the interplay between Tuamotu blockage and refraction by surface currents.

Waves on the surface of a sphere

Modeling the results of section 2 requires tracing rays on the surface of a sphere, including the effects of refraction by surface currents. [START_REF] Backus | The effect of the earth's rotation on the propagation of ocean waves over long distances[END_REF] developed this ray theory, including shallow-water effects and rotation, but without considering ocean currents. In this section we provide an account of the relevant theory required for the model in section 4.

a. The ray equations in spherical coordinates

Let us denote the phase of a wavepacket by S(x, t). Then the frequency and local wavevector are respectively -∂ t S and k = ∇S. The dispersion relation can be written in terms of S(x, t) as ∂ t S + Ω(x, ∇S) = 0. This partial differential equation is the Hamilton-Jacobi equation for a mechanical system with action S and Hamiltonian Ω. The solution is accomplished via Hamilton's equations, which in this context are also called the ray equations. These are evolution equations for the position x(t) and wavevector k(t) of the wave-packet [START_REF] Bühler | Waves and Mean Flows[END_REF].

For the spherical problem at stake, we use latitude ψ and longitude φ, with unit vectors e ψ and e φ . The conjugate momenta are then p ψ = ∂ ψ S and p φ = ∂ φ S. Using the expression for the gradient in terms of latitude and longitude, the wavevector is

k = p ψ R e ψ + p φ R cos ψ e φ , (3)
with R the radius of the Earth. The current velocity is u = u(ψ, φ)e φ + v(ψ, φ)e ψ , and

(1) then gives the Hamiltonian in terms of ψ, φ and their conjugate momenta:

Ω(ψ, φ, p ψ , p φ ) = g R p 1/2 + p ψ v(ψ, φ) R + p φ u(ψ, φ) R cos ψ , (4)
where

p def = p 2 ψ + p 2 φ cos 2 ψ 1/2 . ( 5
)
The spherical ray equations are then obtained from Ω(ψ, φ, p ψ , p φ ) via:

ψ = ∂ p ψ Ω = g R p ψ 2p 3/2 + v R , (6)
φ = ∂ p φ Ω = g R p φ 2p 3/2 1 cos 2 ψ + u R cos ψ , (7)
and

ṗψ = -∂ ψ Ω = - g R p 2 φ 2p 3/2 sin ψ cos 3 ψ - p ψ R ∂ ψ v - p φ R ∂ ψ u cos ψ , (8) ṗφ = -∂ φ Ω = - p ψ R ∂ φ v - p φ R ∂ φ u cos ψ . ( 9
)
An alternate method for deriving ray equations in spherical geometry is given in [START_REF] Hasha | Gravity wave refraction by threedimensionally varying winds and the global transport of angular momentum[END_REF]. In the special case of propagation through a still ocean, u = 0, the conjugate momentum p φ is constant and equations ( 6) through ( 9) reduce to those of [START_REF] Backus | The effect of the earth's rotation on the propagation of ocean waves over long distances[END_REF] and describe great circle propagation.

b. A special solution

An educational solution of the ray equations is obtained by considering the Cartesian case with a uniform current U flowing along the axis of x: see Fig. 5. If the source S is at the origin, and the receiver R is at x R = r cos α e x + r sin α e y , then the ray connecting S to R is a straight line, despite the Doppler shift corresponding to U. Thus, in planar geometry, a uniform current does not bend rays. This straight-line propagation, while simple in principle, is perhaps counterintuitive. Thus it is worthwhile to understand straight-line propagation through a uniform current by explicit solution of the ray equations. In cartesian geometry, the dispersion relation is Ω = √ gk+Uk x , with wave vector k = k x e x +k y e y and total wavenumber k = k 2 x + k 2 y . The cartesian ray equations are

ẋ = ∂ kx Ω = U + 1 2 k x k g k ẏ = ∂ ky Ω = 1 2 k y k g k , (10)
and kx = -∂ x Ω = 0 , ky = -∂ y Ω = 0 . (11)
The wave numbers k x and k y are constant, and the solution of ( 10) is therefore

x = (U + cos β v g ) t , y = sin β v g t , (12)
where (k x , k y ) = k(cos β , sin β) and v g = g/4k is the group velocity. Eliminating t between x and y in ( 12), and requiring that the ray pass through the receiver at x R , one finds that the direction, β, of the wave vector k is given by

v g sin(β -α) = U sin α . ( 13
)
The relevant case in oceanography is U ≪ v g , so that (13) can always be solved for β. In Figure 5 the slight inclination of the wave vector k to the straight-line ray path SR is the small angle β -α.

The importance of this simple solution is that it shows there is no relation between ray bending and the speed of currents.

The difference between α and β is a previously unremarked source of error for directional inferences which suppose that the wave vector k is precisely parallel to the direc-

β-α α x k U U U
Figure 5: A ray from a source S to a receiver R is not bent by a uniform current U. The wave vector k is inclined to the ray-path so that part of the group velocity compensates for advection by the current. For clarity, this schematic shows a large value of the angle β -α between the ray RS and k. Realistic surface currents are weak compared to the group velocity v g and hence β -α is at most one degree. tion of propagation e.g., as assumed by us in section 2, and previously by [START_REF] Snodgrass | Propagation of ocean swell across the pacific[END_REF] and [START_REF] Munk | Directional recording of swell from distant storms[END_REF]. With a typical velocity of oceanic surface currents U = 0.3 m/s and swell with 500 m wavelength, we obtain U/v g ≃ 0.02, hence α -β is of the order of 1 • : the direction β of the wave vector k is therefore a good, but not perfect, estimate of the direction of the straight ray between S and R. By contrast, in the following we show that non-uniform currents induce ray-bending, resulting in the direction of the wave vector k at the receiver being a poor estimate of the direction of the source, with errors often larger than 10 • .

c. Ray bending and the vertical vorticity of currents

Returning to the spherical case, in the absence of currents the solutions of ( 6) through ( 9) are great-circle geodesics (analogous to the straight line in Fig. 5) connecting the source to the receiver. Non-uniform currents will refract or bend the rays away from great-circle paths. In fact, it is the vertical vorticity of currents that is crucial for bending rays away from great circles [START_REF] Kenyon | Wave refraction in ocean currents[END_REF][START_REF] Dysthe | Refraction of gravity waves by weak current gradients[END_REF][START_REF] Landau | of a Course of Theoretical Physics[END_REF].

The connection between ray bending and vertical vorticity is simplest in the case when waves travel much faster than currents. For example, swell with a wavelength 500m has a group velocity v g = g/4k of 14m s -1 , which is much faster than the velocity of typical surface currents (at most 1m s -1 ). In this limit of fast wave-packets, the ray equations (6) through (9) reduce to the simpler and more insightful curvature equation

χ ≃ ξ v g , (14)
which is valid to first order in |u|/v g . Above, ξ is the vertical vorticity of the current,

ξ(x) def = ∂ φ v -∂ ψ (u cos ψ) R cos ψ , ( 15
)
and χ is the geodesic curvature: χ is the curvature of the trajectory projected onto the local horizontal plane. A curve with zero geodesic curvature is a great circle in the present context. Rays are thus deflected by vorticity just as the horizontal trajectory of a charged particle is bent by a vertical magnetic field: vorticity is analogous to the magnetic field. The magnetic-field analogy is the starting point of an alternate derivation of equation ( 14): in the limit of slow currents, the magnitude of k varies very little, and the group velocity of the waves remains almost constant, at the initial value v g(0) . Let us approximate the square root in the dispersion relation (1) by the parabola which has the same slope at the initial value of k:

Ω(x, k) ≃ 2 v 3 g(0) g |k| 2 + u•k + constant . ( 16
)
This quadratic expression is the same as the quadratic Hamiltonian for a negatively charged particle in a weak magnetic field ∇ × A(x), where A(x) is the vector potential. In dimensionless form, the Hamiltonian is:

H(x, p) = 1 2 |p + A| 2 + 1 2 |∇ × A| 2 , (17) ≃ 1 2 |p| 2 + A•p + O(A 2 ) , (18)
where p is momentum. The last expression corresponds to the weak-field limit. The charged particle experiences a Lorentz force, and its trajectory is well-known to have a local curvature proportional to the strength of the magnetic field χ = |∇ × A|/|p| [START_REF] Jackson | Classical Electrodynamics[END_REF]. For the wave-problem governed by the approximate dispersion relation ( 16), this translates into ray curvature |∇ × u|/v g(0) , which is the relation ( 14). The significance of ( 14) is that looking at a map of surface vorticity, one can assess which features will strongly deflect swell, and in which direction the rays will bend. The result also shows that wave propagation through this moving medium is isotropic, despite the direction determined by the velocity field u = u(ψ, φ)e φ + v(ψ, φ)e ψ . We use ( 14) to understand and interpret the results obtained by integration of the exact ray equations ( 6) through (9).

Deflection of swell by surface currents

Maps of surface currents u(x, t) are made available by satellite altimetry and scatterometry: the Ocean Surface Current Analysis in Real-time (OSCAR) dataset gives the surface-current velocity field with a spatial resolution of one-third of a degree and a temporal resolution of 5 days [START_REF] Bonjean | Diagnostic model and analysis of the surface currents in the tropical Pacific Ocean[END_REF]. OSCAR estimates the velocity u(x, t) required to determine the trajectory of swell by integration of the ray equations ( 6) through (9). We interpolate the OSCAR surface velocity field linearly onto a triangular mesh, before integrating the ray equations ( 6) through ( 9) with a forward Euler method.

One-third of a degree does not fully resolve mesoscale vorticity and thus computations based on OSCAR underestimate the deflection of gravity waves by currents. This underestimate indicates unambiguously that surface-current refraction is quantitatively sufficient to explain fluctuations of ±10 • in the direction of the incoming swell shown in Fig. 2(c) and 2(d).

a. A point source in the Southern Ocean

To model the observations of section 2, we begin by considering a point-source in the Southern Ocean, at latitude 58 • S and longitude 160 • E. This source emits surface waves in every direction, and with different wavelengths. For a given initial wavelength, we numerically integrate the exact ray-tracing equations ( 6) through ( 9), starting from the source point and with different initial orientations of the wavevector k. We assume that the emission is isotropic and so we shoot rays with the initial direction of k uniformly distributed on the unit circle, with a 3 × 10 -4 radian step (approximately one minute of arc). Six rays determined by this procedure are shown in Fig. 6. Out of all the emitted rays, an observer receives only the few rays which connect the source to the receiver. Thus we keep rays which reach the receiver within a radius of 30 nautical miles. The thick curve in Fig. 6 is an example of such a ray. Because the ray connecting the source to the receiver in Fig. 6 is not a straight line we conclude that swell is significantly deflected from a great-circle route by the OSCAR currents. For rays incident on the receiver, we define the deflection angle θ as the angle between the great-circle route and the direction of the ray at arrival. We use the convention θ < 0 if the ray at arrival is South of the great-circle route between source and receiver.

If an observer backtracks along the great circle indicated by the direction of arrival of the thick ray in Fig. 6, and determines the range with (2), then the inferred source is over sea-ice and almost on Antarctica. In other words, refraction by OSCAR currents has produced a mirage. Fig. 7 compares the left and right hand sides of equation ( 14) along a ray connecting the southern source to the receiver. The two curves coincide almost to within the line width. This validates the weak-current approximation used in ( 14) and shows that rays are bent from great circles only where there is strong surface vorticity i.e., in localized current systems such as the Antarctic Circumpolar Current (ACC), the subtropical frontal zone and the equatorial current system. The thick ray of Fig. 6 undergoes strong refraction when crossing these three features, and then travels more-or-less on great circles between these current systems. For example, after leaving the Southern Ocean the thick ray in Fig. 6 is on a great circle headed away from the receiver. But refraction by a large equatorial eddy subsequently bends the ray onto a great circle passing through the receiver.

We conclude that fluctuations in incident direction observed at the receiver are not the result of accumulation of many small random deflections. Instead, the model indicates that there are two or three large deflections of a ray that are associated with major hydrographic features of the surface current system.

b. An extended source in the Southern Ocean

For a given surface velocity field, several rays -none of which are great circlesconnect the source to the receiver. An example is shown in Figure 8. We refer to the collection of rays that connect the source to the receiver as a "multipath". As we did with the directional data in Fig. 2, an observer can construct the "average inferred source" at the location of the mean position of the multiple inferred sources. This mean deflection is denoted as θ , where the brackets denote averaging over all rays. For a point-source of swell, the average inferred source is usually significantly displaced from the actual source, with values of the order of ±5 • for swell with initial wavelength λ = 500m. An extreme example is presented in figure 8, where refraction by surface currents is so strong that the average inferred source is on land. However, real storms have a typical extension of several hundreds of kilometers and should not be considered to be point-sources.

To determine the influence of this spatial extension, we model an extended source by a disk of radius 4 • , or 440 km, around a center located at latitude 57 • S and longitude 172 • E. We assume a uniform density of incoherent point sources inside this disk. It is then easier to solve numerically the backward problem: we shoot rays "backwards" from the receiver with a constant step of 10 -4 radian in initial direction (approximately twenty seconds of arc). We keep only the rays that have a nonzero intersection with the swellemitting disk in the Southern Ocean. To compute the direction of the average inferred source, we perform a weighted average of the angular directions of the rays that intersect the extended source. The weight of each ray is equal to the length of its intersection with the swell-emitting disk: the contributions from incoherent sources inside the disk add up on each ray. The insert in figure 8 shows a typical ray pattern obtained with this extended source for swell with wavelength λ = 500m. A wide bundle of rays connects the source to the receiver. We find that the average inferred source is close to the center of the actual source: contributions to the average deflection from the different points inside the emitting-disk average out to almost zero. Going from one 5-day frame of the OSCARdata to the next one, the average deflection θ fluctuates weakly around zero, with a rootmean-square (rms) value of the order of 3 • for swell with 500m initial wavelength: most often, the average inferred source is inside the swell-emitting disc. Assuming gaussian statistics for the average deflection θ , such a low 3 • value of the rms fluctuations means that an average deflection of 10 • or higher has a probability lower than 0.04. This does not rule out the possibility of a large θ occasionally occurring because the OSCAR data underestimate the actual vorticity of surface currents and because 0.04 is not 0: if we analyze 25 events we might hope to see one example of a 10 • average deflection. This is consistent with the analysis from section 2: mirages due solely to the effect of surface currents are rare events. An additional ingredient seems therefore necessary to explain the more frequent occurrence of strong average deflections in the observations by Munk et al., together with the systematic displacement of the inferred sources towards the South of the actual storms. We return to this in section 6.

Fluctuations in the direction of swell

We return now to the ±10 • directional fluctuations in Fig. 2(c) and 2(d). Previously, to estimate the location of the source, we removed these fluctuations by averaging. But in this section we investigate the hypothesis that the frequency dependence of directional fluctuations contains information about the strength of the surface vorticity field.

Let us assume that, during a storm, surface wave packets are emitted somewhat randomly along the different rays that connect the storm to the receiver. These different wave packets will not reach the receiver at the same time, and thus the direction of the incident swell should fluctuate in time, with fluctuations of the order of typical values of θ visible in Fig. 8. Our hypothesis is that propagation of wave packets along the various component rays of the multipath in Fig. 8 results in the ±10 • fluctuations in incident angle shown in Fig. 2(c) and 2(d). To investigate this further we first characterize the fluctuations in direction induced by surface currents, using the twenty years of OSCAR data to perform a Monte Carlo computation. We then compare this prediction to the pitch-and-roll buoy data.

a. Directional fluctuations in model based on OSCAR currents

To gather statistics on the deflection of the rays, we repeat the point-source analysis of section 4 for each 5-day OSCAR surface velocity field recorded between October 1992 and October 2011. From this extensive simulation, we compute the probability density function (PDF) of the deflection angle, θ, which is defined as the angle between the ray at arrival and the great-circle route. This PDF of θ is shown in figure 9(a) using four different values of initial swell wavelength λ. The root-mean-square (rms) deflection is around 17 • for wavelength λ = 250m. The rms deflection decreases with increasing wavelength. Indeed, longer waves are faster and thus less refracted according to the curvature equation ( 14): when considered as a function of θ √ λ, the four PDFs collapse onto a single master curve: see Fig. 9(b). This collapse indicates an rms deflection angle proportional to the frequency of the swell, with θ rms in degrees and f in Hertz.

θ rms = 218f , (19)
The order of magnitude of typical deflection θ can be understood as follows: according to the curvature equation ( 14), when a ray crosses an eddy with size ℓ and typical vorticity ξ 0 (see Fig. 10), the direction of propagation is deflected by an angle

θ ∼ ξ 0 ℓ v g (20) = 4π ξ 0 ℓ g f , (21)
where v g = g/(4πf ) has been used. Out of the three current systems that refract the swell, the ACC has the most intense vorticity (large ξ 0 ) whereas the equatorial current has the largest eddies (large ℓ). The subtropical frontal zone has weaker eddies with small sizes and can be neglected for this rough estimate. Because the ACC is close to the source and far from the receiver, the ACC has a smaller impact on the deflection θ than the equatorial current system: even after strong refraction in the ACC, rays escape this first current system rather close to the source point. If the ACC were the only vortical refractor between source and receiver, then the observer would make only a small error in inferring the position of the source. But the equatorial current system is close to the receiver. Equatorial eddies bend some rays so that they hit the receiver, which results in large values of θ. Hence an order of magnitude of θ is given by the angular shift due to a few, or even one, eddy in the equatorial current system. With λ = 500m, and using ℓ = 400km and ξ 0 = 10 -5 s -1 for a typical eddy in the equatorial current system, one obtains from (21) θ ≃ 9 • . Although Kenyon's intuition that surface currents refract the swell was correct, he focused on the time-averaged ACC and neglected both the equatorial current and the role of the eddies. The latter are crucial to explain the directional fluctuations due to OSCAR currents. Indeed, to determine the role of the average currents, we time-averaged the twenty years of OSCAR data before simulating the same point-source as in section 4a: the time-averaged surface currents produce a negligible deflection |θ| of the rays (typically 1.5 • for λ = 500 m), which shows that the strong fluctuations visible in figure 9 are due to the small-scale eddies, particularly their intense vorticity, and not to the time-averaged currents.

b. Directional fluctuations in the pitch-and-roll buoy data

To gather statistics on the fluctuations of θ in the buoy data shown in Fig. 2, we analyzed 18 swell events recorded between 2004 and 2007. These 18 events were selected because the signal was particularly clean i.e., there was no evidence of multiple sources in any of these 18 events. For each swell event, we removed the mean value and the linear trend in the incident angle signal. We consider the remaining fluctuating θ(t) as a function of the peak frequency f (t) instead of time. We divide the signal into small bins in f , and we compute the root-mean-square θ for each bin. The result is displayed in Fig. 11, together with the predictions from the analysis of the OSCAR data. The rms value of θ indeed shows a linear dependence with frequency, as predicted for refraction by surface currents. More surprisingly, the prefactor of the linear law is very well captured by the OSCAR data analysis. We were expecting to underestimate the effect of surface currents using the coarse OSCAR data so the agreement in Fig. 11 might be fortuitous. But the orders of magnitude are definitely compatible.

The good comparison in Fig. 11 is also hostage to errors in directional measurements at NOAA station 46086. A significant uncertainty is that the accelerometers on NDBC 3m discus buoys (such as buoy 46086) are known to be a little noisy, or at least they are noisier than accelerometers on Datawell Directional Waveriders [START_REF] O'reilly | A comparison of directional buoy and fixed platform measurements of Pacific swell[END_REF]. The discus-buoy noise is known to bias estimates of spread, but does not influence estimates of mean direction [START_REF] O'reilly | A comparison of directional buoy and fixed platform measurements of Pacific swell[END_REF]. To further allay these concerns, the November 2007 case study has been repeated by Sean Crosby (personal communication) using measurements made by a Datawell buoy deployed at NOAA NDBC station 46232 (32.530 • N and 117.431 • W). Datawell buoys have the advantage of measuring accurately both the mean direction and the spread of the swell, using 26-minute-long samples. The Datawell directional measurements are consistent with the NDBC results: directional spread increases in time as the peak frequency increases. The fluctuations in incident angle (or "mean direction" signal) are greater than the Datawell noise level, and their magnitude increases between the beginning and the end of the swell event. These results are consistent with our hypothesis that some part of the observed directional fluctuations is due to ray bending by surface currents.

Discussion and conclusion

Because of refraction by surface currents, a storm and a receiving station are connected not only by the great-circle route, but by a multipath: a bundle of rays with an angular width which is much larger than the angular width of the storm. Wave packets travel on these many rays before reaching the receiver, which leads to strong temporal fluctuations in the incoming direction measured at the receiver. The root mean square fluctuations in the directional signal at pitch-and-roll buoy 46086 are consistent with predictions using ray-tracing through the mesoscale vorticity of the OSCAR dataset.

Most often, the fluctuations average out in time and the mean direction measured by the buoy over several days coincides with the actual direction of the storm. However, mirages do occur. These rare events are observed when most of the swell from a single storm is deflected in the same direction. The mirage effect can be greatly enhanced by topography: if topographic features close to the great-circle route obstruct part of the ray bundle, then the average direction of the incoming swell can be very different from the direction of the storm.

The historical case of [START_REF] Munk | Directional recording of swell from distant storms[END_REF] probably results from such an interplay between surface currents and topography, summarized following [START_REF] Munk | Corrigendum to "Directional recording of swell from distant storms[END_REF] in the simplified schematic Fig. 12. We represent the great-circle route between a storm and the receiver off San Clemente island, together with the extreme rays experiencing a refraction of respectively +10 • and -10 • by mesoscale vorticity in the equatorial current. Now close to the receiver is Cortez bank, a shallow bathymetric structure that selectively blocks rays that travel North of the great-circle route without affecting rays that approach from the South. An observer therefore measures mostly wave packets coming from South of the great-circle route, and infers a source on Antarctica.

Satellite observations provide modern confirmation of long-range propagation of ocean swell [START_REF] Heimbach | Development and application of satellite retrievals of ocean wave spectra[END_REF][START_REF] Collard | Monitoring and analysis of ocean swell fields from space: New methods for routine observations[END_REF]. But the source-location problem remains an issue: recent inferences of storm sources follow [START_REF] Munk | Directional recording of swell from distant storms[END_REF] and backtrack along great circles. This ignores refraction by currents and the resulting increase in the width of the swell beam. It is interesting that through the curvature formula ( 14), vorticity ξ emerges as a key environmental variable that controls refraction of surface gravity waves. NASA's Surface Water and Ocean Topography satellite promises to greatly improve resolution of mesoscale vorticity [START_REF] Fu | SWOT: The Surface Water and Ocean Topography Mission[END_REF][START_REF] Fu | Observing oceanic submesoscale processes from space[END_REF]: better wave forecasts might be an unexpected outcome of this mission. A ray from a source S to a receiver R is not bent by a uniform current U.

List of Figures

The wave vector k is inclined to the ray-path so that part of the group velocity compensates for advection by the current. For clarity, this schematic shows a large value of the angle β -α between the ray RS and k. Realistic surface currents are weak compared to the group velocity v g and hence β -α is at most one degree. . . . . . . . . . . . . . . . . . . . . . . . . 9 6

Swell with λ = 500m emitted by a point-source in the Southern Ocean. This figure uses an azimuthal equidistant projection so that great circles passing through the receiver (and only these) appear as straight lines i.e., the dashed straight lines are great circles. The color scale shows vertical vorticity of OSCAR surface currents, ξ in s -1 . The thick-curve is a ray that connects the source to the receiver. The thin curves are five other rays.

Figure 1 :

 1 Figure 1: Identification of two swell events recorded at NOAA Station 46086 in the San Clemente basin. The "ramps" in a time-frequency plot signal the arrival of swell from a distinct, distant source. Panel (a) shows the peak spectral frequency as a function of time recorded for 14 days in November 2007; we study the event occurring between the 6th and 10th of November. Panel (b) shows the peak frequency as a function of time recorded for 7 days in April 2007; we study the event between the 31st of March and the 3rd of April.

Figure 2 :

 2 Figure 2: Two swell events recorded at NOAA Station 46086 in the San Clemente basin. Panels (a), (c) and (e) correspond to a swell event recorded in November 2007, and panels (b), (d) and (f) to a swell event recorded in April 2007. The top panels show the estimated peak frequency f as a function of time. f increases linearly with time due to dispersive propagation of surface waves. The middle panels show the incident angle measured clockwise from North by the buoy. The incident angle fluctuates around a mean value of 204 • for panel (c) and 224 • for panel (d). The bottom panels show the SWH in meters.

Figure 3 :

 3 Figure 3: Color contours indicate significant wave height (SWH) in meters from the ECMWF ERA reanalysis on October 30th, 2007, at midnight GMT, shown using a South Polar projection. The thick black line is the sea-ice limit. The solid black grid shows great-circle routes from the NOAA station 46086, and lines of constant range from this station. The red spot, which is very close to the region of maximum SWH, indicates the source inferred from swell recorded at Station 46086.

Figure 4 :

 4 Figure 4: A mirage: the swell recorded at NOAA station 46086 seems to originate from New Zealand. The inferred location (the red dot) is displaced by approximately 10 • of arc from the region of maximum SWH. Color contours indicate SWH in meters for the ECMWF ERA reanalysis for the storm of March 23rd, 2007 at 6am GMT. The solid black grid shows great-circle routes from the NOAA station 46086, and lines of constant range from this buoy.

Figure 6 :Figure 7 :

 67 Figure6: Swell with λ = 500m emitted by a point-source in the Southern Ocean. This figure uses an azimuthal equidistant projection so that great circles passing through the receiver (and only these) appear as straight lines i.e., the dashed straight lines are great circles. The color scale shows vertical vorticity of OSCAR surface currents, ξ in s -1 . The thick-curve is a ray that connects the source to the receiver. The thin curves are five other rays. Using great-circle backtracking (red dashed line), the inferred source is far from the true source.

Figure 8 :

 8 Figure 8: A bundle of rays connects the source to the receiver (λ = 500m). Blue rays reach the receiver within 30 nautical miles. The red spot indicates the average inferred source. Insert: same analysis performed with an extended source of radius 4 • , delimited by the red circle. The average inferred source coincides with the center of the circle within 2 • (not shown). This figure uses an azimuthal equidistant projection centered on the receiver and the color scale shows vertical vorticity of OSCAR surface currents, ξ in s -1 .

Figure 9 :

 9 Figure 9: (a) Probability density function of the deflection angle. The typical deflection is around 10 • , which is consistent with the measurements of Munk et al. and the observation of section 2. Waves with larger wavelength are faster and less refracted. (b) The four PDFs collapse onto a master PDF when considered a function of θ √ λ.

Figure 10 :

 10 Figure 10: Bending of a beam of parallel rays by a strong Gaussian vortex ξ = ξ 0 exp(-(r/ℓ) 2 ), with ξ 0 / √ gk = 5/π and ℓk = 40. In this illustration the vortex is unrealistically strong in order to show the ray pattern; realistic vortices produce much smaller deflections.

Figure 11 :

 11 Figure11: Root-mean-square fluctuations in incoming swell direction measured by a pitch-and-roll buoy, averaged over 18 storms, as a function of swell frequency f . Error bars are evaluated using the square root of the number of values in each frequency bin. The dashed-line is the result (19) from the analysis of OSCAR data.

Figure 12 :

 12 Figure 12: Interplay between surface currents and shallow bathymetry. (a) The black line is the great-circle route. The blue and green rays connect the source to the receiving station from Munk et al. (1963). They are deflected by ±10 • by mesoscale vorticity at the equator. The Northern ray hits Cortez bank, and does not reach the receiver. An observer sees only the Southern ray and great-circle backtracking puts the source on Antarctica (dashed red line). (b) Local bathymetry in meters near San Clemente island. Cortez bank blocks the Northern ray.

  Using great-circle backtracking (red dashed line), the inferred source is far from the true source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 7 For the thick ray in Fig. 6 the ray curvature (blue) is almost equal to the surface vorticity divided by group velocity (red dashed). The smooth black curve is low-pass filtered curvature. The strongest refraction occurs when the ray crosses the ACC and equatorial current. Weaker refraction occurs as the ray transits the subtropical frontal zone. . . . . . . . . . . . 13 8 A bundle of rays connects the source to the receiver (λ = 500m). Blue rays reach the receiver within 30 nautical miles. The red spot indicates the average inferred source. Insert: same analysis performed with an extended source of radius 4 • , delimited by the red circle. The average inferred source coincides with the center of the circle within 2 • (not shown). This figure uses an azimuthal equidistant projection centered on the receiver and the color scale shows vertical vorticity of OSCAR surface currents, ξ in s -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 9 (a) Probability density function of the deflection angle. The typical deflection is around 10 • , which is consistent with the measurements of Munk et al. and the observation of section 2. Waves with larger wavelength are faster and less refracted. (b) The four PDFs collapse onto a master PDF when considered a function of θ √ λ. . . . . . . . . . . . . . . . . . . . . 16 10 Bending of a beam of parallel rays by a strong Gaussian vortex ξ = ξ 0 exp(-(r/ℓ) 2 ), with ξ 0 / √ gk = 5/π and ℓk = 40. In this illustration the vortex is unrealistically strong in order to show the ray pattern; realistic vortices produce much smaller deflections. . . . . . . . . . . . . . . . 17 11 Root-mean-square fluctuations in incoming swell direction measured by a pitch-and-roll buoy, averaged over 18 storms, as a function of swell frequency f . Error bars are evaluated using the square root of the number of values in each frequency bin. The dashed-line is the result (19) from the analysis of OSCAR data. . . . . . . . . . . . . . . . . . . . . . . . . . . 18 12 Interplay between surface currents and shallow bathymetry. (a) The black line is the great-circle route. The blue and green rays connect the source to the receiving station from Munk et al. (1963). They are deflected by ±10 • by mesoscale vorticity at the equator. The Northern ray hits Cortez bank, and does not reach the receiver. An observer sees only the Southern ray and great-circle backtracking puts the source on Antarctica (dashed red line). (b) Local bathymetry in meters near San Clemente island. Cortez bank blocks the Northern ray. . . . . . . . . . . . . . . . . . . . . . . . . 19

1

  Identification of two swell events recorded at NOAA Station 46086 in the San Clemente basin. The "ramps" in a time-frequency plot signal the arrival of swell from a distinct, distant source. Panel (a) shows the peak spectral frequency as a function of time recorded for 14 days in November 2007; we study the event occurring between the 6th and 10th of November. Panel (b) shows the peak frequency as a function of time recorded for 7 days in April 2007; we study the event between the 31st of March and the 3rd of April. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The top panels show the estimated peak frequency f as a function of time. f increases linearly with time due to dispersive propagation of surface waves. The middle panels show the incident angle measured clockwise from North by the buoy. The incident angle fluctuates around a mean value of 204 • for panel (c) and 224 • for panel (d). The bottom panels show the SWH in meters. . . . . . . . . . . . . . . . . . . . . . . 4 3 Color contours indicate significant wave height (SWH) in meters from the ECMWF ERA reanalysis on October 30th, 2007, at midnight GMT, shown using a South Polar projection. The thick black line is the seaice limit. The solid black grid shows great-circle routes from the NOAA station 46086, and lines of constant range from this station. The red spot, which is very close to the region of maximum SWH, indicates the source inferred from swell recorded at Station 46086. . . . . . . . . . . . . . . . 5 4 A mirage: the swell recorded at NOAA station 46086 seems to originate from New Zealand. The inferred location (the red dot) is displaced by approximately 10 • of arc from the region of maximum SWH. Color contours indicate SWH in meters for the ECMWF ERA reanalysis for the storm of March 23rd, 2007 at 6am GMT. The solid black grid shows great-circle routes from the NOAA station 46086, and lines of constant range from this buoy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 5

		3
	2	Two swell events recorded at NOAA Station 46086 in the San Clemente
		basin. Panels (a), (c) and (e) correspond to a swell event recorded in
		November 2007, and panels (b), (d) and (f) to a swell event recorded in
		April 2007.

Acknowledgments. This research was supported by the National Science Foundation under OCE10-57838; BG was partially supported by a Scripps Postdoctoral Fellowship. We thank Ryan Abernathey, Fabrice Ardhuin, Oliver Bühler, Fabrice Collard, Sean Crosby, Falk Feddersen and particularly Walter Munk, for many conversations and help with this problem.

The ocean surface current data used in this study was provided by the OSCAR Project Office. The significant wave height data were obtained from the ECMWF data server. The pitch-and-roll buoy data is from the National Data Buoy Center, which is operated by the National Oceanic and Atmospheric Administration.