Gangloff Aurélien
email: aurelien.gangloff@ifremer.fr

Verney Romaric

Doxaran David

Ody Anouck

Estournel Claude

Aurélien Gangloff

Investigating Rhône River plume (Gulf of Lions, France) dynamics using metrics analysis from the MERIS 300m Ocean Color archive (2002 -2012)

Keywords: river plume metrics, sediment dynamics, MERIS, ocean color data, suspended particulate matter, Gulf of Lions

In coastal environments, river plumes are major transport mechanisms for particulate matter, nutriments and pollutants. Ocean colour satellite imagery is a valuable tool to explore river turbid plume characteristics, providing observations at high temporal and spatial resolutions of suspended particulate matter (SPM) concentration over a long time period, covering a wide range of hydro-meteorological conditions. We propose here to use the MERIS-FR (300 m) Ocean Colour archive (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) in order to investigate Rhône River turbid plume patterns generated by the two main forcings acting on the north-eastern part of the Gulf of Lions (France): wind and river freshwater discharge. Results are exposed considering plume metrics (area of extension, south-east-westernmost points, shape, centroid, SPM concentrations) extracted from satellite data using an automated image-processing tool. Rhône River turbid plume SPM concentrations and area of extension are shown to be mainly driven by the river outflow while wind direction acts on its shape and orientation. This paper also presents the region of influence of the Rhône River turbid plume over monthly and annual periods, and highlights its interannual variability.

Highlights

► A new method to extract plume metrics from ocean color satellite images is proposed. ► Use of innovative plume metrics such as centroids and schematic shape (skeleton). ► Analysis of the MERIS-300m database covering the 2002-2012 period. ► Grand Rhone River turbid plume metrics are highly correlated to river discharge. ► Wind regime (Onshore/offshore) drives plume patterns.

Introduction

Understanding particulate matter dynamics over continental margins is a challenging issue for environmental purposes as fine particles are vectors of pollutants such as heavy metals or radionuclides. In the North-Eastern (NE) part of the Gulf of Lions (GoL), located in the North-Western (NW) part of the Mediterranean Sea (Fig. 1), the Rhône River is the main source of sediments [START_REF] Bourrin | Contribution to the study of coastal rivers and associated prodeltas to sediment supplyi in Gulf of Lions (NW Mediterranean sea)[END_REF]. These particles can be trapped within the prodelta located just downstream from the river mouth [START_REF] Radakovitch | 210Pb and Caesium Accumulation in the Rhône Delta Sediments[END_REF][START_REF] Roussiez | Early muddy deposits along the Gulf of Lions shoreline: A key for a better understanding of land-to-sea transfer of sediments and associated pollutant fluxes[END_REF] or advected along the shelf as suspended particulate matter (SPM), directly within the river plume or indirectly, after wave-induced resuspension events, occurring mainly during winter storms [START_REF] Dufois | Impact of winter storms on sediment erosion in the Rhone River prodelta and fate of sediment in the Gulf of Lions (North Western Mediterranean Sea)[END_REF]. This work focuses on the SPM advected directly from the Rhône River mouth into the GoL through the plume dynamics.

The Rhône is one of the major rivers of the Mediterranean basin, the most important in the GoL, flowing from Swiss Alps (Lac Leman) to the Mediterranean Sea, where it divides in two branches, about 50 km upstream the main river mouth: the Grand Rhône and the Petit Rhône (see Fig. 1). The river drains continental waters over a 97800 km 2 catchment area, covering different climatic zones such as alpine, oceanic and Mediterranean ones [START_REF] Maillet | Morphological changes and sedimentary processes induced by the December 2003 flood event at the present mouth of the Grand Rhône River (southern France)[END_REF]. The mean annual discharge of the Rhône River is 1700 m 3 .s -1 , and about 90 % (~ 1500 m 3 .s -1 ) is distributed to the GoL via the Grand Rhône [START_REF] Fanget | Historical evolution and extreme climate events during the last 400years on the Rhone prodelta (NW Mediterranean)[END_REF]. Annual, decennial, centennial and millennial floods water discharges for Rhône River are of about 5000, 8400, 11300 and 14100 m 3 .s -1 , respectively [START_REF] Maillet | Morphological changes and sedimentary processes induced by the December 2003 flood event at the present mouth of the Grand Rhône River (southern France)[END_REF], leading to about 4500, 7600, 10200 and 12700 m 3 .s -1 for the Grand Rhône River, respectively. For the following study, we chose to focus on the Grand Rhône River plume only.

The Rhône River plume dynamics has been studied for several years through observations like sea temperature monitoring with TIROS-N satellite database by [START_REF] Demarcq | La dynamique superficielle du panache du Rhône d'après l'imagerie infrarouge satellitaire[END_REF], showing that the plume response time is short, about 5 h for average river discharges (Q wtypically around 1500 m 3 .s -1 ) and increases with river outflow, above 24 h for river discharge values of 4350 m 3 .s -1 for instance. They also highlighted that the plume is clearly deflected on the right side of the wind, with an average angle of deflection of 50°. However, the thermal plume depth sensed by AVHRR is very thin (O(10m)) compared to turbid plume depth (O(0.1m) to O(1m)), so results are not fully comparable, but provide complementary knowledge about the near surface dynamics. [START_REF] Forget | VHF remote sensing for mapping river plume circulation[END_REF] also showed, through VHF radar analysis, that for low river discharges (500 to 800 m 3 .s -1 ) and calm wind conditions, summer breezes (< 5 m.s -1 ) are able to push the plume front toward the river mouth over a distance of about 10 km in a few hours. [START_REF] Broche | Experimental study of the Rhone plume. Part I: physics and dynamics[END_REF], relying on field and VHF data measurements, reported that flow intensity variations have a weak influence on the plume's shape but act on its vertical thickness and its internal velocity variations. They also highlighted that strong western winds can deflect the plume eastward within a few hours. The Rhône River plume dynamics was also investigated through numerical modelling. [START_REF] Estournel | The plume of the Rhone: numerical simulation and remote sensing[END_REF] showed that the river plume is highly sensitive to wind, yet two types of plumes were brought to the fore, associated to the two predominant wind regimes in the climate of the Gulf of Lions. The typical plume shapes observed are the following: detached from the coast and orientated southwestward, or thin and running along the coast, orientated westward, for northern and southeastern winds, respectively. [START_REF] Marsaleix | A numerical study of the formation of the Rhône River plume[END_REF] added that for simulations forced by northwestern winds of different intensities, the seaward deflection is significant for wind velocities higher than 7 m.s -1 . These studies focused on freshwater river plume dynamics and not on the turbid plume one, which could not necessarily behave identically. Also, they focused on short time periods (daily to monthly analysis). In the present paper, we will investigate the turbid plume dynamics over a 10 years database, allowing event to decennial analysis. Recent works on the studied area attempted to characterize the turbid plume. For high river discharge (~ 4000 m 3 .s -1 ) and calm conditions (northwestern wind with velocity < 3 m.s -1 ), [START_REF] Many | Particle assemblage characterization in the Rhone river ROFI[END_REF] observed a plume expanding up to 30 km off the coast, with a thickness decreasing seaward from 20 m to 5 m, as well as for its depth-averaged SPM concentration, decreasing from 15 down to 5 mg.l -1 . Using satellite ocean colour data, [START_REF] Lorthiois | Daily and seasonal dynamics of suspended particles in the Rhône River plume based on remote sensing and field optical measurements[END_REF] showed that the duration of a flood event has an influence on the turbid plume vertical thickness. Their results revealed that during a strong river discharge event (15 consecutive days with Q w > 2000 m 3 .s -1 ), the plume reached a thickness of ~ 5 m, which was not the case for a shorter event (3 consecutive days with Q w > 2000 m 3 .s -1 ), where the plume thickness was of ~ 1 m.

In the present paper, we propose to investigate the Rhône River turbid plume patterns and dynamics from the decennial scale down to the event scale, in relation with the main forcings prevailing on the area (Fig. 1). These forcings are: (i) the strong continental winds called Mistral (from the North) and Tramontane (from the North-West), which curl drives the coastal currents direction [START_REF] Estournel | Observation and modeling of the winter coastal oceanic circulation in the Gulf of Lion under wind conditions influenced by the continental orography (FETCH experiment)[END_REF]; (ii) the Rhône River discharge and (iii) the general circulation following the slope of the GoL shelf, called the Northern Curent [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF], south the area of interest. The study area, which is located on the inner shelf of the GoL, is a wave-dominated area [START_REF] Dufois | Wave-and current-induced bottom shear stress distribution in the Gulf of Lions[END_REF] with low tidal range of about 0.3 m [START_REF] Fanget | Historical evolution and extreme climate events during the last 400years on the Rhone prodelta (NW Mediterranean)[END_REF].

