Roger Marí
email: roger.mari@ens-paris-saclay.fr

Thibaud Ehret

Jérémy Anger

Carlo De Franchis

Gabriele Facciolo

L1B + : A PERFECT SENSOR LOCALIZATION MODEL FOR SIMPLE SATELLITE STEREO RECONSTRUCTION FROM PUSH-FRAME IMAGE STRIPS

Keywords: Rational Polynomial Coefficients, RPC model, push-frame, SkySat, satellite stereo

We propose a novel method to generate a single image product from a multi-image strip acquired by a push-frame satellite imaging system. The images of the push-frame strips are combined into a large scale mosaic simulating a perfect sensor geometry. The local camera models of the input images are leveraged to produce a new localization model that covers the output mosaic entirely. Among other applications, this simplifies the task of stereo reconstruction enormously: instead of treating multiple stereo pairs of small images, it is possible to reconstruct the entire area covered by the push-frame acquisition using a single pair of mosaics incorporating all the images. We test our method using strips of SkySat L1B scenes and denote the output images as L1B + . To evaluate the quality of the L1B + images and their localization models, the stereo reconstructions obtained with L1B + are compared with those obtained with L1B and with a lidar reference model.

INTRODUCTION

Push-frame satellite imaging systems acquire a continuous strip of small and partially overlapping images as the satellite moves [START_REF] Aati | Optimization of optical image geometric modeling, application to topography extraction and topographic change measurements using PlanetScope and SkySat imagery[END_REF]Avouac, 2020, Planet, 2021). Small satellite (Small-Sat) constellations, such as SkySat and PlanetScope from Planet or Aleph-1 from Satellogic, use a push-frame imaging mode to cover large areas of interest, beyond the limited footprint of their telescope. The resulting image collections are a product of great interest: small satellites can orbit at low altitudes and offer very high resolution, in direct competition with wellestablished large satellite providers. While the latter are capable of covering large areas with a single shot, the derived products are also more expensive. SmallSat constellations can also afford many more units, resulting in shorter revisit times [START_REF] Sandau | Small satellites for global coverage: Potential and limits[END_REF]. However, despite their economic and technical appeal, the highly fragmented nature of push-frame acquisitions may discourage their use for some applications targeting areas of interest of several km 2 [START_REF] Xue | Small satellite remote sensing and applications-history, current and future[END_REF].

We propose a novel method to exploit push-frame image collections from SmallSat acquisitions, which circumvents the difficulties of working with such fragmented data. The method combines the multiple images of a strip acquired by a push-frame system into a single mosaic image, as if it had been acquired by an instrument with perfect sensor geometry (Figure 1). The local camera models of the input images are leveraged to produce a new localization model that covers the output mosaic entirely. The resulting image and localization model can be treated as new product, which we denote L1B + .

Our contributions are:

• A method capable of assembling the frames and camera models of strips of partially overlapping satellite images.

• An evaluation of the method based on an application of major interest: stereo reconstruction. We validate the method using strips of SkySat L1B scenes. The 3D models obtained with the L1B + products are compared with those obtained with L1B and with a lidar reference model.

Figure 1. We combine the multiple images of each strip acquired by the push-frame system to produce large scale mosaics. The resulting perfect sensor images and their localization models can be used to greatly simplify the task of large scale stereo reconstruction from push-frame satellite imagery.

RELATED WORK

SmallSat push-frame imagery is used in a wide range of remote sensing applications, including topography extraction [START_REF] Aati | Optimization of optical image geometric modeling, application to topography extraction and topographic change measurements using PlanetScope and SkySat imagery[END_REF][START_REF] Bhushan | Automated digital elevation model (DEM) generation from very-high-resolution Planet SkySat triplet stereo and video imagery[END_REF], d'Angelo and Reinartz, 2021), super-resolution products [START_REF] Nguyen | Self-supervised multi-image super-resolution for pushframe satellite images[END_REF][START_REF] Anger | Fast and accurate multi-Frame super-resolution of satellite images[END_REF] and various tasks demanding short revisit times, such as monitoring of natural phenomena [START_REF] Cannistra | Highresolution CubeSat imagery and machine learning for detailed snow-covered area[END_REF] or commercial assets (Marí et al., 2021b, d'Autume et al., 2020). This paper focuses on the task of 3D reconstruction, but the proposed method can be beneficial for any of these applications.

Topography extraction from high resolution satellite images is typically performed using stereo pipelines [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF], de Franchis et al., 2014a), capable of producing highly accurate photogrammetric digital surface models (DSMs). These pipelines take as input one pair of images and their geometric camera models, represented with the rational polynomial coefficients (RPC) model. The RPC model represents the image acquisition system by means of a pair of rational polynomial functions that approximate the mapping from 3D space points to 2D image pixels: P : R 3 → R 2 , the projection function; and its inverse L : R 2 × R → R 3 , the localization function. The RPC camera models allow to represent unconventional geometric models such as that of push-broom scanners. Push-broom cameras are not projective as the optical center changes from line to line. However, for satellite images they can be approximated locally by affine models (de Franchis et al., 2014b). This is exploited by stereo pipelines (de Franchis et al., 2014a[START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF], which cut the images into small tiles that can be stereo-rectified and processed independently. This insight motivates our merged camera model for push-frame images and its usefulness in the stereo context.

Image collections acquired by push-frame systems cannot be directly plugged into the state-of-the-art 3D reconstruction pipelines, as the task becomes a multi-view problem in which each image will contribute only to some part of the final reconstruction. The majority of solutions that can be found in the literature address the reconstruction of an area of interest observed by multiple overlapping satellite images as a two-stage multi-view stereo (MVS) problem [START_REF] Facciolo | Automatic 3D reconstruction from multi-date satellite images[END_REF][START_REF] Gong | Point cloud and digital surface model generation from high resolution multiple view stereo satellite imagery[END_REF][START_REF] Bhushan | Automated digital elevation model (DEM) generation from very-high-resolution Planet SkySat triplet stereo and video imagery[END_REF]. In the first stage, a series of local models are computed, resulting from independent executions of a stereo reconstruction pipeline using different input pairs. In the second stage, the local models are fused to obtain the complete reconstruction of the area. This two-stage approach requires significant pre-and post-processing work.

The pre-processing tasks of satellite MVS require the selection of an optimal set of stereo pairs, to minimize the number of local models to be reconstructed and all the derived workload.

In addition to the amount of geographic footprint overlap, the incidence angles of the cameras, the angle between both views and the acquisition dates have proven to be important to select suitable stereo pairs [START_REF] Facciolo | Automatic 3D reconstruction from multi-date satellite images[END_REF][START_REF] Gong | Point cloud and digital surface model generation from high resolution multiple view stereo satellite imagery[END_REF][START_REF] Marí | To bundle adjust or not: A comparison of relative geolocation correction strategies for satellite multi-view stereo[END_REF][START_REF] Gómez | An experimental comparison of multi-view stereo approaches on satellite images[END_REF]. Another classic pre-processing step is the correction of camera positions and orientations, in order to make the multiple views geometrically consistent. Bundle adjustment algorithms have so far stood out as the best practice to correct the inaccuracies of the RPC models [START_REF] Triggs | Bundle adjustment-a modern synthesis[END_REF][START_REF] Grodecki | Block adjustment of highresolution satellite images described by rational polynomials[END_REF], Marí et al., 2021a). The corrected RPC models produced by the bundle adjustment methods minimize the reprojection error of a set of reference points observed across the different images. The differential advantage of adding a camera correction step prior to stereo reconstruction is that local models computed from independent pairs are, in principle, natively registered in the object space.

The post-processing tasks of satellite MVS revolve around the fusion step. Assuming that all local models are aligned, the usual procedure is to discretize the area of interest and apply a mean or median filter on each 2D cell to establish the final altitude. However, the result of the fusion is usually largely improvable and several works have introduced heuristics [START_REF] Facciolo | Automatic 3D reconstruction from multi-date satellite images[END_REF][START_REF] Qin | Automated 3D recovery from very high resolution multi-view images[END_REF] or deep learning refinement strategies [START_REF] Stucker | ResDepth: Learned residual stereo reconstruction[END_REF] to minimize artifacts caused by outliers, vegetation or inexact geometry registration.

In awareness of the complexity of the process, the NASA Ames pipeline [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF] has been expanded with the addition of optional processing tools including bundle adjustment and point cloud registration algorithms.

In the context of 3D reconstruction from push-frame images, we eliminate much of the previous pre-and post-processing work by addressing the root of the problem: the format of the data. The input images are combined to simplify the MVS problem into an ideal scenario with a single input pair. The proposed perfect sensor geometry model for push-frame strips takes inspiration from large push-broom satellites such as in the Pléiades-HR constellation. Due to the complexity of the focal plane, the Pléiades-HR raw products should be considered as 25 sub-products with their local geometrical models [START_REF] Baillarin | Pléiades-HR imaging system: Ground processing and products performances, few months before launch[END_REF]. However, for the sake of usability, the final images emulate the geometry of an ideal push-broom linear array. The Pléiades-HR perfect sensor geometry models are derived from the raw image, the rigorous sensor model of the satellite and a coarse elevation model of the area [START_REF] Baillarin | Pléiades-HR imaging system: Ground processing and products performances, few months before launch[END_REF]. Another existing tool that follows a similar philosophy is the dg_mosaic from the Ames pipeline, which can mosaic multiple subscenes derived from the same parent Maxar pushbroom product and create a new RPC model for it [START_REF] Shean | An automated, opensource pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery[END_REF]. In this work we address a more generic problem, as the input scenes do not originally belong to a common image.

