Gaël Many
email: gael.many@univ-perp.fr

François Bourrin

Xavier Durrieu De Madron

Anouck Ody

David Doxaran

Pierre Cauchy

Glider and satellite monitoring of the variability of the suspended particle distribution and size in the Rhône ROFI

Keywords: Nepheloid layers, ROFI, Gulf of Lions, Rhône River, flood, glider

High-resolution description of the variability of nepheloid layers in the Rhône ROFI -Evidence of the Rhône River as a major supplier of the bottom nepheloid layer -Evidence of a storm-induced nearshore break of the plume stratification -Combination of glider and satellite data to characterize the particle size variability

Introduction

River plumes, which are the primary signature of continental inputs of freshwater as well as dissolved and particulate materials, can influence stratification and circulation patterns, sediment and nutrient budgets, biogeochemistry processes and pollution in coastal waters [START_REF] Dagg | Transformation of dissolved and particulate materials on continental shelves influenced by large rivers: plume processes[END_REF][START_REF] Nittrouer | Continental margin sedimentation: From sediment transport to sequence stratigraphy[END_REF][START_REF] Simpson | Introduction to the physical and biological oceanography of shelf seas[END_REF]. They remain poorly documented regarding in-situ monitoring of their spatio-temporal variability, which is strongly affected by the intensity of river discharge and wind conditions [START_REF] Garvine | A dynamical system for classifying buoyant coastal discharges[END_REF].

These forcings also strongly affect the transport of suspended particulate matter (SPM) from the continent to the shelf sea, the specific particle dynamics (aggregation/disaggregation, settling velocity) of which primarily depend on their concentration, size, and density [START_REF] Curran | Settling velocity, effective density, and mass composition of suspended sediment in a coastal bottom boundary layer, Gulf of Lions, France[END_REF]Manning and Schoellamer, 2013).

Over the past decade, different works on river plume dynamics have been conducted with ocean color imagery of surface SPM [START_REF] Thomas | Satellite-measured temporal variability of the Columbia River plume[END_REF][START_REF] Lihan | Satellite-measured temporal and spatial variability of the Tokachi River plume[END_REF][START_REF] Saldías | Seasonal variability of turbid River plumes off central Chile based on high-resolution MODIS 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][START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF]. These authors described the benefits of such observations in the description of the large-scale spatial and temporal variability of river plumes. Nevertheless, satellite data provide information about neither the vertical structure of river plumes (as water stratification, plume thickness, SPM vertical gradient) nor the characteristics of bottom nepheloid layers (BNL). In-situ monitoring throughout the water column is thus essential to assess the link between surface plumes and BNL that plays a major role in sedimentary transport and deposition [START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF].

The recent use of coastal gliders has made it possible to reduce this lack of knowledge in the water column. Castelao et al. (2008aCastelao et al. ( , 2008b) ) first used a coastal glider to describe the dynamic of a river plume in response to upwelling winds. More recently, [START_REF] Saldías | Optics of the offshore Columbia River plume from glider observations and satellite imagery[END_REF] used glider and satellite observations of salinity to characterize the seasonal variability of the Columbia River plume. To our knowledge, however, there are no studies that have characterized the structure and variability of nepheloid layers in ROFI using SPM measurements from both gliders and satellites. A need exists in the high frequency observations of the 3-D temporal variability of SPM concentrations and size in ROFIs during flood events [START_REF] Nittrouer | Continental margin sedimentation: From sediment transport to sequence stratigraphy[END_REF][START_REF] Simpson | Introduction to the physical and biological oceanography of shelf seas[END_REF]. This work brings together a comprehensive set of hydrological, hydrodynamical, and optical measurements collected from a coastal buoy, a satellite and a glider in the Rhône River ROFI.

It aims at: 1) assessing the temporal variability of the Rhône River ROFI nepheloid layers during a flood event and different wind conditions, and 2) determining the temporal variability of the SPM size.

Regional setting

In the semi-enclosed Mediterranean Sea, the inputs of freshwater and sediment are of paramount importance for the sedimentary budgets at both the local and Mediterranean scales.

The Po and Rhône rivers alone account for approximately 30% of the freshwater discharge to the Mediterranean [START_REF] Ludwig | River discharges of water and nutrients to the Mediterranean and Black Sea: Major drivers for ecosystem changes during past and future decades?[END_REF]Sadaoui et al., this issue). Because of the importance of their freshwater inputs and the large extent of their ROFIs, these rivers are key control factors of the particle dynamics for the Gulf of Lions and the Adriatic shelves, respectively.