Investigating turbid plume patterns and their relation with forcings requires to have a long time-period dataset, representative of the meteorological/hydrological conditions. Satellite imagery has proven to be a valuable tool to assess turbidity or SPM concentration from ocean colour data in estuaries [START_REF] Hudson | Using Satellite Observations to Characterize the Response of Estuarine Turbidity Maxima to External Forcing[END_REF] and coastal zones such as the North Sea [START_REF] Gohin | Annual cycles of chlorophyll-a, non-algal suspended particulate matter, and turbidity observed from space and in-situ in coastal waters[END_REF][START_REF] Nechad | Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[END_REF], the Bay of Biscay [START_REF] Gohin | Satellite-derived parameters for biological modelling in coastal waters: Illustration over the eastern continental shelf of the Bay of Biscay[END_REF], the Danube River [START_REF] Constantin | Estimation of water turbidity and analysis of its spatio-temporal variability in the Danube River plume (Black Sea) using MODIS satellite data[END_REF][START_REF] Güttler | Turbidity retrieval and monitoring of Danube Delta waters using multi-sensor optical remote sensing data: An integrated view from the delta plain lakes to the western-northwestern Black Sea coastal zone[END_REF], the Adour River [START_REF] Petus | Estimating turbidity and total suspended matter in the Adour River plume (South Bay of Biscay) using MODIS 250-m imagery[END_REF][START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF], the Gironde estuary [START_REF] Doxaran | Monitoring the maximum turbidity zone and detecting fine-scale turbidity features in the Gironde estuary using high spatial resolution satellite sensor (SPOT HRV, Landsat ETM+) data[END_REF][START_REF] Doxaran | Dynamics of the turbidity maximum zone in a macrotidal estuary (the Gironde, France): Observations from field and MODIS satellite data[END_REF] or the Gulf of Lions [START_REF] Lorthiois | Daily and seasonal dynamics of suspended particles in the Rhône River plume based on remote sensing and field optical measurements[END_REF][START_REF] Ody | Potential of High Spatial and Temporal Ocean Color Satellite Data to Study the Dynamics of Suspended Particles in a Micro-Tidal River Plume[END_REF] for instance. Ocean colour satellite imagery can be used to study turbid plume metrics such as its extension area in relation with hydrological and meteorological forcings, as it has already been done for the Mississippi River plume [START_REF] Walker | Satellite assessment of Mississippi River plume variability: Causes and predictability[END_REF][START_REF] Walker | Effects of river discharge, wind stress, and slope eddies on circulation and the satellite-observed structure of the Mississippi River plume[END_REF], the Adour River plume [START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF] or the Ebro River plume [START_REF] Fernández-Nóvoa | Analysis of the influence of river discharge and wind on the Ebro turbid plume using MODIS-Aqua and MODIS-Terra data[END_REF].

Based on the 10-year MERIS-300m archive, the main objectives of this study are:

1. to report the inter-annual and seasonal variability of the Rhône River turbid plume in terms of area of extension, averaged SPM concentration and maximum SPM concentration; 2. to highlight turbid plume patterns in response to typical forcings acting on the study area, exploring various metrics, such as its area of extension, its centroid, its shape or its extreme points (southernmost, westernmost and easternmost).

Materials and methods

MERIS-300m database

Entering the water, solar irradiance interacts with its constituents and a part of the signal is scattered backward out of the water column. Ocean colour sensors measure the spectral characteristics of this backscattered signal. This water-leaving radiance (in W.m -2 .nm -1 .sr -1 ) can be normalized by the downward irradiance (in W.m -2 .nm -1 ) just above the surface at different wavelengths. This normalized radiance is called the remote sensing reflectance (R rs , in sr -1 ) and can be expressed in terms of inherent optical properties (IOPs), such as the total backscattering and absorption coefficients, depending on the nature and concentration of a substance. After atmospheric correction processes, it is then possible to retrieve the concentration of substances in the surface water layer sensed by ocean colour sensors (optical depth). In the present study, the medium-spectral resolution imaging spectrometer MERIS-300m seawater reflectance data was used to evaluate SPM concentrations on the study area (Fig. 1). Archived data are available from 2002 to 2012 with a global coverage every 3 days.