METHODOLOGY

Given a push-frame strip S1 of N small footprint images (also known as scenes), we propose a method to generate a single and equivalent large footprint image, denoted S + 1 , along with its perfect sensor localization model. In this work, we assume that the input images are radiometrically calibrated and cloudfree. As shown in Figure 2, the method consists of three main steps, which are detailed in the following subsections.

Correction of camera models

First of all, it is necessary to ensure that the local RPC camera models of all the scenes in S1 are geometrically consistent. That is, 3D points in object space project to corresponding points in each image. Enforcing the consistency will ease the subsequent mosaicing step and it will be fundamental to produce the perfect sensor localization model of the mosaic S + 1 . This is traditionally achieved by bundle adjustment in a multi-view setting. However, since the baseline between the consecutive cameras of the same strip is too small for 3D vision purposes, some additional information is needed. To this end, we employ a secondary strip S2 covering the same area of interest observed by S1 but from another point of view. SkySat stereo or tri-stereo products can provide the secondary strip S2. Given the collection of 2•N scenes in S1 and S2, the bundle adjustment methodology proposed in (Marí et al., 2021a) is applied to perform a relative correction of their local RPC models. Each corrected RPC model results from composing the original projection function with a corrective rotation transform around the estimated camera center. The corrective rotation compensates for the main source of inaccuracies, which is the inexact knowledge of the attitude angles. The reference points used by the bundle adjustment are automatically generated from correspondences of SIFT keypoints [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF].

Note that the camera correction step does not restrict the scope of the presented method to stereo push-frame acquisitions. It is also possible to correct the camera models of S1 without a secondary strip if a basemap with reference DSM or a set of ground control points (GCPs) are available. GCPs are points with known 3D coordinates whose position in the images is also available [START_REF] Grodecki | Block adjustment of highresolution satellite images described by rational polynomials[END_REF].

Image mosaicing

Once the camera models are corrected, each scene Ii in S1 is warped into a common image space using a 2D projective transform Hi. A rough estimate of Hi is first obtained by establishing 2D point correspondences between the i-th scene and the central scene of the strip, i.e. the N 2 -th scene. Instead of using classic feature matching techniques, we take advantage of the corrected RPC models computed in the previous step to initialize Hi (Section 3.2.1). Then, each Hi is refined into Ĥi using an image registration method [START_REF] Baker | Equivalence and efficiency of image alignment algorithms[END_REF]Matthews, 2001, Briandet al., 2018) (Section 3.2.2). Lastly, the output mosaic S + 1 is obtained by averaging all the warped scenes as determined by the corresponding Ĥi. A high-order Spline interpolation is used to perform the warping (we use order 5).

Image warping initialization

For each scene Ii, a regular grid of 10×10 2D points is localized in the 3D space using the i-th RPC localization function Li and the average altitude of the area, havg. The average altitude havg of the area may be only an approximation, for instance a rough estimation can be obtained using the SRTM data [START_REF] Farr | The shuttle radar topography mission[END_REF]. Each grid point is then reprojected into the image space of the central scene of the strip, using the corresponding RPC projection function P N

2

. The reprojection results in a set of 2D correspondences between a point x from each scene Ii and its homologous x, located in the mosaic image space, where the central scene of the strip remains in the center. Equation 1 summarizes the previous procedure:

x = P N 2 (Li(x, havg)) where i ∈ [1, N ].
(1)

The correspondences x ↔ x are then used to fit (using a classic DLT algorithm [START_REF] Hartley | Multiple view geometry in computer vision[END_REF])) homographic transformations Hi ∈ R 3×3 such that x = Hix.

Note that in the case of the strip scenes that do not overlap with the central scene, we simply localize and reproject recursively along the neighboring frames, in the direction of the central scene, until we reach the latter.

Image warping refinement

The correspondences used to compute the transform Hi are inaccurate because the havg value used for the reprojection is not the exact altitude of the points seen in the image. However, the estimated transformations are useful to initialize an image-based registration method. The inverse compositional algorithm (ICA) is used to refine the coefficients of Hi, so that the warps of consecutive frames are precisely aligned. This assures the pixel consistency of the aligned scenes before merging them into the S + 1 mosaic. The benefits of ICA to register push-frame satellite acquisitions has been previously studied by [START_REF] Anger | Fast and accurate multi-Frame super-resolution of satellite images[END_REF][START_REF] Briand | Improvements of the inverse compositional algorithm for parametric motion estimation[END_REF]. Consider two scenes I and I ′ of S1 such that I ′ has to be aligned onto I to construct the mosaic S + 1 . Let H and H ′ be the initial homographies associated to I and I ′ respectively, which are computed as described in Section 3.2.1. Using ICA, we first refine H ′ • H -1 such that the warped version of I ′ using this transformation is perfectly aligned onto I. This defines the refined homography R. Using the refined transform, we define the relative correction factor 3 shows the residual difference after alignment with and without the refinement step. Observe that the refined transformation achieves a much better alignment. In this example, the RMSE without refinement is 19.37 and with refinement 5.71, thus confirming the visual result.

C such that R = H ′ • C. Fig- ure
We then define by recurrence the set of refined transforms Ĥi = Hi • Ci, where Ci corresponds to the composition of all necessary correcting factors from the reference image (the central scene) to the i-th image in S1.

Perfect sensor geometry localization model

After completion of the mosaic S + 1 , the corrected camera models of the scenes that form the mosaic can be used to produce a perfect sensor localization model that follows the RPC standard and is valid throughout the entire S + 1 . The output camera model, denoted RPC + S 1 , is generated by Algorithm 1.

The main idea of Algorithm 1 is to draw a regular grid of 2D points covering S + 1 , which is localized at different heights in the 3D space. By default, we use (N •M )×(3•M ) points, where N is the number of scenes and M = 10. Given a 2D point x ∈ S + 1 , the corresponding 3D point X at height h is obtained as

X = Li( Ĥ-1 i x, h), (2)
where Ĥ-1 i is the inverse warping transform that transforms x to its original small scene space, and Li is the localization function of that scene. We set the range of altitudes [h min , hmax] by taking the maximum and minimum altitudes of the reference points used by the bundle adjustment (Marí et al., 2021a) and adding an extra margin of ±100 meters.

Thanks to the corrected camera models, the different local RPC functions are highly consistent in the object space. This implies that the localization of the grid will result in a reasonably regular point cloud in the object space, without major discontinuities due to RPC inaccuracies. Furthermore, when a point x is seen in two overlapping scenes Ii and Ij, it can be localized in object space, with (2), using either Li or Lj. Since the scenes for each scene Ii in S1 Localize at h the non-visited points of G2D seen in scene Ii, using (2), and mark them as visited 3. Use the

(N•M )×(3•M )×(3•M ) correspondences
G2D ↔ G3D as input for [START_REF] Akiki | Robust rational polynomial camera modelling for SAR and pushbroom imaging[END_REF] to fit RPC + S 1 are registered and the RPC models adjusted, both choices will yield essentially the same 3D points for a range of altitudes centered around the surface. However, for points far from the surface we should start to observe a parallax due to the fact that the scenes are acquired from different positions along the orbit. This parallax can be quantified. Assuming the SkySat parameters (baseline between views ∼1.5 km and altitude ∼500 km), in order to observe parallax of 1 pixel (assuming a resolution of 0.6 m/pixel) an elevation change of about 200 m ≈ 500 1.5 • 0.6 m should be present in the scene. This points to a limitation of the present method for images containing very large elevation changes. However, in our experiments (including the mountainous site of the Morenci mine, see Section 4) we did not observe any artifacts due to parallax.

The previous procedure results into a set of 2D-to-3D point correspondences between S + 1 and the object space, thus the RPC fitting algorithm from [START_REF] Akiki | Robust rational polynomial camera modelling for SAR and pushbroom imaging[END_REF] can be applied to produce the final RPC + S 1 model.

Using the 2D-to-3D correspondences generated earlier with (2) we can define the fitting errors e associated to the RPC + S 1 model as the reprojection distances:

e = ∥ x -P S + 1 ( X)∥2, ( 3
)
where e is in pixel units and P S + 1 is the projection function of the perfect sensor localization model. Figure 4(b) illustrates the usual distribution of errors e across the 3D points used to fit RPC + S 1 . We observe that the error is small in the proximity of the surface, which can be inferred from the reference keypoints (seen as blue dots) used in the bundle adjustment (Section 3.1). This is reasonable as the surface points are registered in the merged product. Larger errors are observed approaching the altitude extrema of the volume, but only in bands that correspond to the overlap of two consecutive scenes. We attribute this to the fact that there is no guarantee that the RPCs of neighboring scenes are geometrically consistent away from the registered surface points. Inside the convex hull that contains all the reference points, e reaches average values ∼0.2 pixels, of the same order as the average reprojection error of the bundle adjustment (Marí et al., 2021a). Note that the 3D points used by the bundle adjustment highlight the part of the volume where the surface observed by S + 1 is located.

EXPERIMENTS

Data

We applied our method to two SkySat L1B stereo acquisitions, each providing two multi-image strips of the same area, with a time difference of a few seconds. One acquisition covers part of the city of Antibes (France) and the other covers the Morenci mine (United States). The two landscapes are very different: urbanized and flat terrain in the case of Antibes, as opposed to the mountainous and bare terrain of the mine.