The Gulf of Lions (NW Mediterranean) is a wide continental shelf mostly influenced by the inputs of particulate matter from the Rhône River, which represents on average between 80 and 95% of the total inputs to the gulf [START_REF] Bourrin | Contribution to the study of coastal rivers and associated prodeltas to sediment supply in the Gulf of Lions (NW Mediterranean Sea)[END_REF][START_REF] Sadaoui | Controls, budgets and variability of riverine sediment fluxes to the Gulf of Lions (NW Mediterranean Sea)[END_REF]. The Rhône River has a mean annual discharge of 1700 m 3 s -1 and typical annual floods with discharge higher than 5000 m 3 s -1 [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]. According to [START_REF] Pont | Evaluation of water fluxes and sediment supply[END_REF], floods deliver 70% of the annual sediment load, highlighting the pivotal role of these events in the delivery of SPM to the coastal area. During floods, freshwater and sediment delivered by the Rhône River form a surface plume and a BNL that can spread over the shelf [START_REF] Aloisi | Origine et rôle du néphéloïde profond dans le transfert des particules au milieu marin. Application au Golfe du Lion[END_REF]. Previous works focused on the remote sensing of the plume extent [START_REF] Arnau | River mouth plume events and their dispersal in the Northwestern Mediterranean Sea[END_REF]Lorthois et al., 2012;[START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF] or used shipborne in situ measurements to characterize the suspended particles in the water column [START_REF] Naudin | River discharge and wind influence upon particulate transfer at the land-ocean interaction: case study of the Rhone River plume[END_REF][START_REF] Thill | Evolution of particle size and concentration in the Rhône river mixing zone: influence of salt flocculation[END_REF][START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF]. The plume is generally advected seaward across the shelf by surface currents driven by continental (i.e., northerly) wind or restricted to the coast during marine (i.e., southeasterly) wind events [START_REF] Broche | Experimental study of the Rhone plume. Part I: physics and dynamics[END_REF]; see also modelling works of [START_REF] Marsaleix | A numerical study of the formation of the Rhône River plume[END_REF][START_REF] Estournel | The Rhone River plume in unsteady conditions: numerical and experimental results[END_REF][START_REF] Ulses | Suspended sediment transport in the Gulf of Lions (NW Mediterranean): Impact of extreme storms and floods[END_REF]. The coastal area generally presents a strong seaward decrease of surface SPM concentrations due to the dilution of the plume and the settling of large particles (such as sand and silt) on the prodelta off the river mouth. This settling is enhanced by the flocculation of fine sediments that enables the settling of "flocs" closer to the coast [START_REF] Naudin | River discharge and wind influence upon particulate transfer at the land-ocean interaction: case study of the Rhone River plume[END_REF][START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF]. Near the seabed, the BNL extends over the inner shelf off the river mouth and presents a seaward decrease of the SPM concentration and size (Durrieu de [START_REF] Durrieu De Madron | Transport de matière particulaire en suspension sur le plateau continental du Golfe du Lion. Situation estivale et hivernale[END_REF][START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF]. On the Rhône prodelta, the BNL is fueled by the settling of large particles from the river plume [START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF]. Resuspension induced by waves, especially during storm events, also contributes to the supply of the BNL (Durrieu de Madron et al., 2008;[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].

Material & Methods

Presentation of the experiment

Time-series measurements of the coastal MESURHO buoy (43°19.130'N; 4°51.980'E, 20 m depth, 2.6 km off the Rhône River mouth) were used to describe the meteorological and oceanic conditions in front of the Rhône River mouth. Additionally, a SLOCUM glider was deployed within the Rhône ROFI during a 14-day period (January, 30 to February, 12 2014).

Seven complete cross-shelf sections (Fig. 1) were carried out from the coast to the shelf edge (40 km offshore -120 m depth).

Nearshore hydrological, hydrodynamical and meteorological observations

Meteorological time-series -Hourly (10 minute burst average) wind speed and direction were measured at the Cap-Couronne Météo-France Station (43°20'.230'N; 5°01.380'E). We derived the wind stress at the sea surface, τ, expressed in N m -2 , following Equation 1:

τ = W f 2 × ρ air × C D Eq. 1
where W f is the wind speed (m s -1 ), ρ air the air density (assumed to be 1.22 kg m -3 ) and C D the drag coefficient of 1.3 × 10 -3 [START_REF] Large | Open ocean momentum flux measurements in moderate to strong winds[END_REF].

Rhône river time-series -The daily Rhône River discharge was measured at the Beaucaire gauging station, 50 km upstream from the river mouth (Compagnie Nationale du Rhônecode V7200015).

MESURHO buoy time-series -Currents and waves were measured using a 600 kHz upwardlooking RDI ADCP fixed at the bottom (27 m) of the MESURHO buoy off the Rhône River mouth [START_REF] Pairaud | MesuRho. Plateforme instrumentée de suivi des paramètres environnementaux à l'embouchure du Rhône[END_REF]. Hourly mean horizontal currents (1 m bin) were estimated from the high frequency horizontal current profiles. Wave directions were estimated by computing near-surface waves orbital, pressure and surface track [START_REF] Teledyne | Waves primer: Wave measurements and the RDI ADCP waves array technique[END_REF].

High-frequency wave measurements were then averaged over 20 min period every 3 h.

ADCP's backscatter index (BI in dB) was derived from the sonar equation (Eq. 2). Because of the sensitivity of acoustic sensors to large particles (as aggregate and zooplankton), BI was only used as relative information concerning the SPM concentration, complementary to optical sensors more sensitive to fine sediments. It can be calculated by adding the echo intensity (EI) received by the instrument, the transmission loss by the sound absorption in seawater (Tl w ), the beam spherical spreading (Tl g ) and the constant A, which incorporates the emitted intensity and the time-and distance-independent coefficients [START_REF] Gartner | Estimating suspended solids concentrations from backscatter intensity measured by acoustic Doppler current profiler in San Francisco Bay, California[END_REF][START_REF] Gostiaux | Extracting meaningful information from uncalibrated backscattered echo intensity data[END_REF].

BI = EI + Tl w + Tl g + A

Eq. 2

Tl w and Tl g depend on the temperature (13°C) and salinity (38) of seawater, the beam angle (20°), the distance to the ensonified SPM volume and the absorption coefficient. We assumed these parameters to be constant in space and time.

Glider CTD and optical observations

The autonomous underwater glider (30 -200 m) used for this work is a Teledyne Webb

Research Slocum G1 [START_REF] Davis | Autonomous buoyancy-driven underwater gliders[END_REF]. Conductivity, temperature and pressure measurements were made by an un-pumped SeaBird 41cp CTD. Conductivity and temperature measurements were corrected for thermal lag effect to avoid salinity bias [START_REF] Garau | Thermal Lag Correction on Slocum CTD Glider Data[END_REF]. Salinity and density were derived following the equation of EOS-80 (UNESCO, 1983). A Wetlabs FLNTU provided turbidity (expressed in nephelometric turbidity units -NTU) based on light backscattering measurements at 700 nm. We used the calibration of [START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF], established during the same deployment and based on in situ SPM gravimetric measurements, to convert NTU in mg L -1 . A Wetlabs BB2FLS provided light backscattering measurements (m -1 sr -1 ) at 532 and 660 nm.

The glider moved in a sawtooth-shaped trajectory (generally with a pitch angle of 25°)

between the surface (1 m depth) and the seabed (2 meters above bottom (mab)). We generally observed an average horizontal speed of 0.4 m s -1 with an average vertical speed of 0.2 m s -1 .

Each section was generally performed in 1.5-2 days. It is thus important to note that a glider section cannot be considered as a snapshot of the coastal sea and can present different processes along the glider track depending on wind and river discharge conditions. When surfacing every 6 dives, the glider took GPS fixes and corrected its navigation toward the next waypoint. Data were first projected along a N-S transect and then interpolated within a 1 m vertical and 100 m horizontal resolution grid.

Optical, hydrodynamical and hydrological derived parameters

Spectral slope estimate -We computed the spectral slope (γ Glider ) of glider particulate light backscattering coefficients (b bp (λ)) at 532, 660 and 700 nm to describe the spatial and temporal variability of the particulate assemblage within the Rhône ROFI [START_REF] Boss | Shape of the particulate beam attenuation spectrum and its inversion to obtain the shape of the particulate size distribution[END_REF].

FLNTU (λ = 700 nm) and BB2FLS (λ = 532 and 660 nm) sensors provided light scattering measurements (β(θ,λ)) at specific angles θ = 142° and θ =124° in the backward direction [START_REF] Sullivan | Coastal thin layer dynamics: consequences to biology and optics[END_REF]. The BB2FLS sensor directly gives output in m -1 sr -1 , while measurements of the FLNTU (factory calibrated in NTU) were converted to extract the volume scattering function thanks to the Beta Scale Factor provided by the manufacturer (i.e., the NTU Scale Factor multiplied by a coefficient of 2.5 × 10 -3 ). We did not correct values of (β(θ,λ)) for loss by absorption due to negligible effects (a few % of the total signal measured).

For each sensor, we estimated the volume scattering function of particles (β p (θ, λ)), by subtracting the molecular scattering of pure water [START_REF] Morel | Optical properties of pure water and pure sea water[END_REF]. The particulate backscattering coefficients (b bp (λ)) (m -1 ) were finally derived following (b bp (λ)) = 2π.X.(β p (θ,λ)), where X is an adjustment factor provided by the manufacturer according to particle type (1.077 for the BB2FLS and 1.132 for the FLNTU).

From [START_REF] Morel | The scattering of light by sea water: Experimental results and theoretical approach[END_REF] and [START_REF] Babin | Light scattering properties of marine particles in coastal and open ocean waters as related to the particle mass concentration[END_REF], the light scattered by suspended marine particles depends on the particle size distribution and is wavelength-dependent. The above authors described this variation as a power law approximation in the case of spherical, homogeneous and non-absorbing particles: b bp (λ)=b bp (λ 0 ).( λ/ λ 0 ) -γ Glider , where λ 0 is the reference wavelength and γ Glider the spectral slope. γ Glider was thus calculated as the slope of the linear relationship between log e (b bp (λ)) and log e (λ), using λ = 532, 660 and 700 nm (i.e., the lower the γ Glider , the higher the proportion of large particles to the total particle assemblage, and reciprocally).

Coastal along-shelf current estimate -The glider "drift", i.e., the difference between the glider's expected surfacing location and the actual new GPS position, was used as an estimate of horizontally depth-averaged current between two surfacings [START_REF] Research | Operations manual Slocum shallow battery glider[END_REF]. To further characterize the current, we estimated geostrophic velocities to determine how much the flow perpendicular to the section was affected by the density field.

The glider's interpolated density field was used to estimate the baroclinic component of the geostrophic velocities. Then, we derived the E-W velocity field (i.e., normal to the gliders' section) by adjusting the depth-averaged geostrophic velocities to the corresponding glider's drift. This fitting was made by estimating the difference between both depth-averaged velocities. The difference was then added to each geostrophic velocity profile to find a final depth-averaged velocity equal to the depth-averaged velocity derived from the glider drift.

This correction includes the barotropic geostrophic and ageostrophic terms, such as currents induced by surface wind stress and bottom friction.

Hydrological derived parameters -The Brunt-Väisälä frequency (N -Eq. 3) was used as an indicator of the stratification. One maximum, indicative of the pycnocline between the colder and fresher plume and the shelf water, appeared on the different sections (Fig. 2). A threshold of N = 0.01 s -1 was used to delimit the pycnocline and defined the lower limit of the plume, which corresponded to the 28.3  0.1 kg m -3 isopycnal (Fig. 2). We used the dimensionless Richardson number (Ri -Eq. 4) to estimate the mixing intensity by comparing the Brunt-Väisälä frequency to the vertical gradient of the horizontal current velocity [START_REF] Pacanowski | Parameterization of vertical mixing in numerical models of tropical oceans[END_REF].

N = √ (-g/ρ) δρ/δz

Eq. 3 Ri = N² / (δu/δz) ² Eq. 4

where g is the acceleration of gravity (m² s -1 ), ρ the potential water density (kg m -3 ) and u the E-W (i.e., along-shelf) current velocity component (m s -1 ) estimated from glider measurements (i.e., geostrophic velocities adjusted to the glider drift). The stratification is stable for N > 0 and vice-versa. The Richardson number, which defines the mixing intensity, is equal to 0.25 in a stable flow and Ri < 0.25 if the turbulence overcomes the density stratification and generates vertical mixing [START_REF] Turner | Buoyancy effects in fluids[END_REF].

Satellite data

Spatial maps of SPM concentrations (250 m resolution) were generated from Moderate Resolution Imaging Spectroradiometer (MODIS) data recorded onboard the Aqua satellite platform. Intermediate Level 1-A (L1A) products were downloaded from the NASA Ocean

Color website (http://oceancolor.gsfc.nasa.gov/cms/). Processing was completed using SeaDAS (version 7.0) software to generate geo-referenced and radiometrically calibrated Level 1-B (L1B) files. L2 remote sensing reflectance products at 645 nm (Rrs(645)) were generated using the l2gen function. The MUMM atmospheric correction [START_REF] Ruddick | Atmospheric correction of SeaWiFS imagery for turbid coastal and inland waters[END_REF] was used because it is the most appropriate to correct MODIS data from aerosol effects in the moderately turbid waters of the Rhône River plume [START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF]. Rrs(645) was converted to SPM mg L -1 applying the regional equation established by [START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF]. The same authors studied the surface plume dynamics and fixed a threshold of 10 mg L -1 to characterize the plume boundary when using this regional equation. Thus, we used a similar value to determine the extent of the surface plume. We derived the particulate backscattering coefficient bbp(λ) at 10 MODIS wavebands between 412 and 678 nm applying the QAA (Quasi-Analytical Algorithm) [START_REF] Lee | Deriving inherent optical properties from water color: a multiband quasi-analytical algorithm for optically deep waters[END_REF](Lee et al., , 2013)). QAA derives the absorption and backscattering coefficients by analytically inverting the spectral remote-sensing reflectance.

The spectral bbp(λ) is modeled assuming that both pure water and particles contribute to the total backscattering coefficient, and assuming that the spectral variations of bbp(λ) follow a power-law function. The spectral slope of the particulate backscattering coefficient, γ MODIS , was therefore obtained by fitting a power-law function to the spectral bbp coefficients.

A filter was applied on γ MODIS to exclude slope values obtained from power-law functions fitted with a R² coefficient lower than 0.9. Excluded values generally corresponded to SPM MODIS > 35 mg L -1 and were located close to the river mouth during peak discharges.

γ MODIS is further used as relative information of the surface spatial variability of the particulate assemblage (i.e., the lower the spectral slope, the higher the abundance of large particles, see 3.4).

Concomitant MODIS and surface glider data of SPM concentrations and backscattering spectral slope estimates were determined using a 6-hour window (3-hours before and after the satellite image was taken, see [START_REF] Saldías | Optics of the offshore Columbia River plume from glider observations and satellite imagery[END_REF]. This window is commonly used. While it increases uncertainty due to the different sampling time, it increases the amount of data used to compare satellite and in situ data.

Results

Meteorological and oceanic conditions

During the experiment, different wind conditions were observed. North-westerly (i.e., continental) winds were measured on January 30, February 2 (Situation 1) and 11 (Situation 3), 2014, with a maximum wind speed of 12 m s -1 . Several south-easterly (i.e., marine) wind events of ~10 m s -1 occurred on February, 5-7 (Situation 2), 8 and 10, 2014, which led to precipitation over the Rhône water catchment (Fig. 3a) and resulted in a rapid increase of the Rhône discharge from 3500 m 3 s -1 (February 5) to 5500 m 3 s -1 (February 11) (Fig. 3b). During these SE wind events, waves mainly came from the south with significant height ranging from 2 to 3.5 m (Fig. 3c). Currents were stronger and generally oriented toward the north (maximum current speed was observed during such an event on February 10). Conversely, currents were weaker and oriented toward the south for NW winds.

ADCP echo intensity showed that SPM concentrations in the water column followed the increase of the Rhône River discharge and wave height conditions with a maximum acoustic intensity at the peak of the flood. During episodic high waves and strong northward currents (February, 10), a drastic decrease of SPM concentrations was observed throughout the water column (Fig. 3f). Dash-line boxes on Figure 3 represent the three situations describe hereafter.

Meteorological and hydrodynamical conditions during these situations are summarized in Table 1. The third situation (February, 10-11) is further qualified as a "back-to-back" storm and flood event due to the marine storm occurring on February, 10 (i.e., marine winds > 10 m s -1 and waves > 3 m) and the subsequent flood event on February, 11. *Due to the "back-to-back" storm and flood events of the third situation both marine and continental wind events are detailed.

Glider cross-shelf CTD measurements

Temperature and salinity sections measured by the glider between January 30 and February 12, 2014 are presented in Figure 4. The results showed the cross-shelf extension of the Rhône River surface plume with fresher (<37.6), colder (<13.5°C) and lighter water (< 28.3 kg m -3 ).

During the first two sections (Fig. 4a-b), the plume was pushed offshore by continental winds.

During the following three sections, under a prevailing S-SE wind (Fig. 4c-e), the plume was pressed against the coast and its thickness increased. The sixth section (Fig. 4f) showed a larger plume, which extended offshore due to continental wind and an increase of the Rhône discharge. Due to the glider direction, we monitored the structure of the plume close to the coast twice within a small time interval (i.e., a few hours during the 10 February 2014).

During this period, the plume separated from the river mouth and deepened to 50 m deep. The seventh section (Fig. 4g) showed the largest plume, which extended up to 45 km offshore, associated with the large discharge of the Rhône River. Finally a notable decrease in temperature was observed near the bottom, forming a distinct layer extending southward with a thickness of 15 m. Nepheloid layers dynamics -Three situations were selected to characterize the SPM dynamics in the ROFI for different river discharges, wind and sea conditions (Fig. 6). The dimensions of the surface plume are derived from both satellite measurements using a threshold of 10 mg L -1 [START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF], and glider measurements using a density threshold of 28.3 kg m -3 (defining the base of the pycnocline).

Under a continental prevailing wind and a moderate river discharge (Situation 1 in February 4, Fig. 6a), the vertical distribution of surface SPM concentrations closely corresponded to the thermo-haline structure of the surface plume. Satellite observations showed that the plume presented a seaward gradient of SPM concentrations from 28 to 10 mg L -1 from the river mouth to 22 km offshore. Glider surface observations showed a similar decrease from 13 to 2 mg L -1 with the largest SPM concentrations (> 5 mg L -1 ) restricted by the 27 kg m -3 isopycnal. A thin BNL was observed over the shelf and presented SPM concentration of 1-2 mg L -1 .

During a marine prevailing wind (Situation 2 on February 7, Fig. 6b), the plume remained over the shelf and thickened to 30 m depth. The entire water column showed a SPM concentration of 2 mg L -1 . The most turbid part of the plume restricted by the 27 kg m -3 isopycnal was pressed against the coast. Satellite observations showed a similar trend and were completed with glider measurements that showed a more turbid BNL (2.5-5 mg L -1 ) 10 km off the river mouth. During the storm event (Situation 3 on February 10, Fig. 6c), glider measurements showed the thickening of the plume down to 50 m depth. At the same time, the BNL presented higher SPM concentrations close to the coast. During the subsequent flood event (Situation 3 on February 11), the plume extended up to 45 km offshore and showed the highest SPM concentrations between 10 and 25 km off the river mouth. SPM properties variability -We used the spectral slopes of glider and satellite-derived backscattering measurements, γ Glider and γ MODIS , as indexes of the particle assemblage properties (Fig. 7). Both indexes ranged from 0.3 to 1.4 and highlighted variable particulate assemblages in nepheloid layers depending on meteorological and sea conditions. From our results, both indexes were clearly linked with the SPM concentration: the higher the SPM concentration, the higher the proportion of large flocs in the particulate assemblage (i.e., low γ). The lower γ values were generally observed in the most turbid area of the plume restricted by the 27 kg m -3 isopycnal.

Under a continental prevailing wind and a moderate river discharge (Situation 1 in February 4, Fig. 7a), both γ indexes showed that the plume close to the river mouth was characterized by increasing values with distance from the river mouth. Larger γ values were also observed in the BNL. During a marine prevailing wind (Situation 2 on February 7, Fig. 7b), when the plume was pressed against the coast, γ indexes decreased in both the inner parts of the plume and the BNL. During the storm and subsequent flood event (Situation 3 on February 11, Fig. 7c), γ generally decreased in the plume that significantly thickened and widened. Close to the seabed, the BNL showed lower γ values close to the coast. The time the satellite image was taken is shown by a black square on each glider section.

From top to bottom, the three situations discussed.

Hydrodynamical and hydrological derived parameters

Along-shelf (E-W) current component, the Brunt-Väisälä frequency (N), and Richardson number (Ri) derived from glider measurements are presented in Figure 8. We present these observations together with wind and wind stress (τ) measurements for the three situations (Fig. 8a).

During the two first situations (left and middle panel in Fig. 8), the results showed the high subsurface (1-5 m) stratification of the plume. The Richardson number highlighted a probable mixing offshore (i.e., Ri < 0.25) between the frontal zone of the plume and outer shelf waters.

During the third situation (Fig. 8c right panel), wind stress showed maximum value of 0.38 N m -2 under a marine prevailing wind (February 10, 2014 at 9:00 AM). At that time, the currents close to the coast were the strongest (50 cm s -1 ) and in opposite direction between the surface and the bottom layers. Moreover, N values were low throughout the water column and Ri values were lower than 0.25. During the subsequent flood and under a continental prevailing wind, high N between 10 and 40 km offshore highlighted the plume thickening (Fig. 8c).

Discussion

Several studies emphasized the necessity of combining satellite and glider platforms and performing high spatial and temporal resolution measurements to properly characterize SPM dynamics in the coastal area, especially during energetic events (Castelao et al., 2008b;[START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF][START_REF] Miles | Glider observations and modeling of sediment transport in Hurricane Sandy[END_REF][START_REF] Saldías | Optics of the offshore Columbia River plume from glider observations and satellite imagery[END_REF]. Based on simultaneous use of satellite and glider observations, we intend here to better understand the variability of the dynamics of nepheloid layers and of particle characteristics in the Rhône River ROFI.

Nepheloid layers dynamics

The extent of the surface plume was derived both from satellite and glider data. We examine here the consistency of the estimates based on these two data sets. The estimate of the plume extent from satellite data was inferred solely from turbidity using a threshold value of 10 mg L -1 defined by [START_REF] Ody | Potential of high spatial and temporal ocean color satellite data to study the dynamics of suspended particles in a non-tidal river plume[END_REF]. It appears that this threshold in fact delimited the core of the plume when compared to the extent of the plume as defined by the glider data. The glider data allowed us to define the buoyant plume both in terms of turbidity and density, the latter being a more physically relevant parameter. Consequently, the plume dimension inferred from the satellite derived turbidity threshold rather matched the core of the plume delimited by the 27 kg m -3 isopycnal, as observed by the glider. This implies that a turbidity limit lower than 10 mg L -1 (typically 1 mg L -1 based on glider measurements (see. Fig. 2)) should have been considered using satellite data to define the actual size of the plume delineated by the 28.3 kg m -3 isopycnal. Hence, care must be taken when using only a turbidity criterion from satellite images to define the plume size.

Previous studies based on models [START_REF] Marsaleix | A numerical study of the formation of the Rhône River plume[END_REF][START_REF] Estournel | The Rhone River plume in unsteady conditions: numerical and experimental results[END_REF] and observations from radars and drifters [START_REF] Broche | Experimental study of the Rhone plume. Part I: physics and dynamics[END_REF][START_REF] Naudin | Effect of mixing on microbial communities in the Rhone River plume[END_REF] highlighted the variability of the thermo-haline structure of the surface plume under different wind and river discharge conditions. Remote sensing of the surface SPM completed these observations and showed the role of the wind-driven circulation in SPM transport, especially during flood events [START_REF] Arnau | River mouth plume events and their dispersal in the Northwestern Mediterranean Sea[END_REF][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 non-tidal river plume[END_REF]. Our observations are in line with these previous works, showing evidence of the agreement between the turbid and freshwater plumes and their rapid temporal variability.

Under continental prevailing wind and moderate river discharge (Situation 1), the plume was deflected and rapidly diluted offshore over the continental shelf, and a BNL developed beneath it. This result suggest a persistent character of the vertical transfer of SPM in the ROFI towards the BNL during continental wind periods, probably because of the large particles sinking from the plume as shown by [START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF].

Under a marine prevailing wind (Situation 2), the plume was pressed against the coast, increasing the turbidity in both the surface plume and the BNL close to the river mouth (Fig. 6b). These observations suggest an increasing transfer of the Rhône River SPM to the BNL. Such a direct supply, which was poorly documented due to sampling difficulties [START_REF] Aloisi | Origine et rôle du néphéloïde profond dans le transfert des particules au milieu marin. Application au Golfe du Lion[END_REF]Lorthois et al., 2012), emphasizes the role of marine wind in the increase of the vertical transfer of particles from the plume toward the BNL. This transfer is believed to contribute to the sedimentation over the inner-part of the shelf where high sedimentation rates (20-40 cm y - 1 ) were estimated [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][START_REF] Drexler | Stratigraphic signatures due to flood deposition near the Rhône River: Gulf of Lions, northwest Mediterranean Sea[END_REF][START_REF] Marion | In situ record of sedimentary processes near the Rhône River mouth during winter events (Gulf of Lions, Mediterranean Sea)[END_REF].

Similar observations were made in other micro-tidal areas, e.g., the Po River ROFI [START_REF] Palinkas | The use of 7Be to identify event and seasonal sedimentation near the Po River delta, Adriatic Sea[END_REF][START_REF] Milligan | Flocculation and the loss of sediment from the Po River plume[END_REF].

During the storm event (Situation 3 on February 10), the occurrence of large waves and a strong westward coastal jet probably enhanced the particle resuspension over the prodelta down to 70 m deep, and at the same time, transported river and resuspended material along the coast. We estimated a bottom shear stress induced by waves and currents of 0.2 Pa, which can resuspend non-cohesive sediment up to 25 µm in grain size according to the Sedtrans05 model (Neumeier et al., 2005). It stresses the role of storm events in the sediment redistribution from the prodelta to the shelf, and these findings are in agreement with previous studies on the stratigraphy off the Rhône prodelta by [START_REF] Drexler | Stratigraphic signatures due to flood deposition near the Rhône River: Gulf of Lions, northwest Mediterranean Sea[END_REF] and [START_REF] Marion | In situ record of sedimentary processes near the Rhône River mouth during winter events (Gulf of Lions, Mediterranean Sea)[END_REF] and on the modelling of sediment transport by [START_REF] Ulses | Suspended sediment transport in the Gulf of Lions (NW Mediterranean): Impact of extreme storms and floods[END_REF] and [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]. Additionally, the low Brunt-Väisälä frequency and Richardson number (see Fig. 8 right panel) in this area indicated a breaking of the plume and mixing throughout the water column close to the coast. Several studies [START_REF] Kirinus | The contribution of waves in mixing processes of the Patos Lagoon plume[END_REF][START_REF] Rogowski | Observations of the Frontal Region of a Buoyant River Plume Using an Autonomous Underwater Vehicle[END_REF][START_REF] Rong | Current-wave interaction in the Mississippi-Atchafalaya river plume on the Texas-Louisiana shelf[END_REF] described such a mixing of a buoyant plume induced by current-wave interactions.

Spatio-temporal variability of SPM properties

We compared the spectral slopes of glider (γ Glider ) and satellite-derived (γ MODIS )

backscattering measurements as indexes of the SPM size. Both indexes ranged from 0.3 to 1.4, consistent with the estimates of [START_REF] Niewiadomska | Submesoscale physicalbiogeochemical coupling across the Ligurian current (northwestern Mediterranean) using a bio-optical glider[END_REF], [START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF] and [START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF] in the northwestern Mediterranean Sea, and showed similar decreasing trends with regards to SPM concentration (Fig. 5). These observations are consistent with previous studies of [START_REF] Boss | Shape of the particulate beam attenuation spectrum and its inversion to obtain the shape of the particulate size distribution[END_REF] and [START_REF] Slade | Spectral attenuation and backscattering as indicators of average particle size[END_REF], where the authors inferred a similar decreasing trend. It is noteworthy that γ Glider followed the same tendency, no matter the depth. γ MODIS was significantly larger than γ Glider , implying that satellite measurements systematically saw finer particle assemblages, especially for large SPM concentrations, than glider measurements. The high stratification of the plume could be the main factor involved in differences observed between satellite, which primarily sensed the topmost meter (see 4.3), and the glider, which sampled below 1 m depth. Higher SPM concentrations might be found in the very superficial plume water due to the strong stratification within the plume (see Fig. 9). As a matter of fact, [START_REF] Ouillon | An inversion method for reflectance in stratified turbid waters[END_REF] stressed such stratification in the Ebro River plume, which presents similar characteristics to the Rhône River plume studied in this work. The density gradient in the top meters of the plume could also limit the mixing with underlying waters and thus limit the action of flocculation. The transition between both surface and subsurface layers within the plume appeared to be decisive in the variability of the SPM concentration and size and requires further investigations using an appropriate strategy. For instance, the glider could surface after each profile instead of "licking" the surface before diving, and also leave its sensors on while transmitting data and receiving mission orders during surfacing. Additionally, whereas [START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF] proved in the same region under similar conditions that γ Glider is representative of the SPM size spectra, γ MODIS has not been calibrated against direct measurements of the SPM size spectra. As a matter of fact, [START_REF] Sahay | Backscattering of light by coastal waters using hyperspectral in-situ measurements: A case study off Veraval, Gujarat, India[END_REF] showed that the derivation of the particulate backscattering coefficients (bbp(λ)), derived from remote sensing through the QAA approach (see 3.5), could present an overestimate of 20% and 43%

for the 470 nm and 700 nm wavelengths, respectively. They based these estimates on in situ measurements and emphasized the variability existing between study areas and water types.

This implies the need for an empirical relationship to improve the estimate of bbp through the QAA approach in the Rhône ROFI using in-situ calibrations of the SPM size, especially close to the river mouth where the highest SPM concentrations (i.e., > 35 mg L -1 ) were measured and corresponded to saturated bbp and non-representative γ MODIS .

Because the spatial and temporal distributions of γ MODIS at the level of the surface plume showed the same trends with that of γ Glider , it can be considered that, apart from the offset, the results are complementary. For marine wind conditions, when the Rhône River plume was pressed against the coast (Fig. 7b), γ MODIS showed lower values that coincide with minimum γ Glider values in the whole water column. This likely results from a higher proportion of large flocs that formed due to increasing SPM concentrations. Such flocs, which have a larger settling velocity, can rapidly reach the bottom layer and contribute to sustaining the BNL.

This observation supports the hypothesis of [START_REF] Naudin | River discharge and wind influence upon particulate transfer at the land-ocean interaction: case study of the Rhone River plume[END_REF][START_REF] Thill | Evolution of particle size and concentration in the Rhône river mixing zone: influence of salt flocculation[END_REF] and [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], who identified the role of the flocculation in the rapid settling of particles close to the coast during flood events. [START_REF] Fox | Flocculation and sedimentation on the Po River Delta[END_REF] and [START_REF] Milligan | Flocculation and the loss of sediment from the Po River plume[END_REF] also stressed the role of the flocculation in the fate and settling of SPM from the surface plume close to the Po river mouth in the Adriatic. During the storm event (Fig. 7c), γ Glider showed that the particle assemblage near the river mouth was finer throughout the water column. This diminution of the particle size is believed to result both from the breakup of flocs by turbulent motions induced by the large waves and a strong along-shore current.

Moreover, wave-induced resuspension of sediment down to 70 m depth could enhance the flocculation close to the seabed and accentuate the SPM transport off the Rhône River ROFI in the BNL (see similar results in [START_REF] Curran | Settling velocity, effective density, and mass composition of suspended sediment in a coastal bottom boundary layer, Gulf of Lions, France[END_REF] and [START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF]).

Conclusion

The TUCPA experiment offered a unique set of observations on the high-frequency spatial and temporal variability of the Rhône River ROFI for different wind conditions and river discharges, including a flood event. The combination of coastal buoy monitoring, glider transects and satellite observations provided a detailed description of the hydrology, hydrodynamics and SPM characteristics in the ROFI.

The results emphasized the high variability of the Rhône River plume extent as well as the persistence of a bottom nepheloid layer. Three typical meteorological situations that are involved in the general transport and sedimentation of Rhône River particles have been detailed. Continental wind enhanced the offshore export of SPM within the surface river plume, the thickness and concentration of which vary with the river discharge. Conversely, marine wind, which leads to wavy and flood conditions, limited the offshore dispersion of the river plume. This situation, thus, increased the surface SPM concentration and favored the formation of large particles by flocculation and their settling toward the BNL. Moreover, the BNL was also fed by the resuspension of sediment due to wave action. During a storm event, the results emphasized the role of waves and currents on the breakage of the plume stratification and resuspension of sediments down to a depth of 70, which favored the export of SPM off the ROFI.

The combination of both satellite and glider SPM measurements was determinant to monitor both surface (< 1 m depth, satellite data) and subsurface (> 1 m depth, glider data) parts of the river plume. The results from both platforms showed a similar pattern (i.e., increasing abundance of large particlestypically flocswith increasing SPM concentrations), but different ranges, probably due to the strong stratification of the plume. The estimates of the index of the particle size distribution, γ, highlighted the larger proportion of flocs within the plume as well as the change of the particle assemblage, which became finer seaward in both the river plume and BNL due to the settling of flocculated particles. Finally, the use of this index permitted the identification of the role of rough weather on the breakup of flocculated particles during the storm event that favored the transport of finest sediments off the ROFI.

Figure 1 -

 1 Figure 1 -RGB geo-referenced OLI image (Landsat-8) of the Rhone River plume (23 February 2014). Glider cross-shelf deployment (30 January -12 February, 2014 -7 sections) is shown in yellow. The position of the coastal MESURHO buoy is shown by a gray dot, 2.6 km off the Rhône river mouth. The Beaucaire gauging station and Cap-Couronne wind station are shown by yellow and orange dots, respectively. Depth contour lines are shown in gray.

Figure 2 -

 2 Figure 2 -Example of vertical profiles (average: solid line; standard deviation: shaded region) of potential density anomaly (red), Brunt-Väisälä frequency (blue) and SPM concentration (green) for different glider sections performed in the Rhône River plume in February 2014. A threshold of 0.01 s -1 is used to define the lower limit of the pycnocline (gray line). Lower limit of the plume in terms of SPM concentration and density anomaly are shown in green and black, respectively.

Figure 3 -

 3 Figure 3 -Time-series from January 30 to February 15, 2014 of: (a) the hourly mean wind speed and direction, (b) the Rhône River daily-mean discharge, (c) the significant wave height and direction, (d) the hourly mean current speed and (e) direction, and (f) the backscatter index. The duration of the glider deployment is shown in yellow above (a). Black triangles show when the glider was close to the MESURHO buoy. The dash-line boxes show the 3 situations discussed. Wind and wave directions indicate the direction from which the wind and waves are coming.

Figure 4 -

 4 Figure 4 -Cross-shelf hydrological sections of temperature (left panel) and salinity (right panel) measured by the glider. Density contour lines are shown in black and white on the temperature and salinity sections, respectively. White arrows in the bottom-left corner of each

Figure 5 -

 5 Figure 5 -Relationship between glider and satellite-derived a) (SPM) mg L -1 and b) γ (dimensionless). Relationship between γ Glider/MODIS and (SPM) Glider/MODIS mg L -1 for the 0-35 mg L -1 range. Satellite (February 11, 2014) Satellite measurements are in gray. Glider measurements are in black.

Figure 6 -

 6 Figure 6 -Left panels: Surface SPM concentrations retrieved from satellite data. The distance from the coast (km) and along the glider track is shown in black. Date and wind speed and direction are specified. The threshold of 10 mg L -1 used for comparison with glider measurements is shown in black. Clouds are shown in white. Right panels: Corresponding glider sections of SPM concentrations. Density contour lines are shown in black. The MESURHO buoy is shown by a black triangle above the first section. The time the satellite image was taken is shown by a black square on each glider section. From top to bottom, the three situations discussed. Note the different color bars used for the satellite and glider observations.

Figure 7 -

 7 Figure 7 -Left panels: Surface estimated γ MODIS from satellite data. The distance from the coast (km) and along the glider track is shown in black. Date and wind speed and direction are

Figure 8

 8 Figure 8a) Time-series of wind speed and direction (m s -1 ). The wind stress is shown in red (N m -2 ). Cross-shelf section of b) the along-shelf current (m s -1 ), c) the derived Brunt-Väisälä frequency (s -1 ) and d) the Richardson Number (dimensionless). Density contour lines are shown in black on the upper panel and in white on middle and lower panels.