Data was collected via the online GIS COOC data portal (http://kalicotier.gis-cooc.org/), providing two reflectance products of MERIS Full Resolution. The first one is based on the atmospheric correction processor MEGS8.1 (third reprocessing of the MERIS Ground Segment prototype processor), standard processor for the MERIS sensor [START_REF] Müller | The Ocean Colour Climate Change Initiative: I. A methodology for assessing atmospheric correction processors based on in-situ measurements[END_REF] which uses a neural network algorithm applied for atmospheric correction specific for the retrieval of case2 water constituents [START_REF] Doerffer | The MERIS Case 2 water algorithm[END_REF]. The second one is based on the SAABIO (Semi-Analytical Atmospheric and Bio-Optical) processor [START_REF] Gernez | Remote sensing of suspended particulate matter in turbid oyster-farming ecosystems[END_REF], using MEGS8.1 with a different bright pixel atmospheric correction and without the case2 regional neural network for computing SPM concentration.

The standard SPM product, available from the MEGS8.1 neural network [START_REF] Doerffer | The MERIS Case 2 water algorithm[END_REF] was tested, as well as three other regional algorithms identified hereafter as SPM1, SPM2 and SPM3. SPM1 and SPM2 were derived from field data collected during the TUCPA (2014) and the PLUMRHO (2015) experiments, as well as SPM3, following the formulation of the generic model developed by [START_REF] Nechad | Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[END_REF]. These algorithms use the Rrs at 620 nm and/or 655 nm and are the following: a linear regression, associating SPM concentration (mg.l -1 ) to Rrs in the 665 nm band [START_REF] Ody | Potential of High Spatial and Temporal Ocean Color Satellite Data to Study the Dynamics of Suspended Particles in a Micro-Tidal River Plume[END_REF], another linear regression, using the mean of the 665 nm and the 620 nm bands, and a non-linear regression, derived from the generic model of [START_REF] Nechad | Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[END_REF] and using the 665 nm band, where is the water leaving reflectance at 665 nm, coefficient C = 0.1728 was theoretically calculated for the MERIS 665 nm band in [START_REF] Nechad | Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[END_REF], and coefficients A = 213 and B = -0.0785 were optimized for our region of interest considering data collected during the two field experiments mentioned above.

The performances of the standard SPM product and the regional SPM algorithms were evaluated using SPM matchups. These matchups are in situ measurements close in time (< 6 h) and space (nearest pixel) to the satellite acquisition and were obtained during field campaigns from 2010 to 2012 (Table 1). The same sampling and filtering protocol was used for each campaign. They cover a range of SPM concentration from 0.8 mg.l -1 to 32 mg.l -1 , with a mean of 8.2 mg.l -1 . The MEGS8.1 standard product showed very good performances with a mean relative error ε r of 40.0 %, a R 2 of 79.4 and root mean square errors RMSE for SPM concentrations ≤ 10 mg.l -1 and > 10 mg.l -1 of 2.56 and 6.45 mg.l -1 , respectively. The semi-empirical relation described by [START_REF] Nechad | Calibration and validation of a generic multisensor algorithm for mapping of total suspended matter in turbid waters[END_REF] for the SAABIO product also showed good statistical results with a ε r of 37.0 %, a R 2 of 79.3 % and RMSE for SPM concentrations ≤ 10 mg.l -1 and > 10 mg.l -1 of 2.02 and 6.52 mg.l -1 , respectively. Both algorithm results were close and, in order to minimize data processing, the standard SPM product was chosen for this study. A step-by-step statistical method (Fig. 2) was applied to raw SPM concentrations maps in order to discard images with too many flagged pixels on Rhône River turbid plume as follows:

1. Over all available images, a statistical analysis of SPM concentration was performed on the area of interest to determine a SPM threshold in order to identify the turbid plume from the background SPM concentration. For all images, SPM concentration distribution was analysed and the percentile 95 (the p th percentile value is such that p % of the data are less than this value and (100-p) % are higher), i.e. 3 mg.l -1 , was chosen to be representative of the presence of the turbid plume. 2. Considering this threshold, a region of plume presence (ROPP) was identified using all available images, corresponding to pixels where at least 5 % of the SPM concentrations were above 3 mg.l -1 . This threshold is discussed later in section 4. 3. SPM concentration images with less than 80 % of valid pixels on the ROPP were discarded from the dataset. 4. A gap-filling method [START_REF] Garcia | Robust smoothing of gridded data in one and higher dimensions with missing values[END_REF], based on the discrete cosine transform and which has already proven its ability to deal with geophysical data analysis [START_REF] Wang | A three-dimensional gap filling method for large geophysical datasets: Application to global satellite soil moisture observations[END_REF] was applied on the remaining images to build cloud-compensated images.

Finally, the dataset used for the present study is composed of the MERIS-300m acceptable data available on the area, i.e. 806 daily images from August 14 th 2002 to April 7 th 2012, then covering about 23 % of this period.

Turbid plume metrics extraction and data analysis

The Rhône River turbid plume is characterized by high SPM concentrations near the river mouth, decreasing with the cross-shore distance due to mixing, dilution and particle settling. Considering the shallow waters of the prodelta area and its waves exposure, deposited particles can be resuspended. As ocean color analysis does not permit yet to distinguish if SPM is either resuspended or directly coming from the river, we decided to remove all data where water depths are below 20 m, that we assume to be potentially impacted by wave-induced resuspension.

In order to delineate the turbid plume from clearer waters, the threshold denoting the presence of the turbid plume (3 mg.l -1 ) was applied to the dataset. A routine was developed to identify the turbid plume boundary and extract different metrics information: the area of extension (km 2 ); the southernmost, westernmost and easternmost positions, the centroid location and the schematic shape, called 'skeleton' in the present paper. Mean and maximum near-surface SPM concentrations within the turbid plume were also computed to provide information on seasonal and inter-annual variability of the turbid plume. Seasonal and inter-annual SPM spatial distributions were calculated considering the average of identified plumes.

In order to extract the area of extension of the Rhône River turbid plume, the number of pixels within its boundary was summed and converted to area unit (in km 2 ). Southernmost, westernmost and easternmost points were isolated directly from the boundary and the centroid corresponds to the average position of the pixels constituting the identified plume, regardless of their weight in terms of SPM concentration. The turbid plume schematic shape (skeleton) was computed considering the concentration gradient and the extreme points locations. The identified plume was divided in two parts: a proximal part (high SPM concentration), and a distal part (less concentrated). To identify each part, the same method aiming to isolate the turbid plume from clearer waters was applied inside the plume. The percentile 90 of the SPM concentrations of the plume was used as a threshold to separate distal and proximal parts, for each image. Centroids for both parts were computed and the turbid plume's skeleton was constructed joining (1) the Rhône River mouth, (2) the proximal plume's centroid, (3) the distal plume's centroid, and (4) one of the extreme points (southernmost, westernmost and/or easternmost), depending on the plume shape. If the plume is stretched out seaward, southernmost point is selected. If not, either the westernmost or the easternmost point is selected, depending in which direction it stretches out the most; and if the plume spreads out almost equally to the East as to the West, both easternmost and westernmost points are considered. Turbid plume metrics are shown on Fig. 2f. To highlight Rhône River turbid plume patterns, relationships between plume metrics and meteorological/hydrological forcings were investigated. We chose to study the plume response to the Rhône River water discharge (Q w , in m 3 .s -1 ), the wind velocity (in m.s -1 ) and wind direction (in degrees). The Rhône River discharge data were obtained from the French freshwater office database (http://www.hydro.eaufrance.fr), covering the study time period with daily measurements at the Beaucaire station (58 km upstream the Grand Rhône mouth). Wind velocity and direction at 10 m above sea level were extracted from three meteorological models provided by Météo-France to cover the entire study time period: the ALADIN model for the 19/06/2002 -20/03/2006 period, the ARPEGE model for the 20/03/2006 -24/11/2011 period, and the ARPEGE-HR model for the 24/11/2011 -08/04/2012 period. Each image was associated with a river discharge value, a wind velocity and a wind direction. To consider turbid plume response time to these forcings, water outflow was averaged over the three days before satellite data acquisition. Concerning the wind, pioneer work made by [START_REF] Demarcq | La dynamique superficielle du panache du Rhône d'après l'imagerie infrarouge satellitaire[END_REF] indicated that the plume response time to wind forcing generally varied between few hours to 12 h and can be exceptionally higher, above 24 h for high water discharge values (Q w > 4000 m 3 .s -1 ). Thus, we chose to associate each image to a wind (direction and intensity) averaged over 12 h before satellite acquisition.

Cloud-covered images were discarded from the dataset and MERIS-300m data does not provide daily acquisitions. Therefore the representativeness of the final dataset was investigated: the Rhône River water discharge (Fig. 3a) and wind conditions (Fig. 3b) over 2002-2012 were compared with the conditions associated to each image. Results showed a good agreement of water discharge ranges with a slightly overrepresentation of dry conditions (+ 5 % for Q w < 1500 m 3 .s -1 ), balanced by an underrepresentation of higher discharge values (-5 % for Q w > 1500 m 3 .s - 1 ). Considering winds, southern and eastern winds are underrepresented, and it is especially the case for strong SE winds (> 10 m.s -1 ). These onshore winds, mainly blowing during the winter period, are known to bring clouds over the study area, which explains the underrepresentation highlighted by the analysis. Fig. 3b clearly identifies two main wind sectors, characteristic of the study area: north-western to northern ones (295° to 15°) and south-eastern ones (80° to 160°), qualified hereafter as "onshore winds". We therefore chose to consider these two sectors to study relationships between plume metrics and meteorological forcings.

Results

Seasonal and interannual variability of the Rhône River turbid plume

The Rhône River turbid plume dynamics is studied through monthly averaged nearsurface SPM concentration maps over the 2002-2012 period (Fig. 4) and the analysis of metrics statistics such as plume area and its maximum and mean SPM concentrations (C mean and C max , respectively -Fig. 5).

Results presented in Fig. 4 show that during November, December and January, the Rhône River turbid plume is the largest and the most concentrated, with averaged concentrations always higher than 15 mg.l -1 near the river mouth (for water depth higher than 20 m) and up to 5 mg.l -1 10 km offshore. Statistics in Fig. 5 show that, during this period of the year, for 25 % of the time, the plume expands over 250 km 2 (over 320 km 2 for December and January) and C mean and C max reach values from 7 to 8 mg.l -1 and 33 to 44 mg.l -1 , respectively. From February to April, the turbid plume covers a smaller area (with median values always below 50 km 2 ) and is less turbid than for the three previous months as averaged concentrations never exceed 14 mg.l -1 near the river mouth (maximum recorded in February). For 25 % of the time, the turbid plume expands on 130 to 160 km 2 , with C mean and C max values ranging from 5 to 6 mg.l -1 and 21 to 24 mg.l -1 , respectively. From May to October, a small (median area ranging from 5 to 32 km 2 ) and diluted plume is observed, featured by low average concentrations: C mean and C max median values range from 4 to 5 mg.l -1 and 6 to 10 mg.l -1 , respectively. This seasonal trend of the turbid plume area, C mean and C max is also observed for the Grand Rhône River discharge. Fig. 5 also shows an interannual variability for the turbid plume area, C mean and C max , which is the strongest for the months where the turbid plume is the largest. This variability can be explained by the variability of the Rhône River discharge and related sediment input, following exactly the same trend. During the summer/autumn period, June is atypical as metrics show unexpected high variability: values of first and third quartiles range from 4 to 200 km 2 and 4.5 to 17 mg.l -1 in terms of area and maximum concentration, respectively. This month is particularly interesting because of its dispersion and that the median value for all these metrics is similar to the one of the month of February for instance, where larger and more concentrated plumes are expected. This dispersion is explained by the strong inter-annual variability of the Rhône River turbid plume metrics, shown on Fig. 6 and Fig. 7. A year-by-year statistical analysis highlights that the years 2007-2008-2010, where large and concentrated plumes were observed, correspond to wet years, with median river discharges higher than 1600 m 3 .s -1 . *outliers are definied to be greater than q 3 + w × (q 3 -q 1 ) or less than q 1 -w × (q 3 -q 1 ), where q 1 and q 3 are first and third quartiles, respectively, and w is the whisker value. *outliers are definied to be greater than q 3 + w × (q 3 -q 1 ) or less than q 1 -w × (q 3 -q 1 ), where q 1 and q 3 are first and third quartiles, respectively, and w is the whisker value.

Turbid plume dynamics: response to hydro-meteorological forcings

As monthly average results showed strong variability, it is interesting to use daily (when available) information of the Rhône River turbid plume metrics to investigate its high frequency variability and to link it to the hydro-meteorological forcings. To this end, the turbid plume metrics were separated in two groups, depending on the prevailing winds direction: offshore winds (Mistral and Tramontane -335° to 15°) and onshore winds (80° to 160°). These two groups represent 71 % of the database.

The positive trend observed at the seasonal scale between water discharge (3 days average) and turbid plume area is confirmed by the high frequency analysis (Fig. 8a). For outflows higher than 1500 m 3 .s -1 , the turbid plume expands on a surface always larger than 10 km 2 , 20 % of the time larger than 100 km 2 and can reach values up to more than 1000 km 2 . For dry conditions (Q w < 1500 m 3 .s -1 ), the turbid plume is often small (44 % under 10 km 2 and 92 % under 100 km 2 ) but can expand very episodically up to 550 km 2 for onshore winds. This same positive trend can be observed on Fig. 9a, linking the circles size (representing plume area) to the colour scale (representing the Grand Rhône River discharge), for offshore winds conditions. For dry conditions (Q w < 1500 m 3 .s -1 ), the turbid plume's centroid is located for 95 % of the time at less than 7 km seaward from the Grand Rhône River mouth (43.2662° N latitude). However, for water discharges higher than 1500 m 3 .s -1 , the centroid is for 64 % of the time located at more than 7 km seaward. Concerning meteorological forcing, Fig. 8a shows that wind orientation does not impact the global relationship linking water discharge and turbid plume area. However, it is noticeable that for high river discharges (Q w > 2000 m 3 .s -1 ) and onshore wind conditions, the Rhône turbid plume does not expand as much as it is the case during offshore wind conditions. Indeed, for onshore winds, the area ranges from 200 km 2 to 400 km 2 (with a mean value of ~ 290 km 2 ), whereas in offshore wind conditions, areas range from about 100 km 2 to more than 1000 km 2 (with a mean value of ~ 400 km 2 ). Fig. 8b highlights that onshore winds (in red) tend to confine the turbid plume along the coast as in these conditions, it never expands further than 26 km offshore (74 % below 10 km) whereas Mistral and Tramontane winds stretch the plume seaward up to the limits of the study area. Wind velocity does not seem to affect directly the southern extension of the turbid plume during offshore wind conditions, unlike for onshore wind ones: the stronger the wind is, the more alongshore the plume is. 2, confirms this observation. The skeleton representation of the turbid plume (Fig. 9b) also indicates that, for onshore winds, the Grand Rhône River turbid plume is located at the north of the study area, with almost isotropic orientations. For Mistral and Tramontane winds, the higher the discharge rate, the more defined the turbid plume is. Indeed, for low river discharges (Q w < 1500 m 3 .s -1 ), the plume orientation is almost isotropic whereas for higher discharge rates (Q w ≥ 1500 m 3 .s -1 ), the orientation is better defined, showing a southwestward to westward deflection of the distal part of the plume, giving to the turbid plume a "coma shape". A positive linear relationship was also derived for mean and maximum concentrations detected by the satellite as a function of the Rhône River discharge (Fig 10). As an example, maximum concentrations range from 5 mg.l -1 to 35 mg.l -1 for average discharge conditions while they reach 35 mg.l -1 to 50 mg.l -1 for river discharges larger than 3000 m 3 .s -1 . Given the correlation coefficients, the river discharge explains 54 % and 72 % of Rhône River plume C mean and C max , respectively. Similarly to the turbid plume area dynamics, wind orientation does not dominantly drive the SPM concentration variations, as offshore/onshore scatters are not associated with two distinct populations. However, it is remarkable that mean and maximum concentrations during onshore wind events are significantly lower than during offshore events, and especially for low to moderate river discharges (lower than the mean river discharge). Nevertheless, these differences should be considered with attention as they could be caused by the exclusion of the very close coastal domain (i.e. where depth are lower than 20 m), where the plume could main develop in case of strong onshore winds. These differences are discussed in section 4.

Discussion

Database constitution

Flood conditions are responsible for 80 % of the total Rhône River particulate inputs in the GoL [START_REF] Roditis | Dynamiques fluviales et milieux de sédimentation du Rhône à l'amont immédiat de son delta[END_REF][START_REF] Rolland | Transfert des radionucléides artificiels par voie fluviale: conséquences sur les stocks sédimentaires rhodaniens et les exports vers la Méditerranée[END_REF] and their impact is thus of critical importance in the study of sediment dynamics. However, it is difficult to collect valuable satellite data during these events as they are often related to important rainfall, and thus cloudy weather conditions. The 806 images constituting the MERIS-300m dataset in this area cover 23 % of the 08/14/2002 to 04/07/2012 period and the missing information is almost equally divided between dry (Q w ≤ 1500 m 3 .s -1 ) and other (Q w > 1500 m 3 .s -1 ) conditions, as shown in the left illustration of Fig. 3a. Nevertheless, the study of the dataset representativeness also reveals that, for high river discharges (e.g. Q w ≥ 2500 m 3 .s -1 ), the higher is the discharge, the less it is represented in the dataset: 22 %, 73 % and 81 % of 2500 < Q w ≤ 3500 m 3 .s -1 , 3500 < Q w ≤ 4500 m 3 .s -1 and Q w > 4500 m 3 .s -1 events are missing, respectively (right part of Fig. 3a). Concerning wind conditions, the MERIS database is shown to underrepresent meteorological conditions associated to onshore winds (Fig. 3b), with 46 % of missing values for winds of velocity higher than 4 m.s -1 (representing 13% of the events from the full dataset and only 7% in the MERIS database used for analysis). This underrepresentation is explained by the advection of clouds over the study area by these winds. Despite these missing data, the database is sufficient to provide a good overview of the Rhône River turbid plume and explore its response to different hydro-meteorological conditions.

Threshold for plume detection

Choosing a concentration threshold for detecting turbid plume from ambient background concentration is highly challenging. Field measurements show that marine, offshore GoL concentrations range from 0.7 to 2.2 mgl.l -1 . We chose, in this study, to derive the threshold from a statistical analysis of near surface SPM concentration maps. We first assumed that, in the Gulf of Lions, the presence of high turbidity at the surface is exceptional, and located either close to the coast and in the neighboorhood of river mouths. Next, for embedded areas centered on the rhone river mouth (Fig. 11a), we estimated the percentile 95 of the SPM concentrations observed by the satellite in each area and for all images available. We observed that this value follows a plateau around 3 mg.l -1 for areas larger than 5000 km 2 , and suddenly increased towards the smallest areas (i.e. below 5000 km 2 ), close to the rhone mouth (Fig. 11b). Hence we considered that this plateau corresponds to the ambient background concentration observed by the satellite in the GoL. We also tested the sensitivity of our analysis to this thereshold, and compared the detected areas for threshold values of 3 mg.l -1 and 4 mg.l -1 (Fig. 12). As expected, a larger threshold results in the detection of smaller plume areas (approximately 23 % lower) but this shift follows a linear trend, and hence does not impact the global analysis and the correlation or relationship between metrics and forcings. The same threshold value of 3 mg.l -1 was used by [START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF] in the study of the Adour River turbid plume metrics, determined from field data collected in 2007 [START_REF] Petus | Estimating turbidity and total suspended matter in the Adour River plume (South Bay of Biscay) using MODIS 250-m imagery[END_REF]. In a second step, pixels where water depth is below 20 m were discarded from the analysis. This mask was applied because for these water depths, numerical modelling showed that bottom shear stress (BSS) values are higher than 0.2 N.m -2 for at least 10 % of the time [START_REF] Dufois | Wave-and current-induced bottom shear stress distribution in the Gulf of Lions[END_REF]. This value is close to the critical shear stress of 0.35 N.m -2 used in the simulations of [START_REF] Dufois | Impact of winter storms on sediment erosion in the Rhone River prodelta and fate of sediment in the Gulf of Lions (North Western Mediterranean Sea)[END_REF] for muddy sediments. Therefore, bottom resuspension in shallow areas can be significant as well as their contribution to the surface turbidity signal. This could bias the surface plume detection, which justifies why a 20m-mask was applied to the ocean color data, and visually checked. This choice induces a systematic underestimation of the plume extension area, and for very dry conditions, the plume could potentially not be detected at all. Indeed, for Mistral and Tramontane wind conditions, an area of 15 km 2 in front of the Grand Rhône River mouth is not taken into account for the plume metrics calculation. This value is higher for onshore wind conditions as the flagged area is of about 230 km 2 (along the coast from Grand Rhône River mouth to Petit Rhône River mouth). A sensitivity study to the bathymetric flag (20 m) was performed over onshore-winds associated SPM concentration maps. Results (Fig. 13) revealed that values are biased (but with a mean bias of 46 km 2 ) and significantly impact metrics for plume covering areas smaller than 50 km 2 (underestimating plume area from 20 % to 100 % for the smallest plume area). These differences in metrics are observed both in presence of low wind, i.e. without waves, or high wind, i.e. potentially with waves resuspending bottom sediments, up to the sea surface [START_REF] Dufois | Impact of winter storms on sediment erosion in the Rhone River prodelta and fate of sediment in the Gulf of Lions (North Western Mediterranean Sea)[END_REF]. In the same way, mean and maximum SPM concentrations for onshore wind conditions can be underestimated of about 1 mg.l -1 and 5 mg.l -1 ,

Turbid plume response to hydro-meteorological forcings

The day-by-day analysis confirms the positive correlation between the Rhône River water discharge and the extension of the turbid plume area. The dispersal scaling relation proposed by [START_REF] Warrick | Dispersal scaling from the world's rivers[END_REF] (P = cA b where P is plume area, A the watershed area and b and c regionally-based constants) was tested on our dataset on the Rhone river plume. Considering a decennial flood value of the Rhône River (about 7600 m 3 .s -1 ), the estimated plume area for the relationship described in this paper led to a 0.63 coefficient, in the range of values (0.68 +/-0.04) found by [START_REF] Warrick | Dispersal scaling from the world's rivers[END_REF] in their study comparing different rivers of the world for "special events". Onshore and offshore wind-generated plumes behave the same in response to the river outflow (red triangles are predominantly grouped with blue circles on Fig. 8a). The three red triangles isolated from the point cloud (Q w ~ 1100 m 3 .s -1 and area greater than 200 km 2 ) correspond to one particular event (February, the 12 th , 13 th and 15 th 2008) and not to three different periods. This event can also explain why the averaged area for Q w < 1500 m 3 .s -1 and onshore wind conditions seems high compared to its corresponding value for offshore wind conditions (Table 2). Indeed, discarding these three days from the average calculation, the averaged area would drop down to 25 km 2 . Considering the systematic under-estimation of turbid plume area by the wave mask, plume areas associated to onshore winds would be proportionally more underestimated as they are more often located along the coast (Fig. 9a). From Fig. 8a, it is noticeable that, even if the two populations follow the same trend, onshore wind associated plume areas rise to a level much below the one for Mistral and Tramontane associated ones. For onshore wind conditions, the turbid plume can expand up to 400 km 2 (with a potential maximum underestimation of about 200 km 2 due to wave-mask) while, for offshore wind conditions, it can expand up to more than 1000 km 2 . It is also visible on the centroid representation of Fig. 9a, where the biggest circles, corresponding to the largest plumes, are smaller for onshore winds conditions than for Mistral and Tramontane ones.

N/NW winds strongly drive the surface circulation towards the south, which explains why the plume extends southern during these conditions. On the contrary, onshore wind drive the surface circulation Northward, and explain the location of the plume close to the coast. For Mistral and Tramontane wind conditions, when turbid plumes are stretched to the south, they are globally deflected to the west under the Coriolis effect and the influence of the Northern Current, flowing westward along the shelf slope. Fig. 9b highlights a different orientation of the plume either in the case of Mistral winds (N) or Tramontane ones (NW): first ones tend to promote a southwestward expansion and the latter a southward one. This is particularly visible during wet conditions (Q w > 1500 m 3 .s -1 ), where turbid plume's centroids associated to northwestern winds are located at the west of northern winds associated ones. Skeletons, representative of the plume's shape, also confirm this observation. This observation agrees with the [START_REF] Ekman | On the influence of the earth's rotation on ocean currents[END_REF] theory on the relationship between wind and surface currents and [START_REF] Demarcq | La dynamique superficielle du panache du Rhône d'après l'imagerie infrarouge satellitaire[END_REF] observations, establishing that the Rhône River plume is deflected by about 50° on the right side of the wind.

Using metrics for comparing coastal system dynamics

Metrics are useful proxies for investigating large scale and dynamic coastal structures such as turbid plumes. They provide time series of key features such as plume area or centroid position, which can be confronted to forcings as in the present study to examine the short term to long term plume dynamics (i.e. in the future, the impact of global change). Metrics provide also valuable information when comparing observations and numerical model results: identical parameters are then calculated from the different datasets and can be dynamically compared, complementary to raw day-by-day satellite data analysis. Finally these metrics can be used to compare the dynamics of different coastal systems, and evaluate and scale their sensitivity to driving forcings such as tide, wind, or discharge (suspended solids or water fluxes).

Studies on river plume dynamics using metrics are still scarce but data could be found for different systems. Here, we propose a comparison of plume metrics data for the Rhône River (present study), the Mississippi River [START_REF] Walker | Satellite-based assessment of the Mississippi River discharge plume's spatial structure and temporal variability[END_REF], the Ebro River [START_REF] Fernández-Nóvoa | Analysis of the influence of river discharge and wind on the Ebro turbid plume using MODIS-Aqua and MODIS-Terra data[END_REF] and the Adour River [START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF]. Primarily, it is appropriate to mention that each group of authors did use a different criterion to discriminate the turbid plume from the surrounding waters. For the Rhône River, it is described in the section 2.2 of this paper. Concerning the Mississippi River, [START_REF] Walker | Satellite-based assessment of the Mississippi River discharge plume's spatial structure and temporal variability[END_REF] decided that the turbid plume is defined as areas where SPM concentration is above 10 mg.l -1 but no further information was found on this choice. [START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF]chose to delineate the turbid plume from the surrounding waters with a SPM concentration threshold of 3 mg.l -1 on behalf of in situ measurements carried out during a field campaign [START_REF] Petus | Estimating turbidity and total suspended matter in the Adour River plume (South Bay of Biscay) using MODIS 250-m imagery[END_REF][START_REF] Fernández-Nóvoa | Analysis of the influence of river discharge and wind on the Ebro turbid plume using MODIS-Aqua and MODIS-Terra data[END_REF] tested a large range of turbid threshold values and selected the one with the maximum Pearson correlation coefficient between plume area and river discharge, meaning that they assumed a priori that the plume area is linearly correlated to the river discharge. Fig. 14 shows the superposition of the responses of the turbid plume area regarding the water and solid river discharges for the Rhône River, the Mississippi River, the Ebro River and the Adour River. The Ebro, the Rhone and the Adour systems are characterized by similar river forcing and turbid plumes can be ordered as follow: for a given discharge, Rhone plume areas are smaller than the Adour plumes, which are smaller than the Ebro plumes. The Mississippi system differs from the 3 others by river discharges about one order of magnitude larger than the three others, both in term of water discharge and solid discharge. However, the Mississippi plume dynamics ranges between the Rhone and the Adour dynamics once extrapolated to the Mississippi discharges (even if these extrapolations must be considered with care as they are not representative of the real dynamics of the Adour and Rhone systems).

Rhône River and Mississippi River

These systems have similarities: they are both wind-driven [START_REF] Murray | Observations on wind, tidal, and density-driven currents in the vicinity of the Mississippi River delta[END_REF][START_REF] Schroeder | Circulation patterns inferred from the movement of detached buoys in the eastern Gulf of Mexico[END_REF], subjected to a westward geostrophic flow [START_REF] Schroeder | Circulation patterns inferred from the movement of detached buoys in the eastern Gulf of Mexico[END_REF] and to weak tidal currents in general, with an average tidal range of about 30 cm [START_REF] Murray | Observations on wind, tidal, and density-driven currents in the vicinity of the Mississippi River delta[END_REF]. However, the vorticity within the near Mississippi Delta region could be stronger, with the presence of two recirculation cells on both side of the delta [START_REF] Cochrane | Low-frequency circulation on the Texas-Louisiana continental shelf[END_REF][START_REF] Rouse | Circulation observations in the Louisiana bight using LANDSAT imagery[END_REF], which may explain why the Mississippi turbid plume area is larger than the Rhone turbid plume area once extrapolated to river discharges above 10000 m 3 .s -1 .

Rhone River and Adour River

The relationship between the river discharge and the solid discharge for these two systems is significantly different, especially for river discharges lower than 1000 m 3 .s -1 , when the Adour River delivers about one order of magnitude more suspended sediments to its delta than the Rhône River (Fig. 15). Hence sediment supply may explain the difference in plume patterns. Moreover, the Adour River turbid plume is significantly larger than the Rhone River turbid plume for similar solid discharges (Fig. 14b). These differences could be explained by the hydrodynamic regime of these two systems. Contrary to the Rhône River, the Adour estuary is a macro-tidal system with a tidal range varying between 2 m and 5 m [START_REF] Stoichev | Dynamics of mercury species in surface sediments of a macrotidal estuarinecoastal system (Adour River, Bay of Biscay)[END_REF]. [START_REF] Petus | Monitoring spatio-temporal variability of the Adour River turbid plume (Bay of Biscay, France) with MODIS 250-m imagery[END_REF] suggested that the tide affects the river plume area, and hence this forcing would contribute to disperse the plume on larger areas for similar river discharge.

Rhone River and Ebro River

Different plume patterns are observed for these two Mediterranean microtidal coastal systems: the Ebro plume signature is characterized by areas about one order of magnitude larger than the Rhone river plumes, for similar discharge values. However, the Ebro system is located close to the Rhone Delta, controlled by similar wind regimes [START_REF] Mestres | Modelling of the Ebro River plume. Validation with field observations[END_REF] and driven by the Northern Current, which is flowing with maximum speeds of several 10 cm.s -1 about 40 km offshore the Rhône River mouth [START_REF] Millot | Circulation in the Mediterranean sea, The Mediterranean Sea[END_REF] and from 0.1 to 0.5 m.s -1 , also about 40 km off the Ebro River mouth [START_REF] Font | Marine circulation along the Ebro continental margin[END_REF]. These similarities may discard the hydrodynamic regime to be responsible of these different patterns. Also the Ebro plume area is one order of magnitude larger than the Adour plume, while the Adour solid discharge is significantly larger than Ebro. We assume that the differences could be induced by the methodological approach applied to delineate the turbid plume boundaries.

Finally, variability in suspended sediment behaviour driven by flocculation processes (i.e. the settling velocity), may contribute to induce differences in turbid plume dynamics. Unfortunately no data related to suspended sediment (floc) dynamics are currently available on the investigated systems. However this parameter should be further studied for a global understanding of coastal dynamics.

The methodology applied to detect turbid plumes and calculate metrics is certainly a source of uncertainty when comparing systems. As an example, for this study, a smaller threshold applied to delineate the Rhône River turbid plume (e.g. 2.5 mg.l -1 ) would lightly shift the plume area results to larger values and decrease the differences observed on Fig. 14. The method applied in the present study, based on statistics to assess the SPM concentration "noise" of the background waters, does not rely on a priori assumptions on the plume response to forcing, and is only based on ocean color data, then applicable to any coastal systems. The only requirement is to have sufficient cloud-free satellite images over the area.

Conclusion

An automated processing of ocean color satellite data was developed to isolate the Grand Rhône River turbid plume by a 3 mg.l -1 SPM concentration threshold and extract different metrics such as its area of extension, its south-east-westernmost points, its skeleton (proxy of its shape) or its centroid. This tool was applied to the MERIS-300m database covering the 2002-2012 period and more than 800 images were processed. Despite the under-representation of strong SE wind conditions, often related to high river discharges, these images allowed monitoring the area of influence of the Grand Rhône River turbid plume. Results highlighted its seasonal and strong interannual variability, correlated to the Rhône River water discharge. Dayby-day analysis strengthened this correlation and demonstrated that wind orientation also plays a major role on turbid plume's orientation. Offshore winds tend to enhance the current-induced westward plume deflection whereas onshore winds restrain its southern extension and confine it to the coast. The metrics calculated in this study, and especially the centroid and skeleton extraction, are innovative tools to investigate river plume dynamics in various coastal systems. They are also essential in the validation of 3D hydrodynamic and sediment transport models. Remote sensing, autonomous platforms such as gliders, providing complementary observations during flood events [START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF][START_REF] Many | Particle assemblage characterization in the Rhone river ROFI[END_REF], and numerical models will be next combined to estimate sediment fluxes from the Rhône River mouth to the shelf.

Figure 1 :

 1 Figure 1: RGB (Landsat 8 OLI image recorded on 23 rd february 2014) map of the study area and geolocalisation of Petit Rhône and Grand Rhône river mouths. Northern Current (N.C.), Mistral (M) and Tramontane (T) winds are indicated by black arrows on the GoL map on the top-right illustration.

Figure 2 :

 2 Figure 2: Summary of image selection and metrics extraction data processing: example of January 8th 2012. MERIS available data, summarized in (a), is used to determine the ROPP, which corresponds to the area where, for at least 5 % of all images, pixel value is over the SPM concentration threshold of 3 mg.l -1 (yellow domain in -b). For selected images (where at least 80 % of the data is available on the ROPP) containing flagged pixels (white ones in -c), a gap-filling method is applied (d). To avoid wave-induced resuspension and thus better identify the Grand Rhône turbid plume, pixels where water depths are below 20 m are flagged (grey ones in -e). Metrics are finally estimated considering a SPM concentration threshold of 3 mg.l -1 . (f): Plume centroid is represented here by a black dot, plume extension by the solid black line (area = 655 km 2 ), southernmost, westernmost and easternmost points are the white dots on the turbid plume's boundary, proximal and distal centroids are represented by the white dots inside the boundary. The skeleton is represented by the white line.

Figure 3 :

 3 Figure 3: Meteorological conditions representativeness of the MERIS dataset used in this study for (a) Grand Rhône River discharges (occurrences in %) and (b) averaged wind directions and intensities over 2002-2012 (left partdaily average)

Figure 4 :Figure 5 :

 45 Figure 4: Monthly averaged plume SPM concentration maps over the 2002-2012 period. SPM concentrations values corresponding to water depths lower than 20 m are flagged (grey color). m and n are months and the number of images used for the averaging, respectively.

Figure 6 :Figure 7 :

 67 Figure 6: Interannual variability of Rhône River turbid plume for the month of June. SPM concentrations values corresponding to water depths lower than 20 m are flagged (grey color). y and n are years and number of images used for the avearging, respectively.

Figure 8 :

 8 Figure 8: Rhône River turbid plume response to (a) water discharge and (b) prevailing winds. Blue and red colours are associated to Mistral + Tramontane winds (295° to 15°) and onshore winds (80° to 160°), respectively. The solid black line in (b) represents the southern boarder of the study area.

Fig. 9a also

  Fig.9aalso shows the confining of the turbid plume during onshore wind conditions as centroids are clearly located northern of the study area compared to Mistral and Tramontane winds. The mean offshore distance of plume centroid for various ranges of water discharges, reported on Table2, confirms this observation. The skeleton representation of the turbid plume (Fig.9b) also indicates that, for onshore winds, the Grand Rhône River turbid plume is located at the north of the study area, with almost isotropic orientations. For Mistral and Tramontane winds, the higher the discharge rate, the more defined the turbid plume is. Indeed, for low river discharges (Q w < 1500 m 3 .s -1 ), the plume orientation is almost isotropic whereas for higher

Figure 9 :

 9 Figure9: Rhône River turbid plume's centroids (circles) and skeletons (solid coloured lines) for Mistral, Tramontane (NW to N) and Onshore winds (E to SE). The more extended is the plume, the bigger is the circle representing its centroid. Colour is associated to (a) river discharge (m 3 .s -1 ) or (b) northern winds direction (degrees). For each illustration, a wind rose reminds the wind orientation associated. In (b), dry conditions (Q w < 1500 m 3 .s -1 ) are separated from wet ones (Q w > 1500 m 3 .s -1 ) for more visibility (indicated on illustrations). Smallest points (areas < 100 km 2 ) are blue to dark blue.

Figure 10 :

 10 Figure 10: Correlations between Rhône River discharge and turbid plume (a) mean SPM concentration (C mean ) and (b) maximum SPM concentration (C max ).

Figure 11 :

 11 Figure 11: (a) embedded areas tested for the statistical analysis leading to the choice of the threshold value for turbid plume extraction and (b) associated pc95 values of SPM concentration (mg.l -1 ) regarding each area. Scatters color corresponds to the color of embedded areas.

Figure 12 :

 12 Figure 12: Sensitivity of turbid plume area extraction to threshold selection for values of 3 mg.l -1 and 4 mg.l -1 .

  this bias into account would bring these values closer to the offshore wind associated values and of regression lines of Fig.10.

Figure 13 :

 13 Figure 13: Turbid plume area underestimation for onshore winds conditions. Color scale shows (a) flagged area proportion (%) and (b) wind velocity (m.s -1 ). The solid black line represents the 1:1 line. Embedded graphs show a zoom over the 0-100 km 2 range.

Figure 14 :Figure 15 :

 1415 Figure 14: Superposition of plume area to river (a) water discharge and (b) solid discharge relationships for four different systems: the Rhône River (red dots -present study), the Mississippi River (black stars -(Walker, 1996)), the Ebro River (green triangles -(Fernández-Nóvoa et al., 2015)) and the Adour River (Petus et al., 2014). Solid lines represent the different relationships and dashed lines show their extrapolation to higher discharge rates than commonly observed.

Table 1 :

 1 Number of matchups (N), ε r (%), RMSE (mg.l -1 ) and R 2 coefficients for tested algorithms. The selected product(Doeffer & Schiller, 2007) is highlighted in green.

	Product	Algorithm	ε r (

%) RMSE SPM ≤ 10 (mg.l -1 ) RMSE SPM > 10 (mg.l -1 ) R 2 (%) N

	MEGS8.1	TSM1 (Rrs 665 )	41.8	2.86	6.69	76.5 42
	MEGS8.1 TSM2 (Rrs 665 & Rrs 620 )	46.8	2.83	5.84	80.8 43
	MEGS8.1	TSM3 (Rrs 665 )	40.1	1.61	7.71	73.7 42
	MEGS8.1	Standard (NN)	40.0	2.56	6.45	79.4 46
	SAABIO	TSM1 (Rrs 665 )	47.9	3.37	9.18	62.6 42
	SAABIO TSM2 (Rrs 665 & Rrs 620 )	53.1	3.39	8.68	65.8 43
	SAABIO	TSM3 (Rrs 665 )	37.0	2.02	6.52	79.3 42

Table 2 :

 2 Mean area (km 2 ) and centroid offshore distance (km) of the Grand Rhône River turbid plume.

		Mistral and Tramontane winds	Onshore winds
		mean area	mean offshore distance	mean area	mean offshore distance
		(km 2 )	of the centroid (km)	(km 2 )	of the centroid (km)
	Q w < 1500 m 3 .s -1	28.6	3.4	48.4	2.9
	1500 ≤ Q w < 3000 m 3 .s -1	234.4	9.5	178.2	4.9
	Q w ≥ 3000 m 3 .s -1	649.0	15.5	376.3	3.5

Acknowledgements

This work has been done as part of the AMORAD project and received a state aid managed by the Agence Nationale de la Recherche (ANR) under the program Investissements d'avenir (AMORAD-ANR-11-RSNR-0002). We also thank the two anonymous reviewers for their comments and suggestions that significantly improved the manuscript.


Luc Le Magoarou

Taha Yassine

Stéphane Paquelet

Matthieu Crussière

DEEP LEARNING FOR LOCATION BASED BEAMFORMING WITH NLOS CHANNELS

Massive MIMO systems are highly efficient but critically rely on accurate channel state information (CSI) at the base station in order to determine appropriate precoders. CSI acquisition requires sending pilot symbols which induce an important overhead. In this paper, a method whose objective is to determine an appropriate precoder from the knowledge of the user's location only is proposed. Such a way to determine precoders is known as location based beamforming. It allows to reduce or even eliminate the need for pilot symbols, depending on how the location is obtained. the proposed method learns a direct mapping from location to precoder in a supervised way. It involves a neural network with a specific structure based on random Fourier features allowing to learn functions containing high spatial frequencies. It is assessed empirically and yields promising results on realistic synthetic channels. As opposed to previously proposed methods, it allows to handle both line-of-sight (LOS) and non-line-of-sight (NLOS) channels.

INTRODUCTION

Machine learning (ML) techniques have been applied successfully to wireless communications in recent years (see [START_REF] Timothy | An introduction to deep learning for the physical layer[END_REF][START_REF] Wang | Deep learning for wireless physical layer: Opportunities and challenges[END_REF] for exhaustive surveys). In particular, channel estimation and beamforming in the context of massive MIMO systems have benefited from ML [START_REF] Ding | Dictionary learning-based sparse channel representation and estimation for fdd massive mimo systems[END_REF][START_REF] Wei | Knowledge-aided deep learning for beamspace channel estimation in millimeter-wave massive mimo systems[END_REF][START_REF] Le | Online unsupervised deep unfolding for massive mimo channel estimation[END_REF]. Although such approaches are very promising, their main drawback is that most need pilot symbols to be sent in order to estimate channel state information (CSI) or beam training procedures [START_REF] Barati | Initial access in millimeter wave cellular systems[END_REF]. Both approaches induce a consequent overhead and limit overall system efficiency.

Is CSI or beam training absolutely needed for a massive MIMO system to operate? Not necessarily, since appropriate precoders could be chosen at the base station based on other sources of information about users. In particular, it has been proposed to determine precoders based on estimated user locations, giving rise to the so-called location based beamforming (LBB) [START_REF] Maiberger | Location based beamforming[END_REF][START_REF] Kela | Location based beamforming in 5g ultra-dense networks[END_REF][START_REF] Abdelreheem | Location-based millimeter wave multi-level beamforming using compressive sensing[END_REF]. For such approaches, the location of users is directly used to determine precoders. This allows to spare a lot of resources for communications and enhances the physical layer security [START_REF] Yan | Location-based beamforming for enhancing secrecy in rician wiretap channels[END_REF]. However, the main drawback of LBB with respect to CSI based approaches is that it assumes existence of a line of sight (LOS) propagation path on which the precoder is based. This LOS path This work has been partly funded by the European Commission through the H2020 project Hexa-X (Grant Agreement no. 101015956).

may not exist in some scenarios, thus limiting the applicability of LBB. Nevertheless, the central element of LBB is the mapping that associates each location to a suitable precoder. This mapping can be learned by a neural network, in which case it is not limited to LOS channels. Contributions. In this paper, a supervised learning method allowing to map the user location to an appropriate precoder is proposed. It requires a labeled database containing user locations and the associated channels in order to train a neural network realizing the mapping. The obtained precoders handle with similar success LOS and NLOS channels. The neural network has a specific structure based on random Fourier features (RFF) [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] that allows it to learn functions containing high spatial frequencies. This is beneficial to the task at hand because of the very fast variations of channel coefficients with respect to the user location due to fading.

Note that the location determination step is not studied in this paper, the location being considered known at the base station. This makes the proposed method totally independent of the specific location determination method that is used. In practice, location determination can be done either with help of a global navigation satellite system (GNSS), using a radar, with video cameras looking at the scene, or by sending pilot symbols (much less than for CSI acquisition) in order to determine directions of arrival. Related work. Location based beamforming has been proposed a decade ago [START_REF] Maiberger | Location based beamforming[END_REF]. In its original formulation, it relies only on the main angle of departure (estimated by sending pilots) in order to synthesize an estimated single path channel which is then used to determine a precoder (without any learning). Such an approach is inherently limited to LOS channels, or at least channels comprising a dominant path. LBB Methods that were subsequently proposed [START_REF] Yan | Location-based beamforming for enhancing secrecy in rician wiretap channels[END_REF][START_REF] Kela | Location based beamforming in 5g ultra-dense networks[END_REF][START_REF] Abdelreheem | Location-based millimeter wave multi-level beamforming using compressive sensing[END_REF] introduce refinements to the original idea but do not get rid of this inherent limitation. The method proposed in this paper overcomes this limitation to LOS channels since it directly learns the location/precoder mapping, but as a counterpart, it requires a training phase to calibrate this mapping. The proposed approach pertains to integrated sensing and communication (ISAC) [START_REF] Liu | Integrated sensing and communications: Towards dual-functional wireless networks for 6g and beyond[END_REF], and more specifically to sensing-assisted communication which aims at leveraging sensors at the base station to enhance communication [START_REF] Ali | Leveraging sensing at the infrastructure for mmwave communication[END_REF].

RFFs [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] are a crucial element of the proposed neural network. Their usage in neural networks stems from the impossibility to learn mappings containing high frequencies with classical neural networks such as the multilayer perceptron (MLP) [START_REF] Tancik | Fourier features let networks learn high frequency functions in low dimensional domains[END_REF]. Coordinate based neural networks whose objective is to learn a mapping starting from a very low-dimensional space particularly benefit from RFFs [START_REF] Tancik | Fourier features let networks learn high frequency functions in low dimensional domains[END_REF]. RFF based neural networks have been applied successfully to 3D scene reconstruction [START_REF] Mildenhall | Nerf: Representing scenes as neural radiance fields for view synthesis[END_REF], where the network learns the mapping linking the 3D location and orientation to the radiance and density (which can be seen as the light channel).

PROBLEM FORMULATION

The method proposed in this paper applies to a wide variety of multi-user massive multiple input multiple output (massive MIMO) wideband systems [START_REF] Rusek | Scaling up mimo: Opportunities and challenges with very large arrays[END_REF][START_REF] Erik G Larsson | Massive mimo for next generation wireless systems[END_REF][START_REF] Lu | An overview of massive mimo: Benefits and challenges[END_REF], operating indifferently in time division duplex (TDD) or frequency division duplex (FDD), where the antennas at the base station are indifferently colocated or not (in which case it is a distributed MIMO system). Let us consider A base station antennas and a single subcarrier, and denote h ∈ C A the downlink channel vector between the base station and any given user, and l ∈ R D its location, where D can be two or three, depending on whether or not the elevation of users is considered relevant. Moreover, the proposed method can be directly transposed to systems comprising several subcarriers, but is presented here for a single subcarrier for ease of exposition.

Based on a labeled database of N downlink channels associated with the corresponding user locations

{h i ;l i } N i=1 , (1)
the objective in this paper is to build a location based precoding function (or simply precoding function), whose function is to focus energy on the intended user location. It is mathematically defined as

P : R D → C A l → w P(l), (2)
where w is the predicted precoding vector (of unit norm). Performance measure. In order to evaluate precoding functions, the normalized correlation between the precoder w and the channel h is used. It is expressed as

η |w H h| 2 h 2 2 . ( 3
)
It is between zero and one (for a perfect precoder), and is tightly linked to the downlink channel capacity (considering a single user), whose expression is

log(1+η.SNR opt ),
for received signal of the form y = √ P w H hs + n where n ∼ N (0,σ 2 ) is additive noise, s is the sent symbol and P the transmit power. In that setting, SNR opt

P h 2 2 σ 2
is the highest achievable signal to noise ratio (SNR) [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF]. In summary, the correlation is a single number between zero and one allowing to determine the maximum achievable downlink spectral efficiency for any transmit power and noise variance considering a given precoder.

PROPOSED SOLUTION

Neural precoding function. Deep neural networks are known to be universal function approximators [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] and have led to great practical successes [START_REF] Lecun | Deep learning[END_REF]. It is proposed here to implement the precoding function P as a deep neural network. However, it has been shown that classical neural networks known as multilayer perceptrons (MLP) [START_REF] Rumelhart | Learning Internal Representations by Error Propagation[END_REF] are in practice unable to learn functions of high frequency [START_REF] Tancik | Fourier features let networks learn high frequency functions in low dimensional domains[END_REF]. This phenomenon is known as spectral bias [START_REF] Nasim Rahaman | On the spectral bias of neural networks[END_REF][START_REF] Basri | The Convergence Rate of Neural Networks for Learned Functions of Different Frequencies[END_REF] and has been characterized with help of the theory of neural tangent kernels (NTK) [START_REF] Jacot | Neural tangent kernel: Convergence and generalization in neural networks[END_REF].

In order to remedy this fundamental weakness, it has been proposed [START_REF] Tancik | Fourier features let networks learn high frequency functions in low dimensional domains[END_REF] to use random Fourier features [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] to help the neural network building high frequencies. This line is followed in this paper, motivated by the fact that the optimal precoder may vary fast with respect to the user's location, due to fading. This results in a neural network whose first layer is fixed and corresponds to random Fourier features expressed as

γ(l) = [cos(2πBl),sin(2πBl)] T ,
where B ∈ R R×D contains R D-dimensional (spatial) frequencies drawn randomly as

B ∼ N (0,s 2 I),
where the variance parameter s 2 controls the frequency range (the higher s 2 the more high frequencies are likely to appear). It is proposed to concatenate this RFF layer with a classical MLP using fully connected layers and rectified linear units (ReLU) activation functions [START_REF] Nair | Rectified linear units improve restricted boltzmann machines[END_REF] of width M (number of neurons per layer) and depth Q (number of layers). The last layer is of width 2A, which corresponds to the real and imaginary parts of the predicted precoder w stacked together. The structure of the proposed neural network is shown on Fig. 1. Simply put, the role of the RFF layer is to pre-build high frequencies in order to ease learning for the subsequent MLP.

In order to calibrate the weights, the cost function to minimize is expressed as

CF 1- 1 N N i=1 |P(l i ) H h i | 2 h i 2 2 . ( 4
)
It measures the misalignment of the predicted precoders with respect to the training channels and is between 0 (if the precoders are perfectly aligned with the channels) and 1 (if the precoders are orthogonal to the channels). Practical setup. Two distinct phases are to be distinguished for the proposed method to operate in a real system:

• Training phase: At first, the database {h i ;l i } N i=1 has to be built. To do so, the base station can send pilots and estimate channels using any classical method. Moreover, it has to access users' locations, using either a GNSS, video cameras or also pilots. The neural network can be trained with the built database in order to calibrate the precoding function P. • Inference phase: Then, once the neural network is properly trained, no more channel estimation is needed. The base stations only has to estimate users locations and is then able to predict appropriate precoders with the calibrated precoding function.

Note that these two phases are not mutually exclusive, and can be intertwined in an online mode, in order to allow the base station to adapt to slow changes in the environment.

EXPERIMENTS

In this section, the proposed method is empirically assessed and compared to concurrent approaches on realistic multipath synthetic channels. Simulation settings. The experiments are carried out using training multipath channels taken from the DeepMIMO dataset [START_REF] Alkhateeb | DeepMIMO: A generic deep learning dataset for millimeter wave and massive MIMO applications[END_REF]. The 'O1' urban outdoor ray-tracing scenario is chosen (with 5 paths per channel, see [START_REF] Alkhateeb | DeepMIMO: A generic deep learning dataset for millimeter wave and massive MIMO applications[END_REF]Figure 2] for a schematic view of the scenario), with a single base station (BS 16) and a subset of all possible user locations (rows R1000 to R3852). The base station is equipped with a square uniform planar array (UPA) with A = 64 half-wavelength separated antennas at a frequency of 3.5GHz. Implementation details. The method is implemented with help of the PyTorch library [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF], so that gradients are computed automatically. The optimization of the precoding function is done according to the cost given in (4) on a database of N = 10000 gradient descent (minibatches of 100 training channels) using the Adam optimization algorithm [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF] for 50 epochs. The number of RFF is set to R = 1000, and the standard deviation used for their computation is chosen by cross validation and fixed so that 1 s = 50m. The MLP following the RFF computation is of depth Q = 4 and width M = 512. Note that complex weights and inputs are handled classically by stacking the real and imaginary parts so that the neural network treats only real numbers.

