Dynamic Coastal-Shelf Seascapes to Support Marine Policies Using Operational Coastal Oceanography: The French Example

: In the context of maritime spatial planning and the implementation of spatialized Good Environmental Status indicators in the Marine Strategy Framework Directive (MSFD), the deﬁnition of a mosaic composed of coherent and standardised spatial units is necessary. We propose here a characterization of seascapes in time and space within the speciﬁc framework of the MSFD in the English Channel and the Bay of Biscay areas. A spatio-temporal classiﬁcation of coastal-shelf water masses is carried out using twelve essential oceanographic and derived variables from operational coastal oceanography using the HYCOM model. Partitioning is computed using a multivariate hybrid two-step clustering process deﬁning a time series of categorical maps representing hydrographical patch classes. Main patch occurrence is analyzed to understand their spatio-temporal dynamics and their oceanographic characteristics. Finally, patch classes are combined with MSFD marine sub-region delimitations to build seascapes, including ecosystem approach management and marine policy considerations. 4-1 and water mass


Introduction
The sustainable management of maritime spaces requires the establishment of ecosystem-based indicators. Such indicators necessitate, before all, providing managers with inventories and maps of marine biodiversity and habitats, compatible with management obligations and human activity planning, at appropriate spatial and temporal resolutions. In this context, different mapping tools and concepts already exist at large spatial scales, such as Ecologically or Biologically Significant Areas (EBSAs) adopted by the Convention on Biological Diversity (CBD). EBSAs are characterized by their ecological or biological characteristics. They were described to support and encourage the use of area-based management tools [1].
However, in the context of marine spatial planning, ecosystem-based management framework, involves the consideration of complex interconnected systems (i.e., ecological, social and economic) at different scales of time and space. Consequently, the main challenge for managing Europe's maritime spaces, whether through the Marine Strategy Framework Directive (MSFD, Directive 2008/56/CE), the Water Framework Directive (WFD, Directive 2000/60/EC), the Habitats and birds (HD, Directive 92/43/EEC), maritime spatial planning (MSP, Directive 2014/89/EU) or other directives and regulations, is to be able to find a relevant balance between local and regional management scales. For example, in the specific case of the MSFD, which aims at achieving the Good Environmental Status New developments in operational oceanography now provide access to oceanic features of significant interest for ocean partitioning, mainly through the ability to compute Essential Ocean Variables (EOVs). Following the definition of Essential Climate Variables (ECVs) expressed by [25], EOVs can be construed as physical, chemical, or biological variables that significantly contribute to the characterization of the ocean. The Global Ocean Observing System (GOOS) has defined several EOVs "in order to avoid duplication of effort, across observing platforms and networks, and to adopt common standards for data collection and dissemination to maximize their usefulness" (https: //www.goosocean.org/index.php?option=com_content&view=article&id=14&Itemid=114). These variables are currently used to support work on climate change, biodiversity conservation and sustainable management of marine resources [26][27][28]. By using models that are specific to areas of interest, particularly coastal areas, operational oceanography thus allows the computation of EOVs and variables derived from EOVs which are considered of interest for ecological systems [29][30][31][32][33]. Consequently, the development of capacities in operational oceanography, particularly the availability of circulation models outputs, offers new opportunities for the calculation of EOVs, the integration of dynamics into an ocean partitioning for the ecosystem assessment asked by marine policies [33,34].
The objective here was therefore to partition coastal-shelf pelagic water masses through hybrid clustering methods and Operational Coastal Oceanography (OCO), using the coastal HYbrid Coordinates Ocean Model (HYCOM [35]), EOVs and derived variables, while considering the dynamics of water masses and temporal variability. The final objective of this work was to compute shelf-coastal patch time series included in MSFD French marine sub-regions as scale of reference (Bay of Biscay, Celtic seas and East English channel-north sea) in order to construct relevant seascapes for marine policies and ecosystem based management. This work is organized as follows: Section 2 describes the ocean modeling framework, the EOVs computation, and the statistical process, Section 3 presents the results of the ocean partitioning, pattern characterization, and the definition of seascapes, and lastly, Section 4 discusses the results in the context of marine policies.

Numerical Model
The three-dimensional HYCOM numerical model of the French Service Hydrographique et Océanographique de la Marine (Shom) has been used in this study [36,37]. It solves the primitive equations as described in [35] on a regular Arakawa-C horizontal grid, with a vertical variable hybrid coordinates system which is isopycnal in the inland ocean, geopotential in the surface layers and terrain-following in shallow waters. Several numerical developments have been made to optimize HYCOM for coastal zones [38,39].