Figure 9 -

 9 Figure 9 -Mean surface profiles of a) density anomaly (red) and b) SPM concentration (blue) and γ (green) estimated from glider measurements in the Rhône River ROFI. The standard deviation of each measurement is shown by shaded areas. Blue and green circles on b) show satellite-based estimates of the SPM concentration and γ. From top to bottom, the three situations discussed.

Table 1 -

 1 Meteorological and hydrodynamical conditions at the MESURHO buoy station for the 3 situations discussed.

	Situation 1	Situation 2	Situation 3*

Acknowledgments

The corresponding author is funded through a PhD grant from the French ministry. This work was funded by the TUCPA project under the CNRS EC2CO DRILL Program, the MATUGLI project under the ANR ASTRID DGA Program (ANR-14-ASTR-0021-01), the AMORAD project under the ANR Program (ANR-11-RSNR-0002), and the MISTRALS-MERMEX Program. LOV work was funded by the European Community's Seventh Framework Program under grant agreement n°606797 (HIGHROC project) and by the French Space Agency through the TTC TOSCA project. The MESURHO station deployment is part of the MESURHO Program involving IFREMER, IRSN, CNRS, CETMEF, and Phares et Balises.

We thank the SOERE MOOSE for supporting and providing long-term observation data in the Gulf of Lions. We especially thank the French glider team, P. Testor and A. Boss for their help on the glider data processing. We also give a special thank you to S. Dykstra and Y.