Baselines. The proposed method is compared with two baselines:

• A classical LBB approach [START_REF] Maiberger | Location based beamforming[END_REF][START_REF] Kela | Location based beamforming in 5g ultra-dense networks[END_REF], for which the precoder is chosen as the normalized LOS channel in the direction of the user, assuming the azimuth and elevation are perfectly known. This approach would be optimal in terms of correlation for LOS channels comprising a single path. This baseline is here to show the potential of the proposed method with respect to prior art. Fig. 2: CDF of correlation between the channel and the precoder for several precoding methods.

• A deep learning approach using a simple MLP (without RFFs). The architecture of the MLP is exactly the same as the one of the proposed method, except the first layer (RFF computation) that is replaced by a simple fully connected layer. This baseline is here to show the interest of RFFs.

Results. The performance of the three compared methods on a test set of 10000 channels (not used for training) is shown on Fig. 2. The cumulative distribution function (CDF) of the correlation is plotted (with an ideal precoder it would be equal to zero for all values of x < 1 and to one for x = 1). The plot clearly shows the great improvement brought by the proposed method. Indeed, in such an urban environment where a LOS path is not available for every location, the classical LBB approach yields a median correlation of 0.708. This is explained by the fact that the approach is inherently biased towards LOS channels. Using a deep learning approach with a simple neural network (without RFF) yields a slight improvement (median of 0.714). This very simple neural network is prone to too much variance. However, the proposed method using RFFs allows to attain a median of 0.946, which is a lot better than what is attained by the two aforementioned baselines, because of a better bias/variance trade-off.