SkySat L1B scenes have a nadir resolution of 0.58-0.86 m/pixel and a total size of 1349×3199 pixels. Each scene is delivered with an RPC camera model. The geometric accuracy of the L1B RPC models is of 30-50 m, with SkySats orbiting at altitudes of 400-500 km [START_REF] Planet | Planet imagery product specifications[END_REF]. Our proposed L1B + mosaics extend the footprint of the original L1B images and incorporate a consistent RPC model. In this paper, we present experiments using N = 3 and N = 5 scenes per strip but the method can generalize to strips with more scenes. Using SkySat L1B scenes, we observed no deformation or misalignment in the output mosaics with strips with a number of scenes up to N = 13.

The SkySat acquisition platform has three staggered sensors, resulting in the push-frame system simultaneously acquiring three multi-image strips. Note that in this work we assemble images from only one of the sensors at a time.

Stereo reconstruction based evaluation

To validate the quality of the L1B + images and their camera models, we evaluate them in the context of stereo reconstruction from two push-frame strips S1 and S2. For this purpose, we used the open-source satellite stereo pipeline S2P1 (de Franchis et al., 2014a), to reconstruct the areas of Antibes and the Morenci mine covered by the SkySat acquisitions, both using the original L1B scenes and the L1B + mosaics, i.e. S + 1 and S + 2 .

As explained in Section 2, for the case of the L1B scenes the 3D reconstruction is a multi-view problem. We use the MVS methodology described in (Marí et al., 2021b) to solve it. Following the selection of P suitable pairs of scenes, S2P is employed to reconstruct P independent local DSMs, using the corrected RPCs of the L1B scenes (Section 3.1). The P local models are lastly fused by taking the mean altitude at each cell of the DSM. In the conducted experiments, P = 5 for N = 3, while P = 9 for N = 5.

In contrast with the above, the L1B + products allow to reconstruct each area with a single execution of S2P, using as input the two perfect sensor images and their localization models.

Discussion

Table 1 lists the mean absolute error (MAE) between the DSMs obtained with the L1B + and L1B products. In addition, we computed the MAE of each model with respect to a groundtruth (GT) lidar model covering a subregion of the observed areas. Since the acquisition dates of the lidar and the SkySat images are not coincident, we manually annotated the parts of the surface models that are expected to be coherent. Figure 5 shows the L1B + images for N = 5, the resulting L1B + DSMs and the absolute difference with respect to the equivalent L1B DSMs. In Figure 5(c), we can see that the absolute difference between L1B and L1B + DSMs is below 0.3 m in the majority of the surface points (the average corresponds to the L1B + -L1B column of Table 1). However, there are parts of the area where this difference increases and approaches values close to 1 m. These traces are a consequence of the fusion of local models that is needed to generate the L1B DSM. In fact, the traces coincide in great measure with the green boundaries in Figure 5(b), which indicate the areas of overlap between the local models used to produce the L1B DSM. The local models do not match perfectly, because the correction of the camera models registers their geometries with an accuracy < 1 m, but a residual remains (Marí et al., 2021a). Consequently, the average altitude retained by the fusion process is subject to a certain degree of bias, especially in these overlapping zones. In the case of Antibes, the areas showing the largest differences follow a pattern of horizontal stripes, because the local geometries consist of overlapping planes (flat terrain), which are stacked along the vertical axis. In the case of the Morenci mine the pattern is more irregular, with curves caused by non-exactly coincident mountain shapes and peaks due to the presence of outliers near the open pit (close to the upper left corner).

In Figure 6, we selected two subregions where the altitude differences between the L1B + and the L1B DSMs exhibit a strong increase. In accordance with the above observations, we can see that such differences are indeed caused by biases in the L1B DSM. By using L1B + products we eliminate the cause of such biases, i.e. the need to register and merge any local models, so the discontinuities disappear in a natural way.

The last two columns of Table 1 indicate that the L1B + DSMs improve the accuracy of the L1B ones, as they exhibit smaller differences with respect to the lidar. The MAE values obtained with L1B + are quite stable too, regardless of whether 3 or 5 images per strip are used. We attribute to the aforementioned misalignment between local geometries the fact that the MAE of L1B DSMs with respect to lidar tends to be slightly larger and more irregular.

Lastly, Table 1 shows that the MAE with respect to the lidar is higher for Antibes. This last observation is mainly explained by the presence of vegetation and the edges of the buildings, which are not as sharp in the photogrammetric DSM. Both vegetation and buildings are absent in the Morenci landscape.

CONCLUSION

We have presented a method to generate large scale images from fragmented push-frame satellite acquisitions. A perfect sensor localization model is generated for the output images. We denote the resulting product L1B + . The method is validated using two SkySat stereo acquisitions of L1B scenes.

The use of L1B + offers several advantages over SkySat L1B scenes. In this paper we focused on the benefits for 3D reconstruction, which becomes significantly faster and simpler. We avoid the need to handle multiple stereo pairs and to merge a collection of local models, a common drawback in 3D reconstruction from push-frame imagery. The L1B + products make it possible to reconstruct areas of interest of several km 2 with a single execution of a satellite stereo reconstruction pipeline. We also notice accuracy improvements in the L1B + derived DSMs, mainly due to the disappearance of any biases caused by the fusion of local models that is necessary with the L1B scenes.

Future work will focus on extending this methodology to assemble images from all three SkySat sensors at once, combining scenes from three multi-image strips instead of one. The generalization capability of the method should also be investigated using push-frame acquisitions from other satellites.

Figure 2 .

 2 Figure 2. Diagram of the presented methodology.

Figure 3 .

 3 Figure 3. Residual difference in the overlap region after the alignment of two consecutive SkySat L1B scenes of the same push-frame strip. Top to bottom: using the initial transform Hi, using the refined transform Ĥi. The same scaling and colormap have been used for both residuals.

Algorithm 1 :

 1 Perfect sensor localization model generation Input : 1 strip S1 of N partially overlapping scenes Ii, with their local RPCi camera models (i ∈ [1, N ]). S + 1 , the mosaic image of S1, and the matrices Ĥi that warp scene Ii onto S + 1 Output: 1 RPC + S 1 camera model covering the entire S + 1 1. Build G2D, in the image space: Draw a regular grid of (N•M )×(3•M ) 2D points on S + 1 2. Build G3D, in the object space: Mark all the points of G2D as non-visited for h in 3•M uniformly spaced altitudes ∈ [hmin, hmax]

Figure 4 .

 4 Figure 4. (a) 3D grid used for RPC fitting, for a strip of 3 scenes. The z coordinate is the altitude of the points in meters, while x and y correspond to their projection in the image plane. Point colors depend on the scene of the strip that was used to localize each 3D point. (b) Error of the perfect sensor localization model across the 3D grid, in pixel units. The reference points used by the bundle adjustment (Marí et al., 2021a) in the prior correction of the local RPCs are shown in blue. (c) Front view of the error distribution.

Figure 5 .

 5 Figure 5. Qualitative results for SkySat stereo reconstruction, using 5 scenes from the input strips, i.e. N = 5. Left to right: (a) L1B + images. (b) L1B + derived DSMs. The size of the reconstructed area is indicated in square kilometers. The red boundary delimits the region for which a lidar model is available. The green boundaries delimit the regions of overlap between two local models used to derive the equivalent L1B DSM. The double-headed arrow indicates the length covered in the experiments with N = 3. (c) Absolute difference between the L1B + and L1B DSMs. The black rectangular outline delimits a subregion of interest inspected in Figure 6. MAE [m]

Figure 6 .

 6 Figure 6. Detailed view of the L1B + DSMs and a subregion of interest. The inspection of the subregions shows that, in absence of outliers, the largest differences between L1B + and L1B DSMs coincide with small altitude discontinuities (circled in green) in the L1B model. The colormap is different in the subregion images with respect to the complete DSM to improve the contrast between local altitude values. The L1B + -L1B images represent the absolute difference between altitudes, with the same colormap of Figure 5(c).

Table 1 .

 1 Quantitative results of the stereo reconstruction based evaluation, using 3 or 5 scenes from the input strips. Left to right: MAE, in meters, between the L1B + and L1B derived DSMs, and MAE of each DSM with respect to a GT lidar model.

		N	L1B + -L1B L1B + -GT L1B-GT
	s107_20210705T131230Z	3	0.325	1.214	1.270
	s107_20210705T131300Z
	(Antibes)	5	0.331	1.204	1.248
	s4_20190127T175119Z	3	0.391	0.781	0.854
	s4_20190127T175154Z
	(Morenci)	5	0.395	0.745	0.780

https://github.com/centreborelli/s2p

ACKNOWLEDGEMENTS

This work was supported by a grant from Région Île-de-France. It was also partly financed by Office of Naval research grant N00014-17-1-2552, MENRT, and Kayrros. This work was performed using HPC resources from GENCI-IDRIS (grants 2021-AD011012453 and 2022-AD011011801R1). We thank Planet for providing the L1B SkySat images.


Elise Vissenaekens

Katell Guizien
email: guizien@obs-banyuls.fr

Xavier Durrieu De Madron

Ivane Pairaud

Yann Lerredde

François Pere Puig

Yann Leredde

François Bourrin

Accuracy of high resolution coastal flow speed simulations during and outside of wind, wave and stratification events (Gulf of Lion, NW Mediterranean)