Leredde for their helpful discussion.


Massimo Tisi
email: massimo.tisi@imt-atlantique.fr

Zheng Cheng
email: zheng.cheng@inria.fr

CoqTL: an Internal DSL for Model Transformation in Coq

Keywords: Model-Driven Engineering, Model Transformation, Interactive Theorem Proving, Coq

In model-driven engineering, model transformation (MT) verification is essential for reliably producing software artifacts. While recent advancements have enabled automatic Hoare-style verification for non-trivial MTs, there are certain verification tasks (e.g. induction) that are intrinsically difficult to automate. Existing tools that aim at simplifying the interactive verification of MTs typically translate the MT specification (e.g. in ATL) and properties to prove (e.g. in OCL) into an interactive theorem prover. However, since the MT specification and proof phases happen in separate languages, the proof developer needs a deep knowledge of the translation logic. Naturally any error in the MT translation could cause unsound verification, i.e. the MT executed in the original environment may have different semantics from the verified MT. We propose an alternative solution by designing and implementing an internal domain specific language, namely CoqTL, for the specification of declarative MTs directly in the Coq interactive theorem prover. Expressions in CoqTL are written in Gallina (the specification language of Coq), increasing the possibilities of reuse of native Coq libraries in the transformation definition and proof. In this paper we introduce CoqTL, we evaluate its practical applicability on a case study, and identify its limitations.