It is also interesting to look at the spatial distribution of the attained correlation. This is depicted on Fig. 3 (for the classical LBB approach) and Fig. 4 (for the proposed approach). On these figures, the colors (blue to yellow) denote the obtained correlation for 50000 test channels, and the base station location is denoted by a red cross. Note that on both figures, black areas on the left correspond to zero-norm channels. Comparing these figures, it is obvious that the classical LBB approaches (Fig. 3) yield good precoders (yellow) for areas where user are in line of sight with respect to the base station, but very bad precoders elsewhere (blue areas). On the opposite, the proposed approach (Fig. 3) yields good precoders everywhere on the considered area, even in locations that are in non-line of sight with respect to the base station. This clearly shows the benefit of the proposed method, which allows to obtain good precoders in areas where classical LBB approach fails, again owing to the incapacity of handling NLOS channels.

In order to assess the applicability of the proposed method, it is interesting to vary the number of training channels N in order to determine which density of training channels is required. To that aim, the results in terms of average and median correlation are shown as a function of N on Fig. 5. From this figure, it can be seen that, obviously, the greater N the greater the test correlation. However, one can see that the proposed method starts to perform well (on test data) for N > 3000. The area in which users can be located spans a total surface of 20520m 2 , which means that the proposed method requires approximately one training channel every 20520 3000 ≈ 7m 2 on average to perform well. Such an empirical result is important in order to determine the time duration of the training phase in which both channels and locations have to be collected in order to build the database used to train the neural network.