Model Set-Up
The numerical model configuration used in this study consisted of a domain covering the North-Eastern Atlantic Ocean region, encompassing the four Atlantic French marine MSFD sub-regions: The Greater North Sea, including the English Channel, the Celtic Sea and the Bay of Biscay (Figure 1). The horizontal grid resolution was set to 1/60°, for a total of 471 × 720 horizontal cells, with a vertical discretization over 40 layers. The bathymetry has been extracted from the Historique, Observation, MOdélisation des NIveaux Marins (HOMONIM [40],) project managed by Météo-France and the Shom.
The forcing current velocities field was derived from both the geostrophic current calculated from the densities extracted from operational PSY4V3R1 products provided by the CMEMS [41] at a spatial resolution of 1/12°and reprojected on the HYCOM grid, and the tidal current calculated from the tidal global barotropic MOG2D 2003 atlas [42] using fourteen tidal constituents (Q1, O1,  P1, K1, 2N2, N2, M2, S2, K2, M3, MN4, M4, MS4 and M6). Vertical turbulent mixing were modeled using the K-Profile Parametrization (KPP) of [43] based on thermodynamic surface forcing and flow characteristics to determine the depth of the surface boundary layer and calculate the diffusivity applied to velocity, temperature and salinity. Bottom friction was distributed over a 50-m thick layer of water using a quadratic law with a constant coefficient equal to 2.5 e −3 . Atmospheric forcings (wind, precipitation, atmospheric pressure and bulk heat fluxes) were derived from the operational forecast model ARPEGE-Climate at 1/10°resolution provided by Météo France [44]. River runoffs (Seine, Loire, Gironde and Adour) were prescribed using the daily measured data provided by the Centre Data in Operational Coastal Oceanography (CDOCO, [45]). Temperatures of fresh waters were calculated as an annual sinusoidal cycle estimated from available historical data. The model was initialized on 1 January 2011. A spin-up of one year has been performed before the five-years long simulation, from 1 January 2012 to 31 December 2016.

Model Outputs: EOVs and Derived Variables
Sea surface temperature, sea surface salinity, sea surface current velocities, sea surface elevation, along with three-dimensional temperature, salinity and density fields were extracted from the model at a one-hour time frequency. Based on these results, twelve variables among EOVs, and other sub-variables (hereinafter referred as subEOVs) of interest have been calculated (Table 1) using a monthly mean approach from January 2012 to December 2016, in order to detect large-and meso-scale physical processes in the regional area. The chosen hydrological and dynamical variables are well known to be of importance for the structuration of marine ecosystems [23,24,32]. Each variable is computed on a monthly basis. This choice is a compromise between the description of high-frequency variations and the computational cost of the statistical process described in Section 2.4.
Hydrological variables have been calculated as in [15]: Sea surface temperature (SST), sea surface salinity (SSS), mixed layer depth (MLD), deficit of potential energy (DPE), SST gradients (GRADSST) and SSS gradients (GRADSSS). SST is used to describe the seasonal cycle in variations of hydrological properties of water masses (Figure 2a). SSS allows to identify and to characterize the variability of river plumes, associated with low salinities (Figure 2b). MLD describes the thickness of the ocean layer exposed to atmospheric forcings. That corresponds to the pycnocline depth, calculated from three-dimensional density fields as the depth where the density increases by 0.03 kg · m −3 compared to a reference value taken at a 10-m depth [46]. If this criterion is not satisfied, the MLD is set to the maximum depth. Figure 2c shows the seasonality of the MLD from winter (left) to summer (right). DPE represents the amount of energy required to homogenize the water column and therefore characterizes the stratification intensity. It is calculated from [15,33]: where η is the sea surface elevation, H 0 + η is the total water depth, ρ is the water density, g is the acceleration of gravity equal to 9.81 m·s −2 andρ is the vertical mean water density.    Figure 2d shows the seasonal variations of DPE from winter (left), when the water column is weakly stratified and high values are mainly related to freshwater inputs associated to river plumes, to summer (right), when the stratification is strong except in coastal areas and in the English Channel area, mixed by strong tidal currents. The spatial distribution of the temperature and salinity tracers shows areas of small scale variations associated to frontal structures that delimit water masses with different hydrological properties [47][48][49]. These frontal structures can be detected using GRADSST and GRADSSS thanks to an arbitrary detection threshold (e.g., Figure 2e,f).
Dynamical variables related to physical mesoscale processes such as eddies [31][32][33] have been calculated: Root Mean Square (computed over one month) of low frequency currents filtered of the tidal component (RMS_FILT), Root Mean Square (computed over one month) of high frequency (tidal) currents (RMS_TIDE), Okubo-Weiss criterion (OW [50]), relative vorticity (RVORT), Eddy Kinetic Energy (EKE) and Mean Kinetic Energy (MKE). These dynamical variables have been calculated at a near-surface reference depth of 10 m rather than at the surface in order to avoid too high-frequency variations, due for example to the diurnal cycle and to accurately capture the turbulent activity and the baroclinic signature of the associated structures. RMS_FILT and RMS_TIDE were calculated using a Godin filtering method based on a weighted average of hourly-sampled data as in Simon [51] on each velocity component. An example of RMS_FILT and RMS_TIDE fields is given in Figure 3a,b, respectively. OW, RVORT, EKE and MKE were calculated from the zonal and meridional components of RMS_FILT. OW and RVORT were used to characterize vortex activities [50,52]. OW allows, thanks to the closed sign change contours, to detect the limits of vortex structures. The sign of RVORT is used to determine the nature of the vortex structure, cyclonic if it is positive and anticyclonic if it is negative. RVORT and OW are calculated as follows: where u and v are the zonal and meridian velocity components, respectively. An example of OW and RVORT fields is given in Figure 3c,d, respectively. MKE characterizes areas of strong current velocities and EKE is representative of the energy associated to turbulent activity. They are calculated following [31][32][33]: whereū andv are the time-mean velocity components. An example of EKE and MKE is given in Figure 3e,f, respectively.