Introduction

Model-driven engineering (MDE), i.e. software engineering centered on software models and MTs, is widely recognized as an effective way to manage the complexity of software development. With the increasing complexity of MTs (e.g., in automotive industry [START_REF] Selim | Model transformations for migrating legacy models: An industrial case study[END_REF], medical data processing [START_REF] Wagelaar | Using ATL/EMFTVM for import/export of medical data[END_REF], aviation [START_REF] Berry | Synchronous design and verification of critical embedded systems using SCADE and Esterel[END_REF]), it is urgent to develop techniques and tools that prevent incorrect MTs from generating faulty models. The effects of such faulty models could be unpredictably propagated into subsequent MDE steps, e.g. code generation.

Deductive verification is a promising approach for quality assurance in MT: correctness is specified by MT developers using contracts (i.e. pre/postconditions), then the semantics of the MT language together with contracts and metamodels are encoded into a deductive theorem prover. Thanks especially to recent advancements in SMT solvers, automatic deductive verification is giving good results in several scenarios [START_REF] Büttner | Verification of ATL transformations using transformation models and model finders[END_REF][START_REF] Büttner | On verifying ATL transformations using 'off-theshelf' SMT solvers[END_REF][START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF][START_REF] Oakes | Fully verifying transformation contracts for declarative ATL[END_REF]. However, because of the general undecidability, interactive deductive verification is inevitable for complex tasks (for instance, automatic deductive theorem provers usually lack support for induction or finding witnesses for existential quantifiers).