CONCLUSION

In this paper, a location based beamforming method was proposed. It is based on a neural network comprising RFFs allowing to learn functions of high spatial frequency. The method was empirically assessed on realistic synthetic channels, yielding very promising results. Indeed, it is able to handle both LOS and NLOS channels whereas previously proposed methods perform well only for LOS channels. This shows the capacity of the proposed structure to interpolate well between training points thanks to RFFs which inject prior information on the learned mapping.

In the future, the method could be extended in several ways. First of all, the distribution of frequencies used to build the RFFs could be optimized to further enhance performance. It would also be interesting to aim at predicting not only an appropriate precoder but directly the channel vector. This would amount to summarize the whole propagation environment in the weights of a neural network, and could serve many applications, such as channel mapping [START_REF] Alrabeiah | Deep learning for tdd and fdd massive mimo: Mapping channels in space and frequency[END_REF][START_REF] Magoarou | Similarity-based prediction for channel mapping and user positioning[END_REF] (in space) or radio environment compression.

Fig. 1 :

 1 Fig. 1: Proposed RFF based neural network (with Q = 4).

Fig. 3 :

 3 Fig. 3: Spatial performance of direction based LBB [7, 8]. Only LOS regions are well covered.

Fig. 4 :

 4 Fig. 4: Spatial performance of deep learning based LBB using RFFs (proposed). Both LOS and NLOS regions are well covered.

Fig. 5 :

 5 Fig. 5: Correlation versus number of training channels.