Keywords: Coastal circulation, modelling, flow speed, uncertainty quantification, Mediterranean, Gulf of Lion

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Ocean currents are the key drivers of dissolved and particulate compound transport. At the global scale, the thermohaline circulation regulates the earth's climate [START_REF] Mccarthy | Ocean impact on decadal Atlantic climate variability revealed by sea-level observations[END_REF][START_REF] Clark | The role of the thermohaline circulation in abrupt climate change[END_REF]. Wind-driven, upwelling currents arrange nutrient transport and mixing and regulate primary production at the regional scale [START_REF] Falkowski | Biogeochemical Controls and Feedbacks on Ocean Primary Production[END_REF]). From regional to coastal scales, ocean currents play an imperative role in sediment transport and pollution diffusion [START_REF] James | Modelling pollution dispersion, the ecosystem and water quality in coastal waters: a review[END_REF]Dufois et al., 2008;Warner et al., 2008;[START_REF] Mansui | Predicting marine litter accumulation patterns in the Mediterranean basin: Spatio-temporal variability and comparison with empirical data[END_REF]. At all spatial scales, vessel navigation and marine population connectivity (from large mammal migration to benthic species' larval dispersal) are affected by ocean currents [START_REF] Cowen | Connectivity of Marine Populations: Open or Closed?[END_REF][START_REF] Briton | Highresolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF][START_REF] Putman | Marine migrations[END_REF][START_REF] Mannarini | VISIR-I.b: waves and ocean currents for energy efficient navigation[END_REF]. These applications are currently simulated with Lagrangian dispersal models which, in contrast to Eulerian models, disregard mixing processes and only account for transport processes. Unfortunately, ocean velocity observations, which are necessary to describe these transport processes, are often limited in either time or space. Satellite-mounted altimeters and radars, land-based radars and Lagrangian drifters can measure the currents over a wide area, but only near the ocean's surface [START_REF] Dohan | Measuring the Global Ocean Surface Circulation with Satellite and In Situ Observations[END_REF][START_REF] Mader | The European HF Radar inventory[END_REF]. Some in situ current meters do provide flow measurement time series along vertical profiles (e.g. Acoustic Doppler Current Profiler, ADCP), but single point measurements are still common [START_REF] Schroeder | Long-term monitoring programme of the hydrological variability in the Mediterranean Sea: a first overview of the HYDROCHANGES network[END_REF]Durrieu De Madron et al., 2019). ADCPs which were previously only deployed at fixed moorings [START_REF] Guizien | Internal tide generation at a shelf break by an oblique barotropic tide: Observations and analytical modeling[END_REF] are now being mounted on the hulls of ships (Système Acquisition Validation Exploitation de Données des Navires de l'INSU -Projet SAVED https://sextant.ifremer.fr/record/6f6e95e9-8e97-48d6-b536-b40f2ad87402/, accessed 04/06/2021) or on autonomous underwater vehicles [START_REF] Dohan | Measuring the Global Ocean Surface Circulation with Satellite and In Situ Observations[END_REF][START_REF] Bourrin | Glider monitoring of shelf suspended particle dynamics and transport during storm and flooding conditions[END_REF][START_REF] Gentil | Glider-Based Active Acoustic Monitoring of Currents and Turbidity in the Coastal Zone[END_REF].

Ultimately, ocean current measuring devices are either deployed on the horizontal or on the vertical plane, which strongly limits their applicability to study transport processes. For this reason, transport processes are mainly studied using current simulations over the entire ocean. Ocean circulation models vary according to the different scales and processes they aim to simulate. Tide models are bidimensional models, predicting sea surface elevation and depth-integrated horizontal flow transport, whose main application is navigation (Le [START_REF] Provost | How can we improve a global ocean tide model at a regional scale? a test on the yellow sea and the east china sea[END_REF]. Global ocean circulation mod-els (OGCMs) are three-dimensional models resolving the ocean dynamics at coarse spatial scales everywhere on earth (1/12 °). They either rely on atmospheric coupling for climate predictions (Siedler et al., 2001;[START_REF] Chassignet | The HYCOM (HYbrid Coordinate Ocean Model) data assimilative system[END_REF][START_REF] Somot | 21st century climate change scenario for the Mediterranean using a coupled atmosphere-ocean regional climate model[END_REF] or on one-way atmospheric forcing for modelling ocean energy, fishery management and ship routing [START_REF] Drévillon | Learning about Copernicus Marine Environment Monitoring Service "CMEMS": A Practical Introduction to the Use of the European Operational Oceanography Service[END_REF]. Coastal circulation models are three-dimensional models forced by atmospheric models, most of the time without air-sea interaction, simulating the ocean flow dynamics and hydrology on a limited area. These models aim to simulate meso-scale to sub-meso-scale ocean processes, like eddies [START_REF] Hu | Numerical simulations of eddies in the Gulf of Lion[END_REF][START_REF] Hu | Study of a mesoscale anticyclonic eddy in the western part of the Gulf of Lion[END_REF], dense water cascading [START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF] and river plumes [START_REF] Marsaleix | A numerical study of the formation of the Rhône River plume[END_REF]. They use a spatial resolution that reaches about 100 m in the horizontal and 1 m in the vertical (Dumas and Langlois, 20O9;[START_REF] Briton | Highresolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF]. Such models are considered capable of describing the processes controlling the transport of dissolved and/or particulate matter in a variety of applications (oil spills, land-sea transfer, ecosystem modelling, population connectivity). Regional circulation models have also been coupled to wave models for sediment transport and beach erosion prediction [START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF]Dufois et al., 2008;Warner et al., 2008). Examples of these models are the Model for Applications at Regional Scale (MARS 3D, [START_REF] Lazure | An external-internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS)[END_REF]Dumas and Langlois, 20O9), the COupled Hydrodynamical Ecological model for REgioNal Shelf seas (COHERENS, [START_REF] Drévillon | Learning about Copernicus Marine Environment Monitoring Service "CMEMS": A Practical Introduction to the Use of the European Operational Oceanography Service[END_REF]), the Regional Ocean Modelling System (ROMS, [START_REF] Moore | The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems[END_REF] and SYMPHONIE [START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF](Marsaleix et al., , 2009a)).

However, circulation simulations are subject to various sources of uncertainties, either linked to the model's implementation or to the model's intrinsic assumptions. The model's implementation includes the spatial and temporal resolution of the baroclinic modes and the precision of the forcing data (atmospheric forcing, river runoff, bathymetry and open-boundary forcing). The sensitivity to the grid's spatial resolution [START_REF] Kirtman | Impact of ocean model resolution on CCSM climate simulations[END_REF][START_REF] Kvile | Sensitivity of modelled North Sea cod larvae transport to vertical behaviour, ocean model resolution and interannual variation in ocean dynamics[END_REF][START_REF] Cai | A Numerical Study of Hypoxia in Chesapeake Bay Using an Unstructured Grid Model: Validation and Sensitivity to Bathymetry Representation[END_REF] and to atmospheric and open boundary forcing [START_REF] Kourafalou | Evaluation of Global Ocean Data Assimilation Experiment products on South Florida nested simulations with the Hybrid Coordinate Ocean Model[END_REF] has been thoroughly illustrated. In addition to uncertainties coming from model implementations, uncertainties can come from the model's intrinsic assumptions, such as hydrostaticity, the Boussinesq approximation, the turbulent closure scheme and air-sea interaction. The hydrostatic assumption that the vertical variation of the pressure is dominated by gravity acceleration (resulting in negligible vertical velocities compared to horizontal ones) is not met during wave events [START_REF] Marshall | Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling[END_REF]Zhang et al., 2014). The Boussinesq approximation (density variations can be neglected except in the terms associated with buoyancy forcing) may not be met in the upper stratified ocean, since water density can vary up to 5%, particularly in coastal areas under riverine influence. Therefore, the Boussinesq approximation can cause inaccuracies in the Eulerian simulated velocity of the same magnitude as the water density variation [START_REF] Mcdougall | On Conservation Equations in Oceanography: How Accurate Are Boussinesq Ocean Models[END_REF]. Turbulence closure is also a vital part of any flow dynamics model as it distributes the total flow energy between the turbulent energy resulting from all velocity fluctuations at the subgrid scale and the mean flow [START_REF] Boussinesq | Théorie analytique de la chaleur mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière: Refroidissement et échauffement par rayonnement, conductibilité des tiges, lames et masses cristallines, courants de convection[END_REF][START_REF] Prandtl | Bericht über untersuchungen zur ausgebildeten turbu-lenz[END_REF]. This splitting of the flow energy is essential to describe transport and mixing processes in the numerical simulations. Turbulence closure is expected to play a more prominent role when energetic transfer happens at scales smaller than the spatio-temporal grid, such as during wind-wave [START_REF] Fisher | Wind-Wave Effects on Estuarine Turbulence: A Comparison of Observations and Second-Moment Closure Predictions[END_REF] or river flooding events [START_REF] Reffray | Secondary flows induced by wind forcing in the Rhne region of freshwater influence[END_REF]. Evaluating model accuracy during selected events when the classical assumptions of ocean models aren't met has been frequent practice in the coastal modelling community over the last two decades [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] Reffray | Secondary flows induced by wind forcing in the Rhne region of freshwater influence[END_REF][START_REF] Petrenko | Circulation in a stratified and wind-forced Gulf of Lions, NW Mediterranean Sea: in situ and modeling data[END_REF][START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF][START_REF] Estournel | High resolution modeling of dense water formation in the north-western Mediterranean during winter 2012-2013: Processes and budget[END_REF] the Gulf of Lion). Nevertheless, to disentangle uncertainties due to model assumption violation from those related to implementation, it is necessary to quantify the uncertainty of the model when the assumptions are valid. To our knowledge, this has never been done together and actually, implementation uncertainties on predictied flow speed have been assessed qualitatively only [START_REF] André | Modelled variability of the sea surface circulation in the North-western Mediter-ranean Sea and in the Gulf of Lions[END_REF][START_REF] Petrenko | Circulation in a stratified and wind-forced Gulf of Lions, NW Mediterranean Sea: in situ and modeling data[END_REF][START_REF] Schaeffer | Influence of high-resolution wind forcing on hydrodynamic modeling of the Gulf of Lions[END_REF] the Gulf of Lion).