Coq is an interactive theorem prover. The user can use Coq to write mathematical definitions, executable algorithms and theorems together with an environment for semi-interactive development of proofs (in the sense that routine proofs can be automatically performed while difficult proofs require human guidance). It has been used to prove non-trivial mathematical theorems, or as an environment for developing formally certified software and hardware (e.g. [START_REF] Leroy | Formal certification of a compiler back-end or: Programming a compiler with a proof assistant[END_REF][START_REF] Gu | CertiKOS: An extensible architecture for building certified concurrent os kernels[END_REF]). While not strictly needed for understanding this paper, we refer the reader to [START_REF] Pierce | Software Foundations[END_REF] for an introduction to the Coq system.

Previous work aiming at simplifying the interactive verification of MTs, has already proposed translations from MT specifications (e.g. in MT languages like ATL) and properties to prove (e.g. in OCL) into Coq. However, the practical applicability of this translational approach is hampered by the fact that the two phases of MT specification and correctness proof require developments in languages (e.g. ATL+OCL and Coq, respectively) at two different levels of abstraction. The proof developer needs a deep knowledge of the translation logic to be able to write meaningful proofs. Any change in the MT code propagates through the translator, and it is difficult to predict the proof steps that will be invalidated. Naturally any error in the MT translation could cause unsound verification, i.e. the MT executed in the original environment may have different semantics from the verified MT. Certifying that the semantics of the MT language has being correctly axiomatized in the back-end theorem prover is a hard task, and very few attempts exist [START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF][START_REF] Ab | A survey of approaches for verifying model transformations[END_REF].

Coq includes Gallina, a functional programming language with pattern matching and rich type system, well suited as a platform for embedding domain-specific programming languages (DSLs) (e.g. [START_REF] Chlipala | The Bedrock structured programming system: Combining generative metaprogramming and hoare logic in an extensible program verifier[END_REF]). In this work, we draw on this aspect of Coq and propose a DSL, namely CoqTL, to turn Coq into a tool for developing certified MTs. We argue that using an internal DSL for the MT specification phase simplifies the iterative process of MT development and proof in MDE. Moreover, expressions in CoqTL are directly written in Gallina, increasing the possibilities of reuse of sophisticated native Coq libraries during the transformation definition and proof.

Our main contributions are:

-We design and implement CoqTL, to our knowledge the first DSL for rulebased MT in Coq (Section 3.2). The language is both functional and declarative in style, its syntax and semantics is inspired from ATL [START_REF] Jouault | ATL: A model transformation tool[END_REF] (hence it should be familiar also to users of other rule-based MT languages, like ETL [START_REF] Kolovos | The Epsilon transformation language[END_REF], or RubyTL [START_REF] Cuadrado | RubyTL: A practical, extensible transformation language[END_REF]). Thus, CoqTL aims to lighten the cognitive load of MT developers trying to build certified MTs in Coq. -We design and implement a transformation engine in Coq that interprets programs written in CoqTL to transform models (Section 3.3). The engine includes an on-the-fly parser that transforms the domain-specific syntax into a Coq data structure to interpret. The parser is transparently invoked (by an extensive use of the Coq Notation mechanism) so that any Coq development environment is able to support the domain-specific CoqTL syntax without requiring ad-hoc modifications. -We show the practical applicability of CoqTL, by using it to specify a sample transformation, prove non-trivial contracts over it and automatically extract a certified implementation.

We make CoqTL publicly available as open source 3 . The repository contains also the example and proofs described in this paper. Paper organization. We motivate our work by a sample transformation in Section 2. Section 3 illustrates the design of CoqTL in detail. In Section 4 we prove theorems on a CoqTL specification. Section 5 compares our work with related research, and Section 6 draws conclusions and lines for future work.

Class to Relational in CoqTL

We consider a very simplified version of the transformation from class diagrams to relational schemas (arguably, the Hello World transformation in the MT community). The example is intentionally very small, so that it can be completely illustrated within this paper. However we believe it to be easily generalizable by the reader to more complex scenarios. The structure of the involved metamodels is shown in Fig. 1. Each class contains a list of named and typed attributes with unique identities. In this simplified model we do not consider attribute multiplicity (i.e., all attributes are single-valued). Primitive data types are not explicitly modeled, thus we consider every attribute without an associated type to have primitive data type. A derived feature identifies which attributes are derived from other values. The simplified structural metamodel of relational schemas is shown on the right part of Figure 1. Tables contain Columns, Columns can refer to other Tables in case of foreign keys.

In Listing 1.1 we use the CoqTL language to specify how to transform class diagrams to relational schemas. A transformation is a Coq Definition. First, we declare that a transformation named Class2Relational is to transform a model conforming to the Class metamodel to a model conforming to the Relational metamodel, and we name the input model as m (lines 2-3).