Database
Each monthly mean physical field was calculated over an area of 471 × 720 points, corresponding to 339,120 points of which 97,678 points were land area and 151,162 points were beyond the 200 m isobath. The initial database was therefore made up of 339,120 points per time step × 12 physical variables × 60 months from January 2012 to December 2016, that being to 2.445 × 10 8 points.

Clustering Algorithm
In order to partition the abiotic component of the regional ecosystem, a preliminary statistical step consisted in a Principal Components Analysis (PCA) on scaled physical EOVs and subEOVs previously defined, for each month. PCA was used as a tool for dimensionality reduction, to reduce collinearity and to construct a composite unitless dataset, prior to the clustering analysis. Eigenvectors corresponding to 90% of the explained variance were selected to reduce the dimensionality of the data while minimizing the loss of information [53]. Partitioning was carried out using a hybrid two step clustering procedure using Hierarchical Agglomerative Clustering (HAC) and k-means. Hybrid approaches consist of the combination of different statistical methods in order to keep their individual advantages for a statistical protocol presenting certain constraints [54]. This two step procedure was chosen in order to (i) determine the number of clusters present at each month and (ii) find common clusters across months ( Figure S1). In the first step of the partitioning process, spatial clusters were assessed on a monthly basis by Hierarchical Agglomerative Clustering (HAC) performed on selected eigenvectors from PCA on oceanographic variables and using Euclidean distance and Ward's linkage [55]. HAC was chosen since its main advantage relies on its stability and because it does not require any a priori knowledge of number of groups. The optimal number of groups was determined using the Calinski-Harabasz criterion [56] within a range between 3 and 20, and over 500 iterations on the monthly dataset. As the HAC was carried out monthly, the clusters detected each month could not be reassigned to the following month. A second step has therefore been implemented to complete the regionalization process running k-means on the whole dataset (i.e., pluri-annual) to refine the monthly clusters labeling from HAC [57]. The optimal number of clusters was, for this "global clustering", again determined using the Calinski-Harabasz clustering evaluation criterion, within a range of k potential clusters between the minimum and maximum number of clusters detected by HAC in the monthly step. This hybrid procedure resulted in a time series of monthly categorical maps composed by 10 clusters over 60 months (2012 to 2016), named "patch class" as defined by [58] (Figures S2 and S3).
All analyses of this section were performed using the "stat" package of the software R (V10) and MATLAB. Ocean circulation model hindcasts, subEOVs and HAC computation were performed using high-performance computing resources of DATARMOR infrastructure (http://www.ifremer.fr/pcdm).

Patch Characterization
Generalized Boosted Regression Modeling (GBM) allows to build a nonlinear relationship between a response variable and explanatory variables. GBM maintains a monotonic relationship between the response (here, patch occurrence over 60 months) and each EOV and subEOV. GBM was used to characterize relative variable importance for each patch a posteriori. GBM models were chosen because Tree-based methods, such as RF or GBM, are more robust to collinearity [59]. Analyses were conducted on the 60-month average EOV climatologies for each patch. As patches were discriminated against the environmental variables, the object here was not to predict the patches with environmental variables, but to determine the most significant variables to characterize the groups. The GBM model calculates the relative influence of each variable in reducing the loss function [60]. We built models using the 75% training data to fit the GBM and the remainder 25% data were used to compute independent estimates of the loss function.
GBM models were constructed using the following parameters: The learning rate (lr) was equal to 0.01 to enhance the model's performance [60]. The maximum nodes per tree were 5 to allow the predictor's interactions. The bag fraction (the subsample fraction) was chosen in the medium range of 0.5, corresponding to the half of the training sample used to implement stochastic gradient descent to minimize overfitting. The relative variable importance, calculated by ranking the individual variables based on their relative importance, was computed using the default setting parameter as described in [61]. All analyses of this section were performed using the GBM package [60] of the software R (V10).
In addition, the categorization of variability for each EOV and subEOV was determined based on the 25th, 50th and 75th percentiles of the EOV and subEOVs distributions over the entire area and over the 60 months, and then compared to mean distributions of EOVs and subEOVs for each patch. Only the ranges of salinity variation were defined according to the following usual characteristics: PSU < 30, 30 < PSU < 34 , 34 < PSU < 35 and PSU > 35 characterizing estuarine waters, ROFI (Region Of Freshwater Influence) waters, shelf waters and full salinity oceanic waters, respectively. The ranges of variation are relative to the study area and allow a typology of the characteristics of the water masses to be established in order to synthesize the information.

Seascape Unit Definition
The concept of seascape implies to take into account the combination of ecological and strategic considerations including administrative delimitations [62]. Seascapes units were designed by computing the spatial median of each patch relative to the frequency of detection of the patch over the time series (60 months). A seascape was defined as the combination of the administrative delimitation of the French marine subregions of the MSFD and the patch mosaic. The seascapes were therefore composed of sub-hydrographical patches simply referred to as "patches". When a hydrological patch was split by administrative boundaries, it did not lose its physical properties and gained a new property linked to the marine sub-region to which it belongs. In addition, some contiguous patches may share boundaries with an existing part or overlap and their intersection as well as the mapping were conducted using the graphical information system QGIS.

Pattern Analysis
The time series of the number of clusters detected by the HAC shows a range of variation between 8 and 13 clusters (Figure 4). The evaluation of the optimal number of clusters for the global k-means on the whole dataset gives an optimal number of clusters equal to ten ( Figures S2 and S3). The ten patch classes are not systematically detected at each time step, indicating a spatio-temporal variability of the area partition illustrated in Figure 5, by the spatial occurrence of each group over the 60 months of the study ( Figure S4). Patch class 1 is distributed on the Bay of Biscay continental shelf, with a percentage of occurrence over 60 months between 10% and 80%, and a broad decreasing offshore-coastal gradient. The maximum stability is observed on the offshore part of the patch. Patch class 2 is a more coastal shelf patch, with a percentage of occurrence between 20% and 80%. A strong variability on the offshore part of the patch is observed, in contrast to Patch class 1. Patch class 3 is an oceanic group, located from the extreme northwest part of the study area to the entrance of the English Channel. This patch presents a strong disparity between the East and the West, reflecting a high variability in its eastern part over time and a very high stability in its western part with a percentage of occurrence over 60 months higher than 80%. Patch class 4 is close to the coast with a moderate percentage of occurrence, reflecting a relative variability in the area. Patch class 5, in the English Channel, shows a high stability in its central part, and variability in its periphery. Patch class 6 has a very small coastal extent and is highly stable along the period. Patch class 7, located just North of the whole area, has a spatial occurrence profile quite similar to Patch class 5, with high stability in the center of the cluster and a more variable periphery. Patch class 8 is located near the main rivers, with a high stability of detection occurrences. Patch class 9 represents the water masses along the Iberian coast, expressing non-recurring sub-patches along the French coast with a strong variability. Finally, Patch class 10 is located in the English Channel. It presents a fairly moderate percentage of occurrence over 60 months, probably related to its situation at the interface between Patch classes 3 and 5.
The analysis of patch classes area distributions according to the seasons reveals heterogeneity among groups ( Figure 6). Patches classes 1, 2, 5, 6, 7 and 8 do not show significant differences in their seasonal area distributions, with a Kruskal-Wallis test p-value greater than 0.05 (not shown). Conversely, patch classes 3, 4, 9 and 10 present a strong heterogeneity in their ranges of area distributions according to the season, with a Kruskal-Wallis test p-value smaller than 0.05 (not shown). Patch class 3 is characterized by a larger distribution area during winter with a gradual decrease in spring and fall, and a minimal spatial extent in summer. Patch class 4 shows a larger extent in winter and spring, then a drastic reduction in summer and autumn. Patch class 9 is large in autumn, and much reduced in other seasons. Patch class 10 is characterized by large areas in summer and autumn, against smaller ones in winter and spring.

Hydrographical Patches Characteristics: Relating Statistical Patches with Abiotic Process
The relative importance of each derived variable, based on GBM model results, is shown in Table 2. For information purposes, GBM evaluation and performance metrics are shown in Table S1. Ranges of categorization are defined in Table 3. As the GBM model was computed using climatology over the 60 months, the model enables an average view of the dominant variables, but does not adequately represent individual months, or a seasonal signal. Combined with the seasonal analysis based on percentile ranges for each EOV and shown in Table 4, results allow the characterization of each patch class which is described below.   Table 4. Seasonal hydrographical typology of the 10 patch classes. Links between typology and the 25th, 50th, 75th range of variability of EOVs distribution are shown in Table 3. Patch class 1 is mainly characterized by DPE and MLD. It shows a seasonal signal, with shelf waters marked by salinity fronts and stratification in summer and autumn, against oceanic well-mixed waters with a lower frontal activity in winter. Spring corresponds to an intermediate situation, with moderately mixed waters. This patch is also characterized by low and moderate dynamical metrics, for example EKE, RMS_FILT and RMS_TIDE.

Patch Number
Patch class 2 is mainly driven by MLD but GRADSST also plays a non-negligible role. It can therefore be summarized as a group with well-mixed shelf waters in winter and moderately stratified waters in other seasons with a strong frontal activity throughout the year, influenced by fresh water outflows.
Patch class 3 is driven by a hydrological regime conducted by DPE, GRADSST and SST. However, the seasonal analysis reveals a strong variability with dominant mixed waters in winter and spring, while water masses are well-stratified in summer and autumn. It is an oceanic group, with a full salinity (>35 PSU).
Patch class 4 is mainly characterized by SSS and GRADSSS. SSS variability shows a strong relation with ROFI waters, with SSS values between 30 and 34 PSU. It is also characterized by an intense frontal activity with strong haline and thermal fronts, including a moderate stratification of the water body.
Patch class 5 is mainly characterized by RMS_TIDE and RMS_FILT. The seasonal analysis also highlights well-mixed water mass patterns due to strong currents, turbulence and tidal flows, but no significant seasonal signal is noticeable. The frontal activity is quite moderate throughout the year. It is an oceanic group with a full salinity.
Patch class 6 is mainly characterized by strong dynamical structures, with high values of OW, RVORT and MKE as relative key variables. The seasonal analysis confirms the GBM results on dynamical structures, and highlights that this patch is also influenced by residual currents and tidal flows, inducing fully mixed waters, without seasonal signal. It is an oceanic patch with a full salinity. Patch class 7 is mainly driven by EKE corresponding to the mean flow coming from the Atlantic to the North Sea. The detailed analysis showed that it is also characterized by a seasonal signal with well mixed waters in winter and spring, against a relatively moderate stratification in summer and autumn, associated with strong tidal flows and residual currents with turbulent processes. This group is mainly a shelf water mass patch except in winter.
Patch class 8 is mainly characterized by the dominance of SSS and GRADSSS variables in the GBM model. This patch is an estuarine patch with reduced salinities and strong density gradients. The behavior of the patch does not show a seasonal signal.
Patch class 9 is mainly dominated by RMS_FILT. It corresponds to oceanic waters, with strong flows and frontal structures . It is also characterized by a seasonal signal, with a strong stratification associated with high frontal activity during summer and autumn against moderate stratification in winter and spring.
Patch class 10, unlike other patches, presents a certain heterogeneity with four relatively dominant variables: DPE, RMS_TIDE, GRADSSS, MLD and SSS, probably due to its position at the interface between Patch class 3 and Patch class 5 ( Figure 5). The seasonal analysis shows a strong variability with strong and moderate mixed layer depths in winter and spring, and a strong frontal activity in summer unlike for other seasons.

Seascape Units
The coupling between water masses partitioning and the administrative delimitations of the MSFD French marine subregions is shown in Figure 7a. Each seascape is characterized by the proportion of each patch in the seascape and the area of each patch (Table 5). In the English Channel-Greater North Sea region, seven patches are identified. Patch 2 EC_NS is a shelf water patch counting for 19% of the total marine sub-region. Patches 4 EC_NS and 8 EC_NS are small patches representing less that 1% of the marine sub-region. Patch 6 EC_NS is a small patch covering 2% of the marine sub-region. Patch 7 EC_NS is exclusively found in this marine sub-region and corresponds to the channel between the North Sea and the English Channel, characterized by a small spatial extent (less than 5%). Patch 5 EC_NS is the third contributing patch in terms of spatial extent (26%). Finally Patch 10 corresponds to a surrounding background (27% of the total marine sub-region area).
In the Celtic seas region, five patches are identified. The region is spatially composed of two main patches, the patch 3 CS is related to the patch 3 detected previously with a spatial extent representing 55% of the total area and patch 10 CS (36%). The patch 2 CS corresponds to the patch class 2 but its spatial extent is much lower than in the other two marine sub-regions (2%). It is the same for patches 5 and 6 in this sub-region. Table 5. Seascape composition metrics, patch proportion in the seascape and patch area.

Seascape Patch Patch Proportion Patch Area (km 2 )
Bay of Biscay S1 BB 0. 53  In the Bay of Biscay area, twelve patches are detected. Patch 1 BB is exclusively found in this marine sub-region and composes the main spatial part of the marine sub-region by occupying 53% of the area. Patch 2 BB is a shelf water mass group (19%). Patch 3 BB is the most southern part of patch 3 described previously counting for 17% of the marine sub-region area. In this subregion, Patch classes 4 and 8 are decomposed into 3 and 2 patches, respectively, according to the rivers influencing them. Thus, patches 4-1 BB, 4-2 BB and 4-3 BB are related to ROFI waters and patches 8-1 BB and 8-2 BB to estuaries. ROFI and estuarine patches represent a small fraction of the marine sub-region. Patches 5 BB, 6 BB and 10 BB are also present in this subregion, but their spatial extent is small (less than 2% of the total area). At least, patch 9 BB corresponds to a part of Iberian coastal waters into French waters, with a small proportion of the total area (1%).
Due to their spatial and temporal variability, some contiguous patches may share boundaries with an existing part or overlap. These overlapping areas are presented here as an interface between the patches (Figure 7b). The Celtic Sea shows a marked interface between patch 3 and 10, confirming the first results described previously. In the Eastern part of the English Channel and in the Bay of Biscay, main interfaces occur in transitional areas, between ROFI and shelf waters on one side and shelf water and oceanic waters on the other side. The Bay of Biscay is also characterized by an overlap between patch 3 BB and patch 1 BB.

Discussion
Seascape computation has required the implementation of added-value products from the OCO system: Calculation of EOVs and subEOVs and the classification of water masses patches. The input of the operational coastal oceanography has been crucial to this work. It indeed has enabled a procedure consistent with the international GOOS work [28], as well as the integration of data from CMEMS for some forcings and model validation [63] in order to build keystone steps for the implementation of MSFD-seascapes. The computation of several derived variables from OCO allowed us to synthesize the information provided by the state variables of the ocean circulation model. As such, they are a relevant tool to describe, for example, (i) stratification in the surface layers (mixed layer depth, intensity of stratification, etc.) from hydrological data, (ii) remarkable energy structures, such as frontal and eddy structures from current data (vorticity and Okubo-Weiss criterion). Our work is in line with previous literature such as [29,30,34], demonstrating that operational oceanography is a sufficiently mature and a major tool for the characterization of abiotic environmental conditions and providing high added-value products useful for marine management.
This work results in the development and the characterisation of a mosaic of dynamic hydrographic patches between 2012 and 2016 along the Atlantic coast of the French metropolitan EEZ. The hybrid combination of HAC and k-means allowed the construction of a time series of patch mosaics, resulting from a multivariate spatio-temporal partitioning. Each mosaic of the time series consists of a categorical map to define a picture of abiotic spatial units in the pelagic domain. From this time series, it was possible to define ten inter-annual main "patch classes" characterized by their spatio-temporal variability and hydrographical characteristics. Some are mainly dominated by hydrological processes such as the continental shelf Bay of Biscay patch class (number 1), the coastal shelf patch class (number 2), the Celtic Sea oceanic patch class (number 3), the river plume patch class (number 4) and the estuary patch class (number 8). Others are mainly dominated by hydrodynamical processes such as the tide-influenced patch class (number 5), the highly turbulent patch class (number 6), the high north flow patch class (number 7), and the coastal Iberian patch class (number 9). The transitional patch class (number 10) is the most variable one, probably reflecting its interface position. In addition, the analysis of the time series revealed differences in stability over time and space between patches, with stable patches (1,2,5,6,7,8) and less stable patches exhibiting significant seasonal variability (3,4,9 and 10).
In the context of the MSFD, specific sub-patches have been detected in each marine sub-region by combining the initial main patch classes with the MSFD regulatory delimitations. This combination corresponds to the definition of a "seascape" composed by sub-hydrographical patches and human consideration. Comparing our results to the literature, it appears that seascapes described and characterized here are consistent with previous works based on ocean partitioning, or works based on specific areas focusing on the processes involved. In the Eastern English Channel-North Sea French marine sub-region, our seascape is consistent with the one of [64], although the scope of the partitioning process and the choice of the derived variables were different. Their number of patches, within a range between seven and nine following the season, is quite similar to our results of seven detected water masses. Their classification showed some stable water masses comparable to the patches 8 EC_NS, 4 EC_NS, 5 EC_NS, 2 EC_NS and 10 EC_NS in the eastern part of the English Channel. They also detected patterns in the same geographical area as patch 7 EC_NS, but with a variability difficult to interpret. Patch 8 EC_NS corresponds to the Seine estuary and patch 4 EC_NS to its river plume. Patch 5 EC_NS is influenced by strong tidal currents in the Alderney Race area with the particularity to be the location of the strongest tidal currents in Europe. The shelf water patch 2 EC_NS corresponds to the coastal flow along the French coast generated by fluvial supplies as described by [65]. Patch 10 EC_NS is an oceanic transitional patch water. Patch 6 EC_NS, characterized by strong EKE and tidal currents, corresponds to many tide-induced gyres well described by [66].
The Celtic sea seascape is composed of five patches. Patch 3 CC in the western part is the biggest one, with a high seasonal dynamics. It reaches its maximum extension in winter. During spring, patch 10 CC enters the eastern part of patch 3. Since patch classes 10 and 3 are neighbours, and their spatial extent dynamics are inverted, it can be assumed that this corresponds to back and forth motion of the water mass from east to west depending on the season. This interface between patches 3 CC and 10 CC corresponds to the tide-induced Ushant frontal region. The Ushant front is set up at the end of spring. It is characterized by a stratified water mass offshore (patch 3 CC) and mixed waters near the coast (patch 10 CC) [67]. The presence of patches 5 CC and 6 CC, mainly dominated by tides and high turbulence is due to strong currents in the Fromveur Strait [68].
In the Bay of Biscay seascape, twelve patches are detected. Our results are consistent with those found by [15], but only for spring. The patch 3 BB corresponds to their group 8 "Northwestern shelf", characterized by oceanic waters and thermal stratification in the springtime. The patch 1-BB is redundant with their patch 7 "Open shelf" along the continental slope, with oceanic waters and moderated stratification. Patch 2-BB corresponds to their patch 6 "Central shelf" characterized by shelf waters (34 PSU < SSS < 35 PSU) with the influence of outflow fresh water, moderate stratification, and low mixing flow during spring. For the other more coastal groups in the Bay of Biscay, the patterns are quite consistent, even if the number of groups differs. We find typical ROFI, patches 4-1 BB, 4-2 BB, 4-3 BB corresponding to groups 5, 4 and 3, respectively in [15], and estuarine water mass patterns 8-1 BB, 8-2 BB corresponding to group 2 in [15]. Patches 4-1BB and 8-1BB are related to the Loire river plume and estuary. Patches 4-2 BB and 8-2 BB correspond to the Gironde river plume and estuary. Patch 4-3 BB is related to the Adour river plume. Some patches present a strong seasonality in terms of spatial extent. For example, ROFI patches 4-1 and 4-2 BB are larger in winter and spring than in summer and autumn, which corresponds to the plumes of the Loire and the Gironde rivers analyzed in [69]. However, a difference can be noted concerning the period of maximum plume extension [69], using satellite turbidity data, showed that the maximum plume extension occurs in winter, whereas in our case it occurs in spring. It may be hypothesized that this difference could be due to the fact that [69] worked only on turbidity and that our study focuses on both hydrological and dynamical metrics, and also because the periods of study are different. It does suggest, however, the opportunity to integrate satellite turbidity data in a future work. Patches 5 BB, 6 BB and 10 BB are also present in this seascape, but are negligible compared to the global reach of the patch class. Patch class 9 on the Iberian coast corresponds to water masses along the continental slope bordering the Spanish coast. It is characterized by a strong seasonal signal with a maximum geographical extension in autumn, corresponding to water masses forming a tongue rising northwards at 2.5°W described by [70].
Our methodology has shown much strength to define seascape in the context of marine policy, but such partitioning methods have also some limitations, in particular on the detection of the patch class 9. The analysis of the global k-means ( Figure S3) detects patch class 9 along the Spanish coast and in the Dover strait at the extreme north-east of the study area. However, it is unlikely to have system similarities between the Spanish coast and the Dover strait. By choosing to limit the study area to the 200m isobath, we included this very narrow zone along the Spanish coast. It appears, however, that this zone should not have been taken into account in the analyses, since it potentially generates a statistical bias, generating a statistical group that does not correspond to the oceanographic reality and is thus difficult to interpret. Furthermore, as shown in the ocean model description, low-frequency boundary conditions are provided by a large-scale ocean model. To prevent any inconsistency between boundary forcing and the interior computation domain, a sponge layer is located in the vicinity of the domain boundaries, where numerical solution is relaxed to boundary forcing. This numerical procedure is nevertheless likely to create artifacts near boundaries, and may lead to wrong patch detections. This may explain the pattern along the northern boundary in the oriental part of the English Channel and in the Celtic Sea. Several methodological improvements to the present work can thus be considered. The first one would be better defined boundaries of the study area, as previously stated. Another one would be to perform partitioning over smaller spatial extents, such as patches defined in this study, but this would raise questions about the balance between the need of a regional assessment and the relevant grain (spatial resolution) of the ocean circulation model. This indeed implies the use of high resolution coastal circulation models, and the need for observation networks adapted to their validation. On the French coasts, it is hoped that ILICO infrastructure [71] will enable significant progress on these issues, linked to operational oceanography in the coastal-nearshore domain.
The partitioning proposed here, points out the importance, for marine policies, to take into account not only hydrological variables related to stratification or temperature, but also dynamical ones related to oceanic currents, in order to better assess seascapes in the pelagic domain. Only based on abiotic characteristics, it allows to represent a specific level of the systems hierarchy, immediately above the ecosystem [72], as a counterpart of the benthic broad habitats described in the EUNIS classification. This strata is thus the abiotic part of the "ecomosaic" in the pelagic realm [72]. Furthermore, the development of a time series of seascapes categorial maps enables the use of usual landscape metrics to characterize seascapes in terms of composition and spatial configuration [58]. Ecologically speaking, this study provides a better understanding of the complexity of coastal-shelf areas in the french EEZ by highlighting the spatio-temporal variability of hydrographical patches. Indeed, seascapes variability induces moving boundaries between patches, i.e., interfaces that may be a proxy for marine ecotones, which are relevant places for ecological phenomenons such as nutrient injection, aggregation of matter, connectivity or trophic level elevation [3,14,24]. Thus, the assessment of a patch mosaic within seascape, based on hydrology and hydrography, could (i) help at defining areas of interest for human activities (e.g., strong current areas for marine renewable energy) and; (ii) allow the detection of, on the one hand, stable areas potentially hosting ecologically mature biocenosis, and, on the other hand, singular and dynamic patches, possibly associated with singular biological communities, which may require specific management measures. We therefore propose that detailed biological informations have to be integrated a posteriori, in a Russian doll's logic, such as in the EUNIS classification, in order to assess finer strata of the hierarchical pelagic system. This would allow to assess the mosaic of biocenosis (small units) hosting in patches defined here (medium units) within the seascape (large units). This will be investigated in future works, in particular by integrating a posteriori satellite ocean color data or outputs from biogeochemical models.

Conclusions
The objective of this study was to compute and analyse a mosaic of abiotic patches in the coastal-shelf system of the french EEZ. The development process was carried out in a pluridisciplinary framework, integrating ocean modelling tools, statistical tools to compute partitioning and water mass characterization, in the context of marine policies. Along this framework, a precise set of specifications appeared to be relevant to implement: (i) It has to be based on hydrological and dynamical ocean variables calculated from outputs of a three-dimensional ocean circulation model such as HYCOM, (ii) partitioning the continental shelf area should take into account its dynamics over space and time and (iii) a hybrid partitioning process has to be computed. We therefore propose mosaics (seascapes), composed of standardised and understandable spatial units (patches) based on oceanic characteristics, and taking into account the combination of ecological and strategic considerations [62], such as marine spatial planning, including cross-border issues. Seascapes are a tool providing common synthetic environmental information which can enhance knowledge on the structure and architecture of ocean ecosystems, can better target marine monitoring on poorly known areas in order to better assess the state of water masses, and can give a setting for the development of operational indicators in the context of the MSFD.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/2077-1312/8/8/ 585/s1. Figure S1: Pipeline of the methodology of the hybrid two-step clustering, Figure S2: Evaluation of the optimal number of clusters for the global k-means using the Calinski-Harabasz clustering evaluation criterion, within a range of k potential clusters between the minimum (8) and maximum (13) number of clusters detected by HAC, Figure S3: Global k-means results on the whole data set, with an initialization with 10 clusters, Figure S4: Distribution of the number of clusters after the second-step procedure, Table S1: Performance and evaluation metrics for the GBM model.