In the present study, we assessed the uncertainties of regional circulation speed simulations performed in the NW Mediterranean Sea with the hydrostatic Boussinesq model SYMPHONIE (S26 version, https://sirocco.obsmip.fr/ocean-models/s-model/download/), implemented at one of the finest spatio-temporal resolution to date for bathymetry, atmospheric data and river data. The simulations, which were performed from January 2010 to June 2013, were compared to hydrodynamic observations available in the area during this period. Uncertainties in flow speed in different locations and periods were assessed when the model's assumptions were valid (reference period in absence of wind, waves and stratification) and when assumptions were violated (strong wind events, wave events and stratification events).

Model performance was systematically assessed by comparing six indicators calculated during each event type and observation station to their statistical distribution outside of these events.

Material and methods

Study area

The Gulf of Lion is located in the northwestern part of the microtidal Mediterranean Sea and has a wide continental shelf with a mean depth of 70 m [START_REF] Aloisi | Carte géologique du précontinent languedocien au 1/250000ième[END_REF]. It is delineated by a steep shelf break, incised by a dense network of submarine canyons (Figure 1). Its coastal circulation mainly results from the interaction between the thermohaline Northern Current, which flows along the shelf break from the northeast to the southwest and the frequent continental winds blowing from the north and northwest (Mistral and Tramontane resp.), which induce winter convection [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF]. The south-easterly and southerly winds, which blow less frequently, occur mainly from autumn to spring and can cause large swells [START_REF] Guizien | Spatial variability of wave conditions in the gulf of lions (nw mediterranean sea)[END_REF]. The Gulf of Lion's coastal circulation is also influenced by the outflow of one of the largest Mediterranean rivers, the Rhône River, and a series of smaller rivers with typical Mediterranean flash-flooding regimes [START_REF] Guizien | Nearshore dynamics of nutrients and chlorophyll during Mediterranean-type flash-floods[END_REF][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]. The size of the freshwater plume from the Rhône River depends on the atmospheric conditions, the strength of the river flow and the sea water circulation [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF][START_REF] Many | Particle assemblage characterization in the Rhone River ROFI[END_REF][START_REF] Many | Glider and satellite monitoring of the variability of the suspended particle distribution and size in the Rhône ROFI[END_REF]. The surface layers in the Gulf of Lion can stratify thermally between spring and autumn and are recurrently destabilised nearshore by coastal upwelling [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF][START_REF] Petrenko | Circulation in a stratified and wind-forced Gulf of Lions, NW Mediterranean Sea: in situ and modeling data[END_REF].

Water current observations

Horizontal velocity measurements were gathered from eight locations in the Gulf of Lion between January 2010 and June 2013 (Figure 1). Observations included the shallow coastal ADCP moorings BeSete, Mesurho, POEM, and SOLA and the deep moorings Planier, Cap de Creus (Creus), Lacaze-Duthiers (LD) and Lion with one or more single point, acoustic Doppler current meters (SP-ADCMs). The time periods for which flow speed data was acquired are given in Table 1. Additional information on the observations, such as equipment specifications, can be found in the appendix Table A.1.

The observations were filtered to remove erroneous data. For the deep stations, if the velocity measurements presented abnormal values (defined as spikes of intensity with respect to the daily average greater than three times the standard deviation), they were replaced by the average of the previous and the following valid value. For the shallow stations, the upper three meters of the water column were not taken into account, to avoid measuring air speed amid sea surface fluctuations. Moreover, all observations were filtered over time to detect unrealistically fast changes in water speed. The maximum change in water speed tolerated was 30 cm/s over one hour. Another filter was applied on the vertical level and the maximum change in water speed tolerated was 10 cm/s over one meter.

Ocean circulation simulations

The free surface ocean model SYMPHONIE (Marsaleix et al., 2009a(Marsaleix et al., ,b, 2012, SIROCCO, https://sirocco.obs-mip.fr/ocean-models/s-model/, accessed 17/05/2021) was set up to perform regional ocean circulation simulations at a very high resolution in the Gulf of Lion [START_REF] Briton | Highresolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF]. The model solves hydrostatic primitive equations with a finite-difference method on a C curvilinear grid under Boussinesq approximation and with an energy conserving numerical scheme [START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF]. Wave-coupling was not activated and turbulent closure scheme was set to two-equation K-ǫ [START_REF] Michaud | Three-dimensional modelling of waveinduced current from the surf zone to the inner shelf[END_REF]. A bipolar, curvilinear, 680x710 horizontal grid was used to mesh the Gulf of Lion yielding a resolution of 80 m at the coast and 2.7 km in the open ocean (Figure 1, [START_REF] Bentsen | Coordinate Transformation on a Sphere Using Conformal Mapping[END_REF]. Generalized σ-coordinates were used for vertical meshing, with 29 vertical levels [START_REF] Briton | Highresolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF]. Simulations were carried out over the period January 2010-June 2013 and were forced by sea-surface dynamical downscaling of the ERA-Interim atmospheric reanalysis by the regional climate model ALADIN-Climate (ALDERA, 12 km horizontal and 3 h temporal resolutions) and by open-sea boundary conditions from the hindcast downscaled simulation NM12-FREE (∼ 7 km horizontal resolution, [START_REF] Hamon | Design and validation of MEDRYS, a Mediterranean Sea reanalysis over the period 1992-2013[END_REF]. Observed daily discharge of nine rivers (Var, Grand Rhône, Petit Rhône, Hérault, Orb, Aude, Agly, Têt, Tech; http://www.hydro.eaufrance.fr/, accessed 17/05/2021) were included as well. The model's internal and external timesteps were 25.48s and 1.59s, respectively. The simulated velocities were extracted four times per hour on minute 0, 20, 30 and 40 to correspond with the times the observations were measured. On the horizontal, the simulated flow speeds were extracted at the grid point closest to the observations' location (less than 132 m apart).

On the vertical, since the simulation's vertical levels did not match the observations' depths, the simulated speeds were interpolated at the same depth as the observations. If the actual water depth was larger than the water depth in the simulation (bathymetric discrepancy), the simulated speeds were in-terpolated at the depth with the same distance from the bottom as the observation.

Statistical indicators

The deviation between observed O ij and simulated M ij current speed at depth i and time j was described by six time-and depth-averaged statistical indicators, calculated as follows: kg/m³, according to [START_REF] Smith | Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature[END_REF]. Practically, during the reference period, wind stress values should not exceed 0.2586 Pa at both Planier and POEM stations, while wind events were defined by wind stress values larger than 0.4041 Pa during more than 12 hr at both stations. Numerous northerly wind events (37, Figure 2) were detected with wind stresses between 0.6903 Pa and 2.4939 Pa, as expected in the Gulf of Lion [START_REF] Guénard | An observational study of the mesoscale mistral dynamics[END_REF].

Bias = Σ N d i=1 Σ Nt j=1 (M ij -O ij ) N d N t (1)
RelativeBias = Σ N d i=1 Σ Nt j=1 (M ij -O ij ) N d N t Σ N d i=1 Σ Nt j=1 M ij O ij (2)
RMSE = Σ N d i=1 Σ Nt j=1 (M ij -O ij ) 2 N d N t (3) HH = Σ N d i=1 Σ Nt j=1 [(M ij -O ij ) 2 Σ N d i=1 Σ Nt j=1 M ij O ij (4) SI = Σ N d i=1 Σ Nt j=1 [(M ij -M ) -(O ij -Ō)] 2 Σ N d i=1 Σ Nt j=1 M ij O ij (5) Correlation = Σ N d i=1 Σ Nt j=1 [(M ij -M ).(O ij -Ō)] 2 Σ N d i=1 Σ Nt j=1 (M ij -M ) 2 .Σ N d i=1 Σ Nt j=1 (O ij -Ō) 2 (
These events were grouped according to their duration into four different Wave events were defined as the occurrence of swell with a peak period larger than 8 s, a significant period larger than 5 s and a zeroth order moment wave height larger than 3 m at any of four stations during at least 12 hours. The four stations were necessary to detect the southerly to easterly swell impacting the Gulf of Lion [START_REF] Guizien | Spatial variability of wave conditions in the gulf of lions (nw mediterranean sea)[END_REF]. Such swell with wave length larger than the resolution of the flow model at the coast (∼ 100 m) exhibit wave steepness (wave heigth to wave length ratio) larger than 1%, which corresponded to vertical to horizontal velocity ruling out the hydrostatic assumption of the flow model. These criteria resulted in the selection of five swell events with different durations: 12 hr (max.: Tp=12.5 s, Tz=8.0 s, Hm0=4.2 m), 15 hr (max.: Tp=11.8 s, Tz=8.0 s, Hm0=5.5 m), 21 hr (max.: Tp=10.5 s, Tz=8.3 s, Hm0=4.1 m), 40 hr (max.: Tp=10.5 s, Tz=7.8 s, Hm0= 5.6 m) and 86 hr (max.: Tp=10.5 s, Tz=7.7 s, Hm0=4.4 m) (Figure 2). On the contrary, the reference period was defined by the absence of swell with the above mention characteristics at the four stations.

Stratification was estimated at each station after computing the Brunt-Väisälä frequency N 2 (z) = -g/ρ 0 dρ 0 /dz with g the gravitational acceleration, ρ 0 the density of sea water and z the depth in the sea water using simulated salinity and temperature profiles to calculate sea water density [START_REF] Fofonoff | Algorithms for computation of fundamental properties of seawater[END_REF]. The threshold to separate stratified and unstratified periods was the maximum value of the Brunt-Väisälä frequency over the entire water column of 0.005 s-² for at least 12 hours. This value was defined according to [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF]. This allowed us to assess the stratification events at each station separately. No stratification events were detected at the stations SOLA, LD, Lion and Planier, while at Mesurho, which was closest to the Rhône river, the water column was almost always stratified.

Since at the aforementioned stations, there was either an absence of stratification events or of reference conditions, there could be no comparison between the two. Therefore, none of these stations were used for testing the effect of the Boussinesq hypothesis violation on the model's performance. The only stations that were considered were BeSete and POEM (shallow stations)

with four stratification events of 249 hr (max. N²=0.0150 s-²), 81 hr (max. The three types of events were decorrelated and wind events could happen any time in the year (Figure 2). Therefore, reference periods were not separated according to the season.

Assessment of model performance during specific events

Each of the aforementioned indicators is expected to vary with the duration, the moment and the location on which they were calculated, either randomly or systematically. Systematic variation indicates a worse model performance. To test the model's performance under specific conditions (such as strong wind, waves or density stratification), the value of each of these indicators was computed during and in absence of such conditions over a same duration. To compare the events to the reference period, reference cumulative frequency distributions (CFDs) were established for each indicator and each station for the same duration as the event to test. To do so, a set of 200 time periods with an equal event duration as the event to test was randomly selected out of the reference period and used to build this reference CFD for the indicator. These 200 time periods each had unique starting moments, but in the case of stations with a short observation period, overlap is possible. A bootstrap procedure was applied to produce 250 repeats of the reference CFD. Those repeats were used to estimate the most probable reference frequency distribution and a confidence interval around it. The most probable reference CFD for the indicator was thus defined by the 50% quartile (median) of the 250 repeats. For the wave and stratification events, the reference CFDs were calculated using the same duration as the event to test. For the wind events, the reference CFD was calculated over a duration equal to the duration of the middle of the class this event belonged to (e.g A wind event of 14 hr would belong to the class of 12-24hr and be compared to the CFD calculated over 18 hr, as this is the middle of the class, see appendix the event was larger than 75%, it was considered to have a higher uncertainty during the event. For the correlation, the uncertainty of the model is the lowest when the correlation is closer to 1. Therefore, there was a bad model performance when the corresponding cumulative frequency was less than 25%. For the bias, the uncertainty is the lowest when bias is close to zero. Therefore, bad model performance was determined by a corresponding cumulative frequency below 12.5% or above 87.5%. The proportion of events during which the model performed worse than during the reference period was calculated per station, per indicator and per event type. Those proportions were averaged across all indicators and stations to assess whether there was a difference in model performance per station, indicator or event type.

Results

The reference period CFD of each accuracy indicator in absence of wind, waves and stratification was computed for durations ranging from 12 hr to 249 hr at each station (Figure 4 for a duration of 42 hr). Overall, the simulated flow speed was underestimated at deep stations during the reference period, with bias median values calculated over 42 hr ranging from -3 cm/s at Lion to -1.2 cm/s in Creus and LD (Figure 4A). At the shallow stations, the flow speed could be either underestimated (BeSete and SOLA, bias median values of -1.2 cm/s and -0.6 cm/s resp.) or overestimated (Mesurho and POEM, bias median values of 3.6 cm/s and 0.6 cm/s; Figure 4A). In both groups of stations, bias values spread was large, with the first and third quartile being -5.4 and 0 cm/s at deep stations and -4.2 and 6.6 cm/s at shallow stations.

After normalizing by the current magnitude in each station, the relative bias was smaller at the shallow stations (with median values ranging from 30% at BeSete and SOLA to 40% at POEM) than at the deep stations (with median values ranging from 35% at Creus to 85% at Lion, Figure 4B). The relative scatter index (SI) was variable amongst the stations, with a similar variability among deep and shallow stations (median values ranging from 65%

at Lion to 93% at Mesurho; Figure 4C). As a result, the HH indicator, which combines the relative scatter and relative bias was larger at deep stations (median values ranging from 95% at LD and Creus to 110% at Lion) than at shallow stations (median values ranging from 75% at SOLA and BeSete to 83% at Mesurho, Figure 4E). Noteworthy, the median HH values were larger than 70% at all stations. In absolute values, the median RMSE was similar at deep and shallow stations, ranging from 2.5 cm/s at Planier and BeSete to 5.6 cm/s at Creus and 5.3 cm/s at Mesurho (Figure 4D). However, the RMSE's third quartile was less homogenous across deep stations, which had values ranging from 3.2 cm/s to 12.6 cm/s, than across shallow stations, with values ranging from 3.2 cm/s to 6.7 cm/s. Although the correlation was low at all stations, it was higher at the shallow stations than at the deep stations (Figure 4F). Median (third quartile) values ranged from 0.03 (0.14, resp.) at BeSete to 0.13 (0.23, resp.) at POEM while median values in deep stations had a median of -0.01 for LD and Lion and were always less than 0.01. Although the CFDs of the accuracy indicators clustered according to the duration of the event, the deviation between the CFDs calculated over 12-24 hr and those calculated over more than 72 hr remained limited (Appendix to the duration of the events, regardless of the station (Figure 6 for wind events). Overall, the proportion of events where the model performed worse during the events than during the reference period was low no matter the event type. The average ratio worse ranged from 25% for the wind events to 35% for the wave events (Table 2, Table 3). For the stratification events, which were only studied at BeSete and POEM, the model performed worse during the events than during the reference period for 25% and 33% of the events on average, respectively (Table 4). However, the assessment of the model's performance varied greatly depending on the indicator, with the HH indicating a 13% ratio worse and the RMSE showing a 45% ratio worse in the wind events for instance (Table 2, Figure 7B). When comparing the model's performance across event types and stations, it was worse during wave events than during wind events at shallow stations (except at the Mesurho station in front of the Rhone River mouth), while no trend could be observed at deep stations (Figure 7A). During the wave events, the model performed similarly across all stations, with all stations indicating that the model was worse than during the reference period less than 33% of the time on average, except at the POEM station, where the ratio worse reached 67% (Table 3).

During the wind events, the model performed slightly better at the shallow stations (ratio worse ranging from 11% to 31%) than at the deep stations (ratio worse ranging from 26% to 37%, Table 2). For both event types, absolute indicators (RMSE and bias) displayed worse model performance than relative indicators (Figure 7B). All indicators except SI displayed worse model performance during wave events than during wind events (Figure 7B).

Discussion

The present study quantified various indicators to describe the deviation between observed and simulated flow speed across shallow and deep stations within a highly dynamic region, during and outside short term events of three types (wind, waves, stratification).

The assessment of ocean model accuracy has largely been implemented by comparing simulated and observed hydrological variables (temperature and salinity; e.g. [START_REF] Gustafsson | Coupling of a High-Resolution Atmospheric Model and an Ocean Model for the Baltic Sea[END_REF][START_REF] Reffray | Secondary flows induced by wind forcing in the Rhne region of freshwater influence[END_REF][START_REF] André | Modelled variability of the sea surface circulation in the North-western Mediter-ranean Sea and in the Gulf of Lions[END_REF][START_REF] Kara | Validation of interannual simulations from the 1/8°global Navy Coastal Ocean Model (NCOM)[END_REF][START_REF] Chelton | Summertime Coupling between Sea Surface Temperature and Wind Stress in the California Current System[END_REF][START_REF] Pairaud | Hydrology and circulation in a coastal area off Marseille: Validation of a nested 3D model with observations[END_REF][START_REF] Renault | Coupled atmosphere-ocean-wave simulations of a storm event over the Gulf of Lion and Balearic Sea: COUPLED ATMOSPHERE-OCEAN-WAVE SIMULATION[END_REF][START_REF] Marzocchi | The North Atlantic subpolar circulation in an eddy-resolving global ocean model[END_REF][START_REF] Seyfried | Modelling deep-water formation in the north-west Mediterranean Sea with a new air-sea coupled model: sensitivity to turbulent flux parameterizations[END_REF][START_REF] Akhtar | Climate modeling over the Mediterranean Sea: impact of resolution and ocean coupling[END_REF] as their dynamics integrates transport (velocity) and mixing (turbulent kinetic energy) in ocean circulation models. However, hydrological variables are little informative about transport and mixing when well-mixed conditions prevail, which is often the case in coastal areas [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF][START_REF] Holt | Modelling the global coastal ocean[END_REF].

The ability of the SYMPHONIE model to simulate flow speed and not only hydrological parameters in the Gulf of Lion has been assessed before but only qualitatively under a variety of coastal processes, such as thermally stratified conditions [START_REF] Petrenko | Circulation in a stratified and wind-forced Gulf of Lions, NW Mediterranean Sea: in situ and modeling data[END_REF], fresh water mixing in the Rhône River prodelta [START_REF] Estournel | The Rhone River Plume in Unsteady Conditions: Numerical and Experimental Results[END_REF], wind driven Eckman flow [START_REF] Davies | A three-dimensional model of wind-driven circulation on the shelf[END_REF][START_REF] Lapouyade | Seasonal variability of the advective transport of particulate matter and organic carbon in the Gulf of Lion (NW Mediterranean)[END_REF][START_REF] Molcard | Wind driven general circulation of the Mediterranean Sea simulated with a Spectral Element Ocean Model[END_REF][START_REF] Schaeffer | Influence of high-resolution wind forcing on hydrodynamic modeling of the Gulf of Lions[END_REF][START_REF] Estournel | High resolution modeling of dense water formation in the north-western Mediterranean during winter 2012-2013: Processes and budget[END_REF], swell events [START_REF] Michaud | Three-dimensional modelling of waveinduced current from the surf zone to the inner shelf[END_REF][START_REF] Mikolajczak | Impact of storms on residence times and export of coastal waters during a mild autumn/winter period in the Gulf of Lion[END_REF] , 2007). This is generally the the case among the few studies that quantified uncertainties on simulated speeds, elsewhere. While comparing Glazur60 simulations of the NEMO model (horizontal resolution of 1/64°hence 1.3 to 1.7 km) to the data of a fixed ADCP mooring located at a hundred meters depth in the eastern part of the Gulf of Lion, a bias of 3.5 cm/s at 90 m and 7 cm/s at 20 m depth was found between simulations and observations over an 11 month integration period [START_REF] Barrier | Strong intrusions of the Northern Mediterranean Current on the eastern Gulf of Lion: insights from in-situ observations and high resolution numerical modelling[END_REF].

Similarly, while evaluating the effect of boundary conditions on simulations using the SoFLA-HYCOM model configuration (1/25°hence 3.5 to 4 km horizontal resolution) at shallow stations around the Strait of Florida, the mean bias and the RMSE calculated between simulations and observations over a one year period ranged from -3.5 cm/s to 8.2 cm/s for the bias and from 5 to 13 cm/s for the RMSE, depending on the model's configuration and the station [START_REF] Kourafalou | Evaluation of Global Ocean Data Assimilation Experiment products on South Florida nested simulations with the Hybrid Coordinate Ocean Model[END_REF]. Despite flow speed simulations not being as precise as ADCP measurements, it is remarkable that the present study's bias and RMSE values were smaller than the values reported in those quantitative studies, despite these indicators were calculated over longer periods in the latter studies than in our study (weeks versus days). Indeed, the systematic bias and the RMSE are expected to decrease with increasing integration duration [START_REF] Dekking | A modern introduction to probability and statistics: understanding why and how[END_REF]. However, comparing bias and RMSE values between simulations and observations in different environments can be misleading regarding model performance and relative indicators should be used.

In the present study, lower relative bias and HH were found at shallow stations compared to deep ones. The better model performance at the shallow stations could be due to the refinement of the horizontal spatial resolution, thanks to the adaptive resolution of the curvilinear grid. Increasing the resolution of model configurations have been tested to improve agreement with other types of observations than flow speeds, sometimes showing predictions improvements [START_REF] Thoppil | Energetics of a global ocean circulation model compared to observations: EN-ERGETICS OF THE GLOBAL OCEAN[END_REF][START_REF] Kirtman | Impact of ocean model resolution on CCSM climate simulations[END_REF][START_REF] Putman | Tracking the long-distance dispersal of marine organisms: sensitivity to ocean model resolution[END_REF][START_REF] Ringler | A multi-resolution approach to global ocean modeling[END_REF][START_REF] Akhtar | Climate modeling over the Mediterranean Sea: impact of resolution and ocean coupling[END_REF][START_REF] Kvile | Sensitivity of modelled North Sea cod larvae transport to vertical behaviour, ocean model resolution and interannual variation in ocean dynamics[END_REF][START_REF] Ridenour | Sensitivity of freshwater dynamics to ocean model resolution and river discharge forcing in the Hudson Bay Complex[END_REF]. In addition to relative bias which indicates goodness of transport predictions, the present study evaluated the correlation between simulated and observed flow speed, an indicator generally disregarded. At all stations, correlation indicated that the simulation failed to reproduce the short term flow dynamics (hours to days). Short term flow dynamics is expected to be driven by atmospherical forcings, especially in the Gulf of Lion, where coastal circulation simulations have been shown to dramatically change with the wind's spatial gradient (Dumas and Langlois,20O9). Hence, the present study simulations were driven by atmospheric field outputs from a reanalysis with assimilated observations and was updated every three hours at the finest resolution available for the area at the time of the simulations [START_REF] Hamon | Design and validation of MEDRYS, a Mediterranean Sea reanalysis over the period 1992-2013[END_REF]. One way to improve the simulations' accuracy is to use the bidirectional atmospheric coupling technique [START_REF] Gustafsson | Coupling of a High-Resolution Atmospheric Model and an Ocean Model for the Baltic Sea[END_REF][START_REF] Chelton | Summertime Coupling between Sea Surface Temperature and Wind Stress in the California Current System[END_REF][START_REF] Schaeffer | Influence of high-resolution wind forcing on hydrodynamic modeling of the Gulf of Lions[END_REF][START_REF] Akhtar | Climate modeling over the Mediterranean Sea: impact of resolution and ocean coupling[END_REF]. Two-ways air-sea coupling performed better than one-way atmospheric forcing during autumn storms, when the sea surface cools rapidly [START_REF] Seyfried | Modelling deep-water formation in the north-west Mediterranean Sea with a new air-sea coupled model: sensitivity to turbulent flux parameterizations[END_REF]. Nevertheless, in the Gulf of Lion, the added value of coupling atmosphere-ocean simulations on modelled wind speed intensity and sea surface temperature was not significant [START_REF] Renault | Coupled atmosphere-ocean-wave simulations of a storm event over the Gulf of Lion and Balearic Sea: COUPLED ATMOSPHERE-OCEAN-WAVE SIMULATION[END_REF]. Interestingly, in the present study, the indicators did not display a worse model performance during strong wind events when atmosphere-ocean interaction increased, than outside those events. In any case, this limitation to reproduce the short term flow dynamics, including in low wind conditions, raises the question of how short term (days) velocity dynamics' inaccuracies alter particle tracking simulations (e.g. used in larval dispersal studies, [START_REF] Briton | Highresolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF].

Similarly, the model's performance was not systematically worse during wave events, although it was slightly worse during wave than during wind events. When comparing a hydrostatic, quasi-hydrostatic and nonhydrostatic model, no difference between the three models was found at large scales with coarse resolution (1°horizontal resolution, [START_REF] Marshall | Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling[END_REF]. However, it is expected that quasi-hydrostatic and nonhydrostatic models should be preferred when the spatial resolution increases as in the present study simulations [START_REF] Magaldi | Hydrostatic and non-hydrostatic simulations of dense waters cascading off a shelf: The East Greenland case[END_REF]. Incorporating the effects of waves on the coastal circulation simulations has been considered previously in the Gulf of Lion and flow speed simulations in the surf zone (0-15m water depth)

were improved by using a fully nonhydrostatic coupled current-wave model [START_REF] Michaud | Three-dimensional modelling of waveinduced current from the surf zone to the inner shelf[END_REF]. However, outside the surfzone, deviations between observed and simulated flow speeds at POEM and Mesurho (same location as in the present study but another time period) were similar regardless of wave forcing.

Another model assumption which could have altered the model's performance is the Boussinesq approximation, which can be violated in thermal or fresh water stratification. In the Gulf of Lion, stratification effect was only studied qualitatively. During summer, incorrect representation of the stratification in the Gulf of Lion led to a misplacement of the NC in the simulations compared to the field observations [START_REF] Petrenko | Circulation in a stratified and wind-forced Gulf of Lions, NW Mediterranean Sea: in situ and modeling data[END_REF]. In contrast, simulations of the Rhône plume compared to radar observations

showed that the SYMPHONIE model can reproduce the spatial variation of the current in front of the river mouth outside of strong wind events [START_REF] Estournel | The Rhone River Plume in Unsteady Conditions: Numerical and Experimental Results[END_REF]. Comparing with the rare quantitative studies from other areas is equally unconclusive as only absolute indicators were computed (bias (4-15cm/s) and RMSE (6-18cm/s) over two week period of salinity stratification in an estuary in the USA Yang and Khangaonkar, 2009). In the present study, testing model performance alteration due to stratification was limited to few fresh water input events in two stations only as in other stations, the water column was either never or always stratified. In these few events, model performance was not significantly worse. However, outside specific events, indicators were systematically larger at the continuously stratified

Mesurho station than at the other shallow stations, suggesting stratification effect should be further tested.

In conclusion, a quantitative validation of simulated current speeds was performed over a three-year period using in situ flow speed observations from eight fixed moorings (four shallow and four deep). Multiple absolute, and more importantly, relative indicators were calculated to evaluate the performance of the model. In absence of wind, wave or stratification events, the model performed better at shallow stations than at deep stations in predicting the mean flow speed (lower relative bias). In contrast, scatter index was equally large at all stations and correlation over short duration periods was always low, indicating discrepancies between simulated and observed flow speed dynamics. Overall, the model did not perform notably worse during wind, wave or stratified events than outside of events. However, the model's performance was lower during wave events than during wind events at shallow stations.

Acknowledgements

Some current observations were obtained within the framework of the Overview of all ocean current observation stations with the number of bins, the bin size, the depth, the time step, the type of equipment, the latitude, the longitude and the source. For the ADCPs, only the maximum depth is indicated (*).

  6) with Ō the observed current speed averaged over depth and time and M the simulated current speed averaged over depth and time. The bias (equation 1) is the difference between the simulated and observed mean. It indicates systematic under-(negative value) or overestimation (positive value) of the simulated flow speed. The relative bias (equation 2) is the absolute bias normalized by the square root of the mean of the product of observed and simulated flow speed. The root mean square error (RMSE, equation 3) is the square root of the quadratic mean of differences between simulated and observed velocities. It adds to the bias as a measure of random deviation and indicates the accuracy of simulations. The Hanna & Heinold index (HH, equation 4, Hanna and D., 1985) normalized the RMSE by the mean of the product of the observed and simulated flow speed. It indicates the relative uncertainty from the mean flow and avoids biasing when the model underestimates the currents (negative bias, Mentaschi et al., 2013). The scatter index (SI, equation 5) is the quadratic mean of the difference between simulated and observed flow speed fluctuations around the mean, normalized by the mean of the product of observed and simulated flow speed. It indicates if the simulated flow speed fluctuates more or less around the mean than the observed flow speed. The correlation index (equation 6) is the product of simulated and observed fluctuations around the mean flow speed, normalized by the product of the standard deviation of the simulated and observed flow speed. It varies between -1 and +1. Values close to 1 indicate co-variation (-1 indicates opposed variation) in the dynamics of simulated and observed flow speed, while values close to 0 indicate the dynamics of simulated and observed flow are different. 2.5. Definition of the reference period and the three types of specific events In order to separate uncertainties due to model implementation and hypothesis violation, the six indicators were assessed separately during reference periods defined by low wind conditions, no swell and absence of significant stratification, and during events with either strong wind conditions (turbulence closure or atmospheric forcing reliability), swell (hydrostatic hypothesis violation), or in stratified conditions (Boussinesq approximation violation). Importantly, the indicators were integrated over a same duration during reference periods as the event duration. Wind conditions over the entire Gulf of Lion were estimated using the wind stress used to force the ocean circulation simulations at the closest atmospheric model grid point from the Planier and POEM stations (Figure 1). Low and strong wind conditions correspond to wind speed lower than 40 km/hr and larger than 50 km/hr, respectively, separating negligible effects from significant impacts in Beaufort scale. Wind speed thresholds were converted into wind stress values to be detected in the atmospheric forcings (using τ = C D ρU 2 with τ the wind stress in Pa, U the wind speed in m.s -1 , C D a drag coefficient of 0.00171 and ρ the air density of 1.225

  N²=0.0084 s-²), 194 hr (max. N²= 0.0202 s-²) and 143 hr (max. N²=0.0124 s-²) at BeSete and three stratification events of 74 hr (max. N²=0.0953 s-²), 79 hr (max. N²=0.0310 s-²) and 103 hr (max. N²=0.0211 s-²) at POEM (Figure 2).

Figure A. 2

 2 Figure A.2 for more information. Reference CFDs were used to determine the corresponding cumulative frequencies of each indicator/station/event by assessing the event's indicator value compared to the reference CFD (Figure 3, additional schematic in appendix Figure A.2). Those corresponding cumulative frequencies were used to assess the model's performance, by comparing its value to a threshold value. For RMSE, HH and SI and relative

Figure A. 3

 3 Figure A.3). The median correlation at BeSete varied between 0.02 and 0.06

  and dense water cascading[START_REF] Marsaleix | Energy conservation issues in sigmacoordinate free-surface ocean models[END_REF][START_REF] Estournel | High resolution modeling of dense water formation in the north-western Mediterranean during winter 2012-2013: Processes and budget[END_REF]. Only one study assessed quantitatively the uncertainty on simulated speeds in the Gulf of Lion. It compared another SYM-PHONIE configuration than the one of the present study (horizontal resolution ranging from 300 m to 7 km, with and without wave coupling) to part of the dataset used in our study that is a two month period which included several wave events in (February toMarch, 2011[START_REF] Mikolajczak | Dynamique de l'eau et des apports particulaires originaires du Rhône sur la marge continentale du Golfe du Lion[END_REF]. The bias was 4 cm/s at the Mesurho station and -4 cm/s at the POEM station whilst the RMSEs were 10 cm/s and 8 cm/s, respectively. The present study compliments previous assessments of the SYMPHONIE model in the Gulf of Lion, whilst extending them in space and time and using six quantitative indicators. Using data from multiple years and stations, particularly shallow versus deep ones enabled us to assess model implementation uncertainties in the present study configuration. Focussing first on a reference period (unstratified, with low wind conditions and no wave), when model assumptions are expected to be valid, bias and RMSE on simulated speeds during the reference period were larger than the measuring device accuracy (about 1 cm/s Instruments

Figure 1 :

 1 Figure 1: The Gulf of Lion. Main bathymetrical contours (20, 50, 200, 1000 m) of the Gulf of Lion including the dipolar model grid (680 × 710; with one blue line every 10 cells; North pole (44.2°N, 5.3°E); South pole (42.37°N, 2.82°E); grid point (170; 710) corresponding to (47°N, S°E); and the reference latitude for Mercator projection was 52°N). Further information on the grid can be found in Briton et al., (2018). The locations of the fixed moorings with current meters are in red: BeSete, Creus, LD (Lacaze-Duthiers), Lion, Mesurho (Measuring buoy at the mouth of the Rhône River), Planier, POEM (Observational Platform of the Mediterranean Environment/Plateforme d'Observation de l'Environnement Méditerranéen), SOLA (SOMLIT Observatory of the Arago Laboratory/SOMLIT Observatoire de Laboratoire Arago) and with wave buoy in black: Banyuls, Espiguette, Leucate, Sète.

Figure 2 :

 2 Figure2: Timetable with selected events (Black). Green is the reference period. For the wind events, the white zones are zones with intermediate wind. The wind and wave events are common to all stations. For the stratification event, striped line (--) is the reference period for Besete and the full line (-) is the reference period for POEM. In the white zone, no observational data was available for these two stations. The dashed vertical lines (:) indicate the seasons and the letter triplets are the first letters of the months in that season.

Figure 3 :

 3 Figure 3: Corresponding cumulative frequency example. The corresponding cumulative frequency of the indicator value during the event can be read on the y-axis of when placing the indicator value calculated during the event (orange X) on the cumulative frequency of the indicator values during the reference period (blue line).

Figure 4 :

 4 Figure 4: 42 hr reference period CFD. The indicators' cumulative frequencies integrated over 42 hr at all stations during the reference period. Shallow stations are depicted with a dashed line, deep stations with a solid line. A) Bias, B) Relative bias, C) SI, D) RMSE, E) HH, F) Correlation.

Figure 5 :

 5 Figure 5: CFD of the correlation between modelled and observed flow speeds at BeSete during the reference period for different durations. Blue: 12-24hr, green: 24-72hr and grey: more than 72 hr.

Figure 6 :

 6 Figure 6: Corresponding cumulative frequency of the wind's correlation. Scatter plot of the wind event duration in relation to the corresponding cumulative frequency of the correlation between modelled and observed current speed. Events with a corresponding cumulative frequency below 25% are considered worse during the event than during the reference period.

Figure 7 :

 7 Figure 7: Mean percent worse per station and indicator for wind and wave events. Histograms of the mean percent of wind/wave events worse during the events than during the reference period. A) Per station, B) Per indicator.

Figure A. 1 :

 1 Figure A.1: Frequency histogram of the durations of the wind events.

Figure A. 2 :

 2 Figure A.2: Scheme on how to compare the uncertainty of the model during the event to the uncertainty of the model outside of the events.

Figure A. 3 :

 3 Figure A.3: Cumulative frequency distribution of the indicators calculated between modelled and observed flow speeds at BeSete during the reference period for different durations. Blue:12-24hr, green: 24-72hr and grey: more than 72hr.

  classes(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)[START_REF] Drévillon | Learning about Copernicus Marine Environment Monitoring Service "CMEMS": A Practical Introduction to the Use of the European Operational Oceanography Service[END_REF](24)(24)(25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48)(48)(49)(50)(51)(52)(53)(54)(55)(56)(57)(58)(59)(60) appendix Figure A.1).

	Wave conditions over the entire Gulf of Lion were assessed using obser-
	vations over the period January 2010-June 2013 at four stations (Banyuls,
	Espiguette, Leucate, Sète, Figure 1) of the In Situ National Data Archiving
	Center of Waves (Centre d'Archivage National des Données de Houle In Situ,
	http://candhis.cetmef.developpement-durable.gouv.fr, accessed 01/06/2021).

Table 2 :

 2 Proportion of events worse during the event than during the reference period per indicator and per station for wind event type.

	Station	Creus LD Lion Planier Mesurho POEM SOLA BeSete	Mean
	/									per
	Indicator									indicator
	RMSE	67	40	40	20	67	100	0	100	54
	Bias	67	40	40	60	33	100	0	0	43
	HH	0	20	20	20	0	100	0	0	20
	SI	0	20	0	40	0	0	100	0	20
	Relative bias	0	20	20	20	33	100	100	0	37
	Correlation	0	40	60	40	33	0	0	100	34
	Mean	22	30	30	33	28	67	33	33	35
	per station
	Nr of events	3	5	5	5	3	1	1	1

Table 3 :

 3 Proportion of events worse during the event than during the reference period per indicator and per station for wave event type.

	Station	POEM BeSete	Mean
	/			per
	Indicator			indicator
	RMSE	33	25	29
	Bias	33	50	42
	HH	33	25	29
	SI	67	25	46
	Relative bias	0	0	0
	Correlation	33	25	29
	Mean	33	25	29
	per station
	Nr of events	3	4

Table 4 :

 4 Proportion of events worse during the event than during the reference period per indicator and per station for stratification event type.

Table A . 1 :

 A1

Creus X X X X X X X X X X X X X X X X X X X X X X X X LD X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Lion X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Planier X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X Mesurho X X X X X X X X X X X X X X X X X X X X POEM X X X SOLA X X X BeSete X X X X X X X X X X X X X X X X X X X X X