Then, the transformation is defined via two rules in a mapping style: one maps Classes to Tables, another one maps non-derived Attributes to Columns. Each rule in CoqTL has a from section that specifies the input pattern to be matched in the source model. A boolean expression in Gallina can be added as guard, and a rule is applicable only if the guard evaluates to true for a certain assignment of the input pattern elements. Each rule has a to section which specifies elements and links to be created in the target model (output pattern) when a rule is fired. The to section is formed by a list of labeled outputs, each one including an element and a list of links to create. The element section includes standard Gallina code to instantiate the new element specifying the value of its attributes (line 13). The links section contains standard Gallina code to instantiate links related to the previous element (lines [START_REF] Picard | Coinductive graph representation: the problem of embedded lists[END_REF][START_REF] Pierce | Software Foundations[END_REF][START_REF] Poernomo | Correct-by-construction model transformations from partially ordered specifications in Coq[END_REF].

For instance in the Class2Table rule, once a class c is matched (lines 6 to 8), we specify that a table should be constructed by the constructor BuildTable, with the same id and name of c (line 13). While the body of the element section (line 13) can contain any Gallina code, it is type-checked against the element signature (line 12), i.e. in this case it must return a Table .  In order to link the generated table t to the columns it contains, we get the attributes of the matched class (line 17), resolve them to their corresponding Columns, generated by any other rule (line 18), and construct new set of links connecting the table and these columns (line 19). While this is standard Gallina code, we use for this example an imperative style with a monadic notation (←, similar to the do-notation in Haskell) that makes the code more clear in this case 4 . The resolveAll function will only return the correctly resolved attributes. In particular derived Attributes do not generate Columns (i.e. they are not matched by Attribute2Column), so they will be automatically filtered out by resolveAll. The result of this Gallina code (i.e. the constructed links) are type-checked against the link signature (i.e. in this case they must have type TableColumns, as specified at line 16).

In the Attribute2Column rule we can notice the presence of a guard. When the Attribute is not derived, a Column is constructed with the same name and identifier of the Attribute. If the original attribute is typed by another Class we build a reference link to declare that the generated Column is a foreign key of a Table in the schema. This Table is found by resolving (resolve function) the Class type of the attribute.

CoqTL naturally enables deductive verification of model transformations.

Users can write Coq theorems that apply pre/postconditions (correctness conditions) to the model transformation. For example, Listing 1.2 defines a theorem stating that if all elements contained by the input model have not-empty names, by executing the Class2Relational MT, all generated elements in the output model will also have not-empty names. Interactively proving this simple theo-

( ∀ (o : Table), In o ( allModelElements rm)→ length (getTableName i)>0).
Listing 1.2. Name definedness theorem for the Class2Relational transformation rem in Coq takes 56 lines of routine proof code (this short proof can be even automated by using modern automatic theorem provers [START_REF] Büttner | On verifying ATL transformations using 'off-theshelf' SMT solvers[END_REF][START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF]).

To illustrate more complex theorems we want to prove that our transformation preserves unreachability. (Un)reachability is an important property for several models, e.g. one may typically need to demonstrate that error states in generated state machines are not reachable. In our simple Class2Relational example, one can inductively define reachability for classes (similarly for tables), i.e. a class is reachable from itself, and two classes are reachable if they are transitively linked by attributes. We can define an unreachability preservation theorem as follows: if a certain class is not reachable from a given class, their corresponding tables will not be reachable from each other. Interactively proving this theorem in CoqTL needs more than a thousand lines of proof code. The major difficulty comes from choosing the right induction strategy, and to our knowledge, the automatic proof of similar theorems is not addressed by existing work. The full proof in Coq is available on the paper website.

The Design of CoqTL

CoqTL is an internal DSL for model transformation in Coq. In this section we will describe the three main parts of its design:

-(Section 3.1) Metamodels and models are encoded as graph structures that can be automatically translated from/to EMF. -(Section 3.2) Transformation specifications are encoded as a data structure wrapped up in a user-friendly domain specific syntax. -(Section 3.3) A transformation engine interprets transformation specifications against input models.

Metamodels and Models

Our encoding of metamodels in Coq is similar to analogous encodings in related work, based on inductive data types. As an example, Listing 1.3 shows the basic definitions for encoding the Relational metamodel of Fig. 1. Since this interface is the main means to access source and target models, we aim at providing the simplest native representation. Each metaclass is represented by an inductive data type, with a single constructor whose arguments are the attributes of the metaclass. References between metaclasses are represented as separate inductive types, with a constructor requiring the source and target elements as arguments. Optional or multivalued attributes and references are respectively represented using the option and list Coq types in the appropriate constructor argument (e.g. at line 15).

Constructing any model requires providing a list of model elements and one of links, as specified by the Model type in the CoqTL library (shown in lines 23-26 in the listing). These lists are typed by generic ModelElement and ModelLink types, that are meant to be the sum types for elements and links of the specific metamodel. For defining the type of Relational model, we first define the two sum types RelationalModelElement, sum of Table and Column, and RelationalMod-elLink, sum of TableColumns and ColumnReference (for simplicity here we omit the definition of sum types, that relies on dependent types). The RelationalModel type is obtained by parametrizing Model with these sum types.

We create accessors for every attribute and reference of each metaclass. Notice that while attribute accessors need only to inspect the element passed as argument to retrieve the attribute value (e.g., getTableId and getTableName at lines 37-44), reference accessors need to pass through the list of links to find the ones connected to the element in parameter. Thus, reference accessors need to have the whole model as extra parameter (e.g., getTableColumns in the listing).

Listing 1.3 shows a small portion of the encoding of the Relational metamodel in Fig. 1. The full encoding takes over 300 lines of Gallina code, and includes a reflective API. Briefly, metamodel classes are reified in a Relational-MetamodelClass type (with values corresponding to Table and Column), that is used as argument to reflective functions. The reflective API can be used for obtaining the metaclass of an element, checking that an element is an instance of a metaclass, and casting a generic element to/from a specific metaclass. Similar functions are available for links.

While our representation allows us to encode any model instance, in our current prototype we do not directly implement several features that are found in modeling frameworks like EMF. Bidirectional references currently have no special treatment: both sides are encoded as separate references, that need to be separately assigned in the transformation code. No direct support is provided for metaclass inheritance: the instance of a superclass can be provided as parameter of a subclass constructor, but the two instances (of superclass and subclass) need to be managed separately. Constraints for reference multiplicity or strong containment can only be encoded via extra pre/postconditions. Finally, differently from EMF, identifiers are considered as normal attributes and elements are considered equal when all their attributes are.

Automatic translators to/from EMF are still under development, and only partial implementations are provided on the CoqTL website. The way we implement the concrete syntax of CoqTL relies on the Notation facility of Coq. A notation is a symbolic abbreviation to denote some expressions, and is one of the main commands that modifies the way Coq parses and prints the representation of expressions.

For example, the first notation shown in Listing 1.4 implements the production rules link-def and link-decl in Grammar 1.1. After the declaration of this notation, when the expression on the left-hand-side is matched, it is expanded in memory to the right-hand-side. A notation allows also the specification of associativity and precedence levels, to solve parsing ambiguities. Notations can be seen as a very limited compiler, that compiles in one pass without memory. For this reason they strongly limit the classes of DSLs that can be implemented. In the implementation of CoqTL every notation is simply translated into an appropriate constructor, encapsulating the values matched by the notation (line 3). ( right associativity, at level 60).

Listing 1.4. A few notations for CoqTL

Whenever the notation is matching the declaration of some variable that needs to be visible to the rest of the code, we introduce a lambda expression as an argument of the constructor. This is shown in the second notation in Listing 1.4, that implements the output-elem production rule in the grammar. The created element elname needs to be visible in the following links section, so we store the content of this section in an anonymous function with elname as input (line 9).

The constructors used in our notations, like BuildOutputPatternLinkDefinition in Listing 1.4 build a representation of the abstract syntax of the CoqTL program. Hence CoqTL is a deeply embedded DSL for the rule structure part. CoqTL has however shallow embedding of expressions, to allow the direct use of the Gallina language for guards and output patterns (gallina-expr in the grammar).

Gallina has several characteristics that make it suitable as an expression language for CoqTL. It is:

-Expressive. Gallina is based on a formal language called the Calculus of Inductive Constructions, combining a higher-order logic and a richly-typed functional programming language. -Easy to learn. In our experience, the learning curve of the language is low if the user had some exposure to functional languages. -Accompanied by sophisticated libraries. Reusing functions in those libraries during the MT specification is also important for the proof phase, that can exploit the theorems and lemmas provided by the library for those functions.

Finally, CoqTL provides auxiliary functions meant to be used in Gallina expressions for guards and output patterns. The most important is the function resolve (and its corresponding multivalued version, resolveAll ) for element resolution. As illustrated at lines 18 and 36 in Listing 1.1, its signature requires the following arguments: 1) the result of the matching phase of the current transformation (match Class2Relational m), 2) the label associated to the required output element, useful for rules with multiple output elements ("col"), 3) the type of the expected result, useful for type checking (Column), 4) the source pattern to resolve (or the list of source patterns in case of resolveAll ). Notice that the matching phase is provided as a new application of the transformation in a specific match mode. While this choice affects the global efficiency of the transformation, it simplifies the development of proofs, because it avoids having a concept of transformation traces as side effects of the transformation execution.

While the expressiveness of CoqTL has important limitations, we are currently able to manually translate to CoqTL a significant subset of ATL transformations: 1-to-1 model transformations, in standard (non-refining) mode, written in the declarative subset of ATL, without lazy rules.

Transformation Engine

Algorithm 1 illustrates in pseudocode how the transformation specifications are interpreted by our transformation engine. This algorithm has been influenced by the execution algorithm of ATL [START_REF] Jouault | ATL: A model transformation tool[END_REF] (notably in the distinction between a match/instantiate and an apply function), but is very different, having the objective to simplify the proof development, at the cost of sacrificing execution efficiency.

Our transformation engine is implemented in an execute function (called for instance in Listing 1.2) that takes as input a transformation specification R and an input model I (which contains elements I e and links I l ). The output is elements O e and links O l , which form an output model.

First, the transformation engine records the maximum size (m) of input patterns among all the rules in the transformation specification. This value is used to calculate all the potential pattern instances P that the input model can produce to be matched against the transformation specification, i.e. all the subsets of I e whose size is less or equal to m are enumerated.

Next, the engine iterates on each potential pattern instance p, and seeks for a rule r in R that matches it (i.e. if model elements in the pattern instance have the types defined in the input pattern of the rule) and satisfies the guard of that rule. If a rule r is found for the pattern instance p, then the instantiation phase of r will be invoked to construct the corresponding output elements of p and add them to the output model. Finally the apply phase is invoked, i.e. to construct the corresponding output links and add them to the output model. if r = None then 6:

Oe ← Oe ∪ instantiate(r, p) 7:

O l ← O l ∪ apply(r, p) 8:

end if 9: end for Notice that Gallina expressions for output links are only evaluated during the apply phase. The developer may include in these expressions calls to the resolve or resolveAll functions, whose evaluation requires the execution of the instantiate phase. As mentioned in the previous section, in our solution the user passes to resolve the result of the transformation execution in match mode (i.e. match function at (lines 18 and 36 in Listing 1.1). The algorithm implemented in match is identical to Algorithm 1, without line 7 (that would make the whole computation recur indefinitely). Multiple executions of the transformation for element resolution slow down the execution, but simplify the proofs, since no explicit traces are necessary as applications of instantiate and apply with identical inputs can be trivially checked for equality. Possible optimizations are however the subject of future work.

Finally the application of the transformation by the execute function can be automatically extracted by Coq into a separate executable program in several languages (e.g., OCaml, Haskell).

Proving theorems with CoqTL

In this section we show that CoqTL can enable practical verification for MTs. We formulate 4 theorem proofs over the model transformation presented in Section 2. Some measures are shown in Table 1, to give the reader an idea of the complexity of the proofs: lines of code (LoC) and number of user-developed lemmas.

As a first theorem we prove that Class2Relational preserves id positivity, i.e. if all identifiers in the source model are positive, then they also are in the target model. In the first and second row we show two proofs for this theorem. In the second proof we obtain a reduction of about 60% LoC, thanks to the use of a generic lemma for transformation surjectivity, provided in the CoqTL library. This shows that CoqTL enables the design and proof of generic theorems that make interactive verification more efficient and concise.

Transformation surjectivity states that for all elements contained in the output model there has to exist a rule and a matching input pattern that created them. Our design choices in CoqTL enable this kind of theorems: during the proof we can refer to syntactic elements of the transformation (e.g. rules, input/output patterns) by their type in the abstract syntax (e.g., OutputPatternLinkDefinition in Listing 1.4), and quantify over them. Moreover we use the reflective model API mentioned in Section 3.1 to reason on metamodel-agnostic properties. The surjectivity lemma is also used in the third and fourth proof. In the third row we prove the name definedness property shown in Listing 1.2, separately for all element types in source and target models. Finally by the fourth row, it is clear that the unreachability preservation theorem (Section 2) is difficult to prove, and shows the need of further work in proof engineering for MTs.

One road we want to follow is providing a complete library of generic lemmas for CoqTL such as transformation surjectivity, to shorten proofs on CoqTL. Some recurring proof patterns could be factorized into domain-specific automatic proof tactics, aware of the CoqTL representation and properties. Another line could be investigating a set of domain specific guidelines to construct proofs for MT verification. For example, to prove that if two Tables are reachable, the Classes that generated them are reachable too; we induct on the definition of reachablity. However other induction strategies, e.g. on the structure of the model, may be more efficient.

Related Work

There are many automatic theorem proving approaches for MTs (e.g. [START_REF] Büttner | Verification of ATL transformations using transformation models and model finders[END_REF][START_REF] Büttner | On verifying ATL transformations using 'off-theshelf' SMT solvers[END_REF][START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF][START_REF] Oakes | Fully verifying transformation contracts for declarative ATL[END_REF]). However, interactive theorem proving is inevitable for more serious verification tasks. In this section, we focus on recent advancements of MT verification based on interactive theorem proving. To our knowledge, none of the existing works designs and implements DSLs for MT within interactive theorem provers.

Yang et al. interactively verify that a particular model transformation, i.e. from AADL to TASM language, is semantic preserving [START_REF] Yang | From AADL to timed abstract state machines: A verified model transformation[END_REF]. The approach is based on providing a translational semantics of both languages as timed transition systems in Coq and then reasoning on their equivalence. CoqTL could be used to simplify this kind of work.

Most previous works focus on giving a translational semantics of a MT language towards the target theorem prover. Generally they do not investigate a way to formally ensure that the semantics of the MT language has been axiomatized correctly in the back-end theorem prover. Calegari et al. encode ATL MTs and OCL contracts into Coq to interactively verify that the MT is able to produce target models that satisfy the given contracts [START_REF] Calegari | A type-theoretic framework for certified model transformations[END_REF]. In [START_REF] Stenzel | Formal verification of QVT transformations for code generation[END_REF], a Hoare-style calculus is developed by Stenzel et al. in the KIV prover to analyze transformations expressed in (a subset of) QVT Operational. UML-RSDS is a tool-set for developing correct MTs by construction [START_REF] Lano | A framework for model transformation verification[END_REF]. It chooses well-accepted concepts in MDE to make their approach more accessible by developers to specify MTs. Then, the MTs are verified against contracts by translating both into interactive theorem provers.

Kezadri et al. defines the Coq4MDE framework to formally embed some key aspects of MDE in Coq [START_REF] Hamiaz | A formal framework to prove the correctness of model driven engineering composition operators[END_REF]. We have a similar abstraction of metamodels as graphs. While our understanding is that Coq4MDE is capable of embedding MT languages and enabling MT verification, no specific work has been proposed. We expect an evaluation in the future to compare the complexity of MT verification between the two works.

Poernomo and Terrell follow the classical approach in type theory to formally specify MTs as ∀∃ types in interactive theorem provers [START_REF] Poernomo | Correct-by-construction model transformations from partially ordered specifications in Coq[END_REF]. Their approach does not target any specific MT languages. In addition, although their work does not propose a generic MT engine as we presented here, a corresponding executable MT program can be extracted once the MT is proved. The approach is further extended by Fernández and Terrell on using co-inductive types to encode bidirectional or circular references [START_REF] Fernández | Assembling the proofs of ordered model transformations[END_REF]. We also plan to investigate how co-inductive types can cooperate with our encoding and proofs (e.g. guardedness issues of corecursive functions might arise because the syntactic criterion applied by the Coq system is too rigid [START_REF] Picard | Coinductive graph representation: the problem of embedded lists[END_REF]).

Conclusion

In conclusion, we present CoqTL, to our knowledge the first DSL in Coq for MTs and their verification. CoqTL is both functional and declarative in style, providing a familiar environment for transformation developers in Coq. Its underlining transformation engine, implemented in Coq, allows CoqTL programs to be interpreted against input models to compute output models. We show the practical applicability of CoqTL, by proving non-trivial contracts over a sample transformation.

Our future work would focus on the issues we identified in different points of our discussion. We want to develop a theorem library on top of CoqTL to facilitate MT verification, including of transformation-agnostic lemmas such as transformation surjectivity and domain-specific proof tactics to automatize recurring proof steps. We aim to investigate whether there are domain-specific guidelines to construct proofs for MT verification. We want to improve interoperability between CoqTL and common MDE tools such as EMF, for industry readiness.

Fig. 1 .Listing 1 . 1 .

 111 Fig. 1. A simplified structural metamodel for class diagrams (left), and relational schemas (right)

  Theorem Table_name_definedness : ∀ ( cm : ClassModel) ( rm : RelationalModel), (* transformation *) rm = execute Class2Relational cm → (* precondition *) ( ∀ (i : Class), In i ( allModelElements cm)→ length (getClassName i)>0) → (* postcondition *)

(

  *** Metamodel classes and references ***) Inductive Table : Set := BuildTable : (* id *) nat → (* name *) string → Table.

  Inductive Column : Set := BuildColumn : (* id *) nat → (* name *) string → Column. Inductive TableColumns : Set := BuildTableColumns: Table → list Columns → TableColumns. Inductive ColumnReference : Set := BuildColumnReference: Column → Table → ColumnReference. (*** Model (from CoqTL library) ***) Inductive Model ( ModelElement: Type) ( ModelLink: Type): Type := BuildModel: list ModelElement → list ModelLink → Model ModelElement ModelLink. (*** Relational Model ***) Inductive RelationalModelElement : Set := ... (* sum type for elements *) Inductive RelationalModelLink : Set := ... (* sum type for links *) Definition RelationalModel := Model RelationalModelElement RelationalModelLink. (***

1 . 1 .

 11 transformation ::= header ':=' '[' rule-list ']' header ::= 'transformation' 'from' id 'to' id 'with' id 'as' id rule-list ::= rule ';' rule-list | rule rule ::= 'rule' id 'from' input-pattern 'to' output-pattern input-pattern ::= elem-decl-list 'when' gallina-expr elem-decl-list ::= elem-decl ',' elem-decl-list | elem-decl elem-decl ::= 'element' id 'class' id 'from' id output-pattern ::= '[' output-list ']' output-list ::= output-elem ';' output-list | output-elem output-elem ::= 'output' string 'element' elem-def 'links' '[' link-def-list ']' elem-def ::= elem-decl ':=' gallina-expr link-def-list ::= link-def ';' link-def-list | link-def link-def ::= link-decl ':=' gallina-expr link-decl ::= 'reference' id 'from' id Grammar The concrete syntax of the CoqTL language

(

  * Output Link Definition *) Notation "'reference' reftype 'from' tinstance ':=' refends" := ( BuildOutputPatternLinkDefinition tinstance reftype refends) ( right associativity, at level 60). (* Output Pattern Element *) Notation "'output' elid 'element' elname 'class' eltype 'from' tinstance := eldef 'links' refdef" := ( BuildOutputPatternElement eltype elid eldef ( fun elname ⇒ refdef))

Algorithm 1

 1 Algorithm of the execute function 1: m ← maxArity(R) 2: P ← allPatterns(Ie, m) 3: for each p ∈ P do 4: r ← findRule(R, p) 5:

Table 1 .

 1 Theorem proofs on Class2Relational

	Theorem	LoC No. Lemmas
	positive ids	180	4
	positive ids surj	75	1
	name definedness	89	2
	unreachability preservation 1161	17

CoqTL (online). https://github.com/atlanmod/CoqTL

the intuitive semantics of ← is: if the right-hand-side of the arrow is not None, then assign it to the variable in the left-hand side and evaluate the next line, otherwise return None

Acknowledgements. We thank Rémi Douence for his valuable help during the development of CoqTL.