A Methodology and Tool for Mapping the Risk of Cumulative Effects on Benthic Habitats

The implementation of the European integrated marine policy poses many scientific challenges. Among them, the knowledge and understanding of interactions between anthropogenic pressures and ecological components is an important issue, particularly to help define Good Environmental Status, environmental targets and monitoring programs of the Marine Strategy Framework Directive (2008, MSFD). Assessment of cumulative effects of different pressures is a particularly complex issue requiring modeling tools and methods, as well as accurate data sets on human activities, anthropogenic pressures and ecological components. The results of these assessments are also uncertain and highly dependent on the calculation methods and assumptions, as well as on the data sets used. Within this context, we developed a technical and methodological approach to map the risk of cumulative effects of different pressures on benthic habitats. These developments were initiated as part of the implementation of the MSFD in France to contribute to the diagnosis of the marine environment. Here we provide a demonstrator to illustrate the feasibility for mapping the risk of cumulative effects of different pressures on benthic habitats, as well as the confidence index and the variability associated with this analysis. The method is based on a spatial analysis using a mapping of benthic habitats and their sensitivity to pressures, as well as the distribution and intensity of human activities and associated pressures. We collected and prepared relatively accurate and consistent data sets to describe human activities and benthic habitats. Data sets are embedded into a grid that facilitates the management and analysis of the data and exploitation of the results. The demonstrator consists of a relational database using the Spatial Query Language (SQL) language as well as data analysis scripts using the R language. The first demonstrator operations validated the main methodological and technical choices and helped to identify future developments needed to facilitate the appropriation and integration of these approaches in the implementation of public policies for the management of the marine environment.

The implementation of the European integrated marine policy poses many scientific challenges. Among them, the knowledge and understanding of interactions between anthropogenic pressures and ecological components is an important issue, particularly to help define Good Environmental Status, environmental targets and monitoring programs of the Marine Strategy Framework Directive (2008, MSFD). Assessment of cumulative effects of different pressures is a particularly complex issue requiring modeling tools and methods, as well as accurate data sets on human activities, anthropogenic pressures and ecological components. The results of these assessments are also uncertain and highly dependent on the calculation methods and assumptions, as well as on the data sets used. Within this context, we developed a technical and methodological approach to map the risk of cumulative effects of different pressures on benthic habitats. These developments were initiated as part of the implementation of the MSFD in France to contribute to the diagnosis of the marine environment. Here we provide a demonstrator to illustrate the feasibility for mapping the risk of cumulative effects of different pressures on benthic habitats, as well as the confidence index and the variability associated with this analysis. The method is based on a spatial analysis using a mapping of benthic habitats and their sensitivity to pressures, as well as the distribution and intensity of human activities and associated pressures. We collected and prepared relatively accurate and consistent data sets to describe human activities and benthic habitats. Data sets are embedded into a grid that facilitates the management and analysis of the data and exploitation of the results. The demonstrator consists of a relational database using the Spatial Query Language (SQL) language as well as data analysis scripts using the R language. The first demonstrator operations validated the main methodological and technical choices and helped to identify future developments needed to facilitate the appropriation and integration of these approaches in the implementation of public policies for the management of the marine environment.

INTRODUCTION
Worldwide marine and coastal ecosystems are hosting more and more human activities using space, mineral and living resources. This leads to an increase in anthropogenic pressure, causing a significant decrease in biodiversity, habitat loss and significant changes in ecological functions (Korpinen and Andersen, 2016;Willsteed et al., 2017). The understanding and assessments of the cumulative effects of anthropogenic pressures are considered as major issues to inform strategic planning and marine ecosystem conservation and management (Foley et al., 2017;Stelzenmüller et al., 2018).
Cumulative effects assessment (CEA) has been defined as 'a systematic procedure for identifying and evaluating the significance of effects from multiple sources/activities and for providing an estimate of the overall expected impact in order to inform management measures' . Cumulative effects assessment can contribute to the diagnosis of the marine environment and provide strategic information to support decision-makers (Korpinen et al., 2013;Goodsir et al., 2015;Depellegrin et al., 2017). Over the past decade, numerous studies have laid the methodological bases of CEA Stelzenmüller et al., 2010;Micheli et al., 2013). This field of investigation is particularly dynamic and is still very exploratory. The methodological and semantic framework, as well as shared principles, conceptual models and tools, are still under development and raise complex scientific questions (Korpinen and Andersen, 2016;Stelzenmüller et al., 2018).
As mentioned in Korpinen and Andersen (2016) about half of the cumulative pressures and impacts assessment studies follow the same method or a similar method to that presented by Halpern et al. (2008). These assessments use mainly three components: the mapping of pressure intensities, the mapping of ecological components to be studied and a sensitivity index (weighting factor) for each ecosystem -pressure couple. The calculation method is most often additive, that is to say, the effects of each pressure on each component of the ecosystems are added. Significant methodological challenges and uncertainties are present at each stage of this cumulative effects modeling process (Halpern and Fujita, 2013;Korpinen and Andersen, 2016).
The ecosystem data sets have uncertainties about the presence of habitats, species and communities, their natural dynamics and ecological status. The human activity data sets have uncertainties about the distribution and intensity of uses in time and space. Moreover, the relationships between activities and pressures and between pressures and ecological components also include uncertainties and hypotheses. Finally, the methods of calculating the cumulative effects of different pressures also rest on many assumptions and uncertainties such as the relative contribution of activities to each pressure, the combined action of pressures, the response of ecosystem components, and especially the variability of their sensitivity according to their ecological status and local conditions. These issues in addition to the different stages of data preparation such as data integration in a grid, the use of specific typologies, data transformation and normalization, have a great influence on the final results. Uncertainty analysis is becoming an increasingly important aspect of modeling processes and is now considered as a mandatory step in cumulative effects assessment studies (Stock and Micheli, 2016;Gissi et al., 2017).
The Marine Strategy Framework Directive (2008/56/EC of the European Parliament and the Council, MSFD) aims to achieve Good Environmental Status (GES) by 2020 (European Commission, 2008) following an ecosystem-based management. The development of methods and tools has been identified as one of the major challenges to contribute to implementation of the MSFD and help decision-makers Borja et al., 2017). As described in the Directive, the ecological assessment must 'cover the main cumulative and synergistic effects.' However, for the past 10 years and the beginning of the MSFD implementation, anthropogenic pressures have been evaluated one by one and there are no validated tools and methods for evaluating cumulative effects.
This article presents a concrete application of the concepts presented above to provide methods and demonstrator tool able to map the Risks of Cumulative Effects (RCE) of different pressures on benthic habitats. These developments are consistent with recently developed concepts to support an integrated approach for assessing the status of benthic habitats under the MSFD (Elliott et al., 2018). In this context, the developments presented here are intended to contribute to the diagnosis of the marine environment, usable for the implementation of the MSFD, the Maritime Spatial Planning Directive and for marine protected areas management.
The proposed RCE assessment is based on Halpern et al. (2008) method. We enriched this approach with a human activities-pressures relationship matrix based on La Rivière et al. (2016) recent work that we complemented with the consultation of scientific and institutional experts. Moreover, we adapted and enriched this method with a new matrix of benthic habitats' sensitivity to anthropogenic pressures (La Rivière et al., 2016). The different data describing human activities, pressures and benthic habitats as well as relationship matrices (activity-pressure and habitat-pressure) were associated with a spatially explicit confidence index. Finally, we developed two approaches to map the confidence and the variability of the RCE depending on different calculation assumptions . The intermediate data produced on human activities, pressures and benthic habitats constitute a homogeneous set of data, contributing to a better interpretation of the RCE.
Results are presented and provide an initial assessment of the capabilities of the demonstrator tool for locating, explaining and prioritizing threats and data gaps in order to inform management and monitoring priorities based on the intensity and distribution of RCE. Gaps, limitations and scientific, methodological and technical challenges are discussed with a view to fostering the acceptability, appropriation and use of these approaches and methods.

MATERIALS AND EQUIPMENT
The demonstrator tool essentially uses free and open source software and widely used languages -i.e., R version 3.5.1 (R Core Team, 2018), Structured Query Language (SQL), Geographic Information System (GIS).
The software used for the preparation, storing and processing of the data is presented in Figure 1. Data formatting, visualization and mapping were performed in QGIS 3.2.1 (QGIS Development Team, 2018). Data storage and management were performed in PostgreSQL 10 database (The PostgreSQL Global Development Group, 2018) using pgAdmin 4 software (pgAdmin Development team, 2018) and PostGIS 2.4 extension (PostGIS Project Contributors, 2018). The RCE analysis was performed with R programming language 3.5.1 (R Core Team, 2018) using RStudio 1.1.4 software (RStudio Team, 2016). The configuration of the access to the database, the results export management and the backup options of the analysis were carried out in an MS Excel file that is read by the R script. This configuration file makes the developed tool accessible to users not familiar with R and allows users to have an overview of the information to fill out and check before launching an analysis.

METHODS
Since the 1990s international organizations have promoted the DPSIR framework (Drivers-Pressures-State-Impact-Response) as an integrative systemic framework to assemble bio-physical knowledge in support to the design and implementation of environmental policies (OECD, 1993;European Commission, 1999). This framework is widely used, particularly in the area of aquatic ecosystem management (Borja et al., 2006;Lewison et al., 2016). It is explicitly referred to in the MSFD regulation and commonly used to report on the regular updates of the initial assessment of the state of marine waters, pressures and impacts, and good environmental status objectives. It has also been criticized for failing to integrate properly the socioeconomic dimensions of environmental problems, leading to the proposal of enhanced frameworks such as DAPSIWR that considers key elements regarding "actions and activities" (A) controlling pressures and the impacts on "human wellbeing" (W) . Characterizing and assessing "impacts" in operational terms remains a challenge. On purpose, and to avoid misinterpretation, the CEA uses the term "effect" rather than impacts under a risk assessment perspective.
The term "Pressure" used here as a synonym of "anthropogenic pressure" refers to the mechanism through which a human activity can have an effect on a habitat. Pressures can be physical, chemical or biological. A same pressure can be caused by a number of different activities (La Rivière et al., 2016).
The term "Risk" used as a synonym of "risk of impact" or "vulnerability" is the combination of the likelihood of a habitat being exposed to a pressure (threat) and impacted by this pressure, depending on the ecological sensitivity to that pressure. "Risk" puts the scope of results into perspective, as they are not a quantitative expression of the measurable effects on biological communities.
The term "Effects" refers to the consequence of a pressure on a habitat where a change in its biotic and/or abiotic characteristics occurs (La Rivière et al., 2016).

Methodology Overview
Statistical and spatial data sets on human activities, pressures, and benthic habitats were synthesized, structured, and mapped into a regular square grid of 1 min of degree in latitude per 1 min of degree in longitude resolution. In the database, each grid cell was associated with a unique identifier, which made it possible to link these datasets. The RCE index was calculated in each grid cell taking into account the intensity of the human activities present, the pressures generated by these human activities, the surface area of the benthic habitats present as well as their sensitivity to pressures. The main steps of the analysis are presented in Figure 2. The descriptive data of human activities, as well as the calculated data describing the intensity of physical pressures, were stored and organized with classifications compatible with the indicative list of human activities and anthropogenic pressures contained in the amended Commission Directive (EU) 2017/845. The intensities of the physical pressures were estimated on the basis of intensities of the underlying human activities. As there were no data sets available where the actual pressure had been measured, the different human activity intensities were considered as proxies to estimate the intensities of the physical pressures. The descriptive data of the benthic habitats were stored and organized with the European Union Nature Information System (EUNIS) classification (ETC/BD-EEA, 2012).
The calculation of the RCE involved two relationship matrices. First, the activity-physical pressure matrix identified the link between human activities and anthropogenic pressures. Secondly, the benthic habitat-physical pressure sensitivity matrix gave the level of ecological sensitivity of benthic habitats to anthropogenic physical pressures. These two matrices made it possible to relate all human activities, pressures and benthic habitats in a consistent way. They formed a linkage framework which was inspired by previous studies Knights et al., 2015).
The human activity-physical pressure matrix identified with a binary link (0 or 1) the different human activities that could contribute to generate each physical pressure. We built this matrix using the typology of human activities adopted under the MSFD and the physical pressures typology from La Rivière et al. (2016). The human activities typology was enriched with more precise subcategories allowing them to be linked to real data describing the activities. We informed the relationships between activities and pressures in a 1-day workshop (December 2016) with the different scientific and institutional teams involved in the MSFD implementation. The 35 Workshop attendees were divided into 4 groups who each worked independently on one of the four parts of the matrix. The four parts of the matrix completed during the workshop are presented in the Supplementary Table 1. Finally, a working group composed of 4 thematic experts of the French Office for Biodiversity analyzed, completed and validated the matrix obtained during the workshop. The source of each activity-pressure relationship was documented and associated with a confidence index that expressed the level and the quality of expertise mobilized for the evaluation as described in the Supplementary Table 2. An extract of the activity-pressure relationship matrix is presented in Table 1.
The benthic habitats-physical pressures sensitivity matrix was based on Marine Evidence-based Sensitivity Assessment (MarESA) for English Channel and Bay of Biscay habitats (Tillin et al., 2010;Tillin and Tyler-Walters, 2014;Tyler-Walters et al., 2018) and on French benthic habitats' sensitivity assessment (La Rivière et al., 2018) for the Mediterranean habitats. Within the framework of the RCE methodological developments, the MarESA sensitivity assessment is used until the French benthic habitats' sensitivity assessment is completed and accessible for English Channel and Bay of Biscay habitats. The French benthic habitats' sensitivity assessment is based on the methodology developed previously under the MarESA project, which allowed us to easily combine the sensitivity assessments of the two projects. The sensitivity is the combination of a habitat's capacity to tolerate a pressure (resistance) and the time needed to recover after an impact (La Rivière et al., 2016). The resistance is defined as the ability of a habitat to tolerate a pressure without significantly changing its biotic or abiotic characteristics (La Rivière et al., 2016). The resilience is defined as the time a habitat needs to recover from the effect of a pressure, once that pressure has been alleviated (La Rivière et al., 2016). For each benthic habitat and each anthropogenic pressure, a sensitivity score (derived from the assessment of a resistance and a resilience score) and its confidence index were evaluated by a detailed literature review, compilation of evidence, complemented when necessary by expert judgment (La Rivière et al., 2016). Sensitivity assessments were made against single pressures, therefore cumulative pressures could not be considered, and resilience was only considered if the pressure had been alleviated or reduced to a magnitude that no longer caused an impact. It was therefore assumed that the sensitivity score provides theoretical information that cannot truly reflect the field and real conditions in which pressures and habitats interact. Table 2 shows an extract of the sensitivity matrix.

RCE Index Calculation
We adapted the cumulative effects model proposed by Halpern et al. (2008). In relation to this work, the main adaptation concerned the use of more precise activity and pressure typologies and the use of a matrix of relationship between activities and pressures. In our analysis, the activities were not integrated into the RCE calculation directly by considering them as pressures. Each activity could contribute to several pressures. The calculation of the RCE index for each grid cell was divided into three main stages: (i) computing the intensity of each anthropogenic pressure, (ii) computing the sensitivity of the grid cell to each anthropogenic pressure, (iii) the calculation of the risk of cumulative effects. Several steps of the calculation of the RCE index required a normalization of some variables x (e.g., activity, pressure, or sensitivity) into a new variable x norm following the equation below: Where x max and x min are, respectively, the maximum and minimum value of x in the study area. The aim of this normalization was to allow comparison of various variables expressed in different units by rescaling from 0 to 1, but also to avoid over-representation of extreme values and correct frequency distribution bias (Micheli et al., 2013;Depellegrin et al., 2017).
(i) Computation of the intensity Pj of each anthropogenic pressure j where A i norm is the normalized intensity of the activity i and γ Ai/Pj is the binary link of the activity -pressure relational matrix (0 indicates that the activity i does not generate the pressure j and 1 indicates that the activity i generates the pressure j).
(ii) Computation of the sensitivity SensPj to each anthropogenic pressure j Where S k is the surface of the benthic habitat k in the grid cell, S is the total surface of the grid cell and µ j,k is the sensitivity of the habitat k to the pressure j according to the benthic habitatsphysical pressures sensitivity matrix.
(iii) Computation of the risk of cumulative effect index (RCE). Prior to this stage, for each pressure j, the intensity Pj and the sensitivity of the grid cell SensPj were normalized with the normalization function described above (x norm ). Note that this imply that a normalized pressure intensity of 0 and 1 correspond, respectively, to the minimum and maximum intensity of that pressure in the study area and for the period considered.

Confidence Index Calculation
A first approach to spatially represent the uncertainties of the model was to develop a single confidence index that qualifies the reliability of the final RCE in each cell (Stock and Micheli, 2016). The calculation of this index first required assigning four intermediate confidence indices in each grid cell to: -the mapping of each benthic habitat k (CI_H k ) -the intensity of each human activity i (CI_A i ) -each relationship of the activity-pressure relational matrix (CI_γ Ai/Pj ) -the sensitivity scores of each link of the habitat-pressure sensitivity matrix (CI_µ j,k ).
The different methods for assigning these confidence indices are presented in Supplementary Tables 2-5. Based on these 4 intermediate indices, a global confidence index, reflecting the reliability of the RCE score, was computed for each cell. Similarly, to the RCE computation described above (2.2), there are three main stages: (i) computation of a confidence index for each pressure Pj (CI_Pj), (ii) computation of a confidence index for the sensitivity of each grid cell to each pressure Pj (CI_SensibPj) and (iii) the computation of the final confidence index (CI_RCE). The final index ranges between 0 and 1, 1 reflecting a maximum confidence.

Confidence Index for Each Pressure P j (Ci_Pj)
The confidence index associated with the pressure P j (CI_P j ) was calculated in 3 steps. First, CIL_AP_P j reflects the confidence in the link between all human activities and the pressure Pj. Two cases were distinguished: if one or more activities are present in a given grid cell, CIL_AP_Pj is the mean of CI_γ Ai/Pj weighted by the intensity of each activity, and divided by 5 to obtain an index ranging from 0 to 1. If no activity is present in the grid cell, the link between human activities and the pressure Pj will not have any effect on the RCE: we thus assigned a maximum value to CIL_AP_Pj.
Secondly, CIQ_AP_P j represents the confidence in the activity data sets associated with the pressure P j . Here as well, two cases are distinguished. If at least one activity is linked to the pressure Pj, CIQ_AP_Pj is the mean of the CI_Ai of the activities associated with the pressure Pj. If no activity is associated with the pressure Pj, the intensity of each activity Ai will have no impact on the RCE: we thus assigned a maximum value to CIQ_AP_Pj.
The confidence index in each pressure Pj (CI_Pj) is the product of these two intermediate indices and ranges from 0 to 1 as well:

Confidence Index for the Sensitivity to Each Pressure Pj (CI_SensibPj)
Computing a confidence index for the sensitivity of each grid cell to each pressure Pj first required us to assign a confidence index to the sensitivity scores of the k habitats of the grid cell to the pressure Pj. This index was computed as the mean of the CI_µ j,k weighted by the surface S k of each habitat k (and divided by 5, CI_µ j,k maximum value, to standardize the index between 0 and 1).
Secondly, we assigned to each grid cell an index that reflected the confidence in the mapping of benthic habitats that were sensitive to the pressure Pj (CIQ_HP_Pj). When at least one habitat was sensitive to the pressure Pj, CIQ_HP_Pj was computed as the mean of the CI_Hk, weighted by the sensitivity of each habitat to Pj (µ jk ) (and divided by 5, CI_Hk maximum value, to standardize the index between 0 and 1). In the case when no habitat of the grid cell was sensitive to Pj, the habitat surface had no impact on the RCE calculation, so we assigned a maximum score for CIQ_HP_Pj.
CI_ SensibPj was calculated as the product of these two indices:

Calculation of the Risk of Cumulative Effects' Confidence Index
The confidence index in the pressure Pj (CI_P j ) and the confidence index in the habitat's sensitivity to the pressure Pj (CI_ SensibPj ) allowed us to calculate a confidence index for the risk of effect of the pressure Pj in the grid cell (CI_RE_P j ). For this step, it is necessary to distinguish 3 distinct cases: -If the pressure Pj is present in the grid cell (Pj > 0) but the grid cell is not sensitive to Pj (µ jk = 0), then the risk of effect is null and its value depends solely on the sensitivity of the grid cell to the pressure (the intensity of the pressure has no impact on the result). In this case, it was considered that CI_RE_P j only depends on CI_ SensibPj . -If Pj is null in the grid cell (P j = 0) but the sensitivity of the grid cell to the pressure j is different from zero (µ jk = 0), then the risk of effect is null and its value depends solely on Pj (the sensitivity of the grid cell to the pressure will have no impact on the risk of effect). In this case, it was considered that CI_RE_P j depends only on CI_P j . -In the other two possible scenarios (both the sensitivity and the pressure are null or non-null), then the value of the risk of effect will depend on both µ jk and P j .
The confidence index in the RCE for a given grid cell (CI_RCE) was computed as the mean of the confidence indices of the risk of effect associated with each pressure Pj (CI_RE_Pj) weighted by the risk of effect associated with each pressure Pj (RE_Pj). In other words, the CI_RCE would mostly depend on the confidence indices in the pressures that contribute the most to the RCE. In the case of a null RCE ( RE_P j =0), we considered that all the Pj contributed equally to this result so the CI_RCE was simply computed as the mean of the CI_RE_Pj.

Uncertainty Analysis
To quantify the uncertainties associated with our model, we investigated the effect of the variation of several sources of uncertainty on the RCE scores. Based on a literature review (Halpern and Fujita, 2013;Korpinen and Andersen, 2016;Stock and Micheli, 2016;Gissi et al., 2017), we identified seven sources of variation -hereinafter referred to as "factors." These factors, presented in Table 3, can be divided in two categories: -Model assumptions (factors X1, X2, X6 and X7) for which the level of uncertainty, as defined by Walker et al. (2003), corresponds to scenario-uncertainty. This means that we know the range of possible assumptions, but current knowledge is not yet sufficient to formulate the probability of any one particular assumption being true. -Variation of several parameters of our model depending on the value of the confidence index associated with these parameters (factors X3 to X5). As proposed in recent studies (Stock and Micheli, 2016;Gissi et al., 2017), we performed Monte Carlo (MC) simulations to explore the range of possible RCE scores taking into account these factors.

X1: Multiple Pressure Effect Model
In most cumulative impact studies, human pressures are assumed to simply add up while it has been demonstrated that antagonist (i.e., the effect is lower than the sum of the pressures) and synergistic (i.e., the effect is higher than the sum) effects also occur in nature (Crain et al., 2008). We investigated the multiple pressure effects by applying randomly one of these three models in each MC simulation. Figure 1 shows the theoretical result of these three models.

X2: Activity Intensity to Pressure Intensity Model
The second factor concerned the method of calculating the intensity of physical pressures from the intensity of activities. The linear model which was applied in our default RCE calculation setting assumed that the intensity of the pressure generated by an activity (normalized between 0 and 1) equals the intensity of this activity (normalized between 0 and 1): The "pessimistic" approach considered that the intensity of the pressure increases very rapidly and then reaches a plateau with the equation: The "optimistic" approach considered that the intensity of the pressure increases slowly and is only significant at high values of A i : where P j A i is the intensity of the pressure P j generated by the activity intensity A i .
The last approach considered that the relationship between a pressure intensity and the activity that generates it is represented by a logistic curve: The Supplementary Figure 2 illustrates the four models used in the MC simulation. The type of activity-pressure response was specific to each activity-pressure pair. Therefore, during each MC simulation, each activity-pressure pair was randomly associated with one of the four options (linear, pessimistic, optimistic or logistic). In this way, for a given simulation, each activity-pressure pair was associated with the same function in all the grid cells.
X3: Value of the Sensitivity Score µ jk Instead of having a fixed sensitivity score, we associated the sensitivity scores µ jk with a statistical distribution and generated in each MC simulation a random value µ jk taken from this distribution. The aim was to take into account the confidence in each sensitivity score, so that the sensitivity scores with a low confidence index were more variable than the ones with a high confidence index. Because the sensitivity score µ jk ranges strictly between 0 and 5, a beta distribution was chosen to model each µ jk in order to maintain the generated sensitivity scores within finite borders. Based on the recommendations of Ferrari and Cribari-Neto (2004), we generated random values µ jk with the following model: Where ϕ jk is a precision parameter, which is proportional to the confidence index in the sensitivity score: The Supplementary Figure 3 shows the results obtained with 20,000 runs of this calculation process for some examples of µ jk and CI_µ j,k pairs.

X4: Value of the Activity -Pressure Relationship Index γ Ai/Pj
Following the same principle as X3, this factor consisted in varying the values contained in the activity-pressure relational matrix (γ Ai/Pj ) depending on the value of the confidence index (CI_γ Ai/Pj ). Since the activity-pressure relationship could only take two values (0 or 1) we used a simple application of the binomial distribution. The confidence index (CI_γ Ai/Pj ) values were between 0 and 5. Let γ Ai /Pj be the new γ Ai/Pj relation, the probability for this link to have the same value as the original link is then: At each simulation, a new activity-pressure relational matrix for the whole grid was randomly generated from the original matrix.

X5: Value of the Activity Intensity A i
Following the same principle as X3, this factor consisted in varying the intensity of activities (A i ) according to the confidence index (CI_A i ) describing the quality of the activity data set. For each simulation, we applied for each cell the same algorithm as factor X3.

X6: Spreading of Pressure Effect From Point Source
As proposed in Stock and Micheli (2016) we investigated the spreading of pressure intensity around the point source. The distance at which the activity exerted a given pressure was specific to an activity-pressure pair. As we did not have sufficient information to assign a fixed effect distance to each activitypressure pair, we generated at each simulation an activitypressure matrix containing a random number between 0 and 3. This number corresponds to the grid cell effect distance over which each activity exerts each pressure for this simulation as described in the Supplementary Figure 4. The function developed performed for each activity and in each grid cell the following calculation: If the activity was present on the grid cell (A i =0), the value remained the same.
If the activity was not present (A i =0), the function searched if the activity was present among the 8 adjacent grid cells. If this was the case, then the new value of the activity in the grid cell (A i ) was equal to half of the highest value of the activity intensity in the adjacent grid cells, otherwise the value remained at 0.
Each new application of this function allowed us to add one more grid cell in the distance function. In this case, the decrease was not linear. Let f(x) the function that gives the value of the activity intensity at a distance of x grid cells from the original grid cell, then we have: For a number of mapped benthic habitats that did not have a sensitivity assessment in the original sensitivity matrix, we calculated a sensitivity index (µ j,k ) from habitats of the same nature for which we had a sensitivity index (µ j,k ). These sensitivity indices were calculated by aggregating the sensitivity scores of these "child" habitats with two methods as described in the Supplementary Figure 5. Using the median option, we chose the median of the sensitivity value of the "child" habitats. Using the precautionary option, we chose the highest sensitivity value of the "child" habitats. At each simulation, we chose randomly between these two options.

Data Preparation
Human Activities Data As previously described, the estimation of the intensity of the physical pressures was carried out by taking into account the intensity of the human activities producing these pressures. For this purpose, relatively accurate data describing the location and intensity of human activities were collected and prepared.
The source and description of the 21 data sets describing the intensity of human activities prepared for the RCE analysis are presented in Table 4. The typology of human activities and pressures was consistent with the typology used for the implementation of the MSFD. The confidence index of the human activity data sets was evaluated as described in the Supplementary Table 4. All human activity data sets were subject to an inter-annual average covering the 2011-2016 period (2005-2013 period for the dumping of dredged material), with the exception of aquaculture, coastline artificialization and buoy mooring places. For these three activities, the data sets corresponded to the most recent data, without temporal data series. For these three data sets it was assumed that they were relatively stable and representative of the current situation. Consequently, it was assumed that all the human activities data sets represented an average situation over the recent 2011-2016 period.
Despite the good overall quality of the data collected and prepared, these data provided an incomplete estimate of the distribution and intensity of these marine activities. This was particularly true for fishing activities that are described here with data sets from the European Vessel Monitoring System (VMS) which only applies in France to vessels over 12 meters long. Therefore, all coastal fishing fleets with vessels less than 12 meters were not included in this analysis. In order to remove some outliers due to VMS data processing, fishing effort values less than or equal to 24 h per year and per grid cell were considered as zeros. For the anchoring activity, the Automatic Identification System (AIS) data set used is representative for the shipping and passenger vessels and for large yachts. However, many small pleasure craft are not equipped with AIS and are therefore highly underestimated in this data. For the aquaculture activities, the cultivated biomass by concession area was estimated from the local regulations defining the cultivation methods and conditions. These estimates gave the maximum potential biomass. They did not take into account local aquaculture contexts like unexploited or partially exploited concessions. The intensity of each human activity was calculated in the grid cells as shown in Figure 3.
The intensity of the activities in each grid cell was calculated on the basis of the relative area of activity present in each grid cell, with the exception of fishing effort data and anchored vessel numbers data which were collected directly in a regular square grid of the same resolution. The grid cells for which the presence and/or intensity of the activity are unknown take the value "No Data." The raw data sets and gridded data sets describing the human activities used in the analysis, with the exception of the fishing activities data sets, are archived and available on the Ifremer Sextant marine and coastal Geographic Data Infrastructure web portal (see the 13 references Contin, 2018a,b,c,d,e,f,g,h,i,j,k,l,m). European benthic habitat map (Populus et al., 2017) formed the basis of the data used. This Broad-scale Seabed Habitat Map (BSHM) is used as part of the implementation of the MSFD and, in some European countries, for the creation and evaluation of the marine protected area network (Andersen et al., 2018).

Benthic Habitats Data
In coastal areas, and in general wherever possible, this data set was replaced by available data with more precise spatial and typological resolution. The aim was to complete the EuSeaMap data set with data providing habitats mapped at least at the EUNIS 4 level. This step was performed using GIS tools to compile, harmonize and carry out geometric union of 150 data layers from 27 main data sources, as presented in the Supplementary Table 6.
The main steps in the construction of the multisource benthic habitats map are presented in Figure 4. The formatting of the source data (geometric formatting, reprojection, typology conversion and confidence index assignment) was carried out in the first step. The conversion of habitat typologies into the EUNIS 2007 typology was carried out using the HabRef v.4 database (Clair et al., 2017) which describes the links between the different typologies as shown in Supplementary Table 7. As presented in Supplementary Table 3, the data sources for which a habitat typology conversion was required were given the value 0 for the criterion (H_typ) contributing to the calculation of the confidence index of benthic habitat mapping (CI_H k ).
The great heterogeneity of the metadata and confidence assessment methods of the different data sources did not allow the use of pre-existing confidence indices. Some data, such as for example EuSeaMap 2017 and Carthamed, are provided with confidence indices based on complex and robust methods. Other data are provided without any confidence assessment. Moreover, it was not possible to use other methods such as the MESH project (Mapping European Seabed Habitats) confidence assessment guidelines (Foster-Smith et al., 2007) because only the data producers are able to evaluate the parameters needed for this confidence assessment method. In this context, we have developed a relatively simple method to assess a confidence index of all benthic habitat data sources as described in the Supplementary Table 3.
The multisource map generated at the end of the second step is composed of 2,077,704 polygons that cover about 375,650 km 2 with 239 different EUNIS habitat codes. The Supplementary Figure 6 presents the main features of this map. Then, the multisource benthic habitats map was integrated in the regular square grid (third step) by a SQL query presented in the Supplementary Appendix 1. Each grid cell contains the list of habitats intersecting the grid cell and the percentage of the sea surface of the grid cell covered by each habitat.
The raw and gridded versions (Quemmerais-Amice, 2020a,b) of the multisource benthic habitats map are archived and available on the Ifremer Sextant marine and coastal Geographic Data Infrastructure web portal.

RESULTS
The results presented here give an average estimate of the Risk of Cumulative Effects (RCE) of different pressures based on the intensity of human activities and associated physical pressures that took place, overall, between 2011 and 2016. The interpretation of the results therefore relates to this period. The RCE and Confidence index analyses presented here were carried out on grid cells with an average depth between 0 and −200 m. These grid cells cover the maritime space between the coastline and the slope of the continental shelf. Prior to these analyses, grid cells of which less than 50% of the area was covered by benthic habitats with a sensitivity assessment were excluded. Figure 5 shows the Risk of Cumulative Effects distribution within the French continental shelf.
Six percentage of the grid cells were excluded from the RCE calculation (n = 4591 with RCE = null, gray mesh in Figure 5). At least 50% of the surface of these grid cells was covered by benthic habitats mapped with EUNIS levels lower than 4 and/or by benthic habitats of higher EUNIS level but whose sensitivity had not been evaluated. The identification of these grid cells made it possible to locate and target the complementary studies to be carried out to map these areas and/or to complete the sensitivity assessments. The RCE was 0 in 39% of the grid cells (n = 30137 with RCE = 0, green grid cells in Figure 5). Considering the human activities and the physical pressures taken into account in this analysis, these zones did not seem to present any risk of cumulative effects for the period studied. 55% of the grid cells had a RCE value greater than 0 (n = 42842, grid cells mapped with a gradient of red). The highest RCE values (top quartile in deep red) were scattered across all the study area, both in very coastal and offshore areas. The sum of the RCE values of the grid cells included in the top quartile (n = 10710, 13% of the grid cells) represented about 58% of the total risk over the entire study area. Considering the human activities and physical pressures taken into account in the analysis, these grid cells corresponded to the areas where the RCE value was the greatest for the period studied. As shown in Figure 6, the French area included in the Celtic Seas marine subregion seems to be the least affected by the RCE, with 63% of the grid cells having RCE values between 0 (45%) and the bottom quartile and only 3 % of grid cells in the top quartile of RCE values. In contrast, the French area included in the Western Mediterranean Sea subregion is highly contrasted, with 44% of the grid cells having a RCE of 0 and 50% of the grid cells having RCE values in the fourth or top quintile. As shown in Figure 5, these grid cells concern all coastal areas and almost all of the Gulf of Lion.
An analysis of the activities and pressures contributing to the RCE throughout the study area provided a better description of management issues. As shown in the Figure 7, the bottom trawl was the most widespread activity and was present in the largest number of grid cells (approximately 50% of the grid cells), which is consistent with a previous study (Lorance et al., 2009). This activity seemed to contribute, via the generated pressures, to approximately 70% of the sum of the RCE over all the grid cells in the study area. This global vision revealed that only 2 activities (listed in descending order: bottom trawl and dredge) contributes to approximately 80% of the sum of the RCE over all the grid cells in the study area.
Similarly, the surface abrasion pressure was present in about 60% of the grid cells and contributed to about 27% of the sum of the RCE over all the grid cells in the study area (Figure 8). Only 4 physical pressures (listed in descending order: surface abrasion, reworking of sediment, light deposition and change in suspended solids) contributes to approximately 80% of the sum of the RCE over all the grid cells in the study area.
This statistical analysis, if carried out on specific areas (territorial sea, MSFD marine sub-region or marine protected areas) may reveal more local issues and inform decision-makers on more precise management issues in order to act on the pressures and activities that contribute the most to the risk.
The mapping of the RCE confidence index (Figure 9) showed a very constrained situation. The map can make it possible to localize the areas having an insufficient confidence index, and then to select the data associated with these grid cells (confidence indices are stored in the Database) to identify the data that best explains the confidence index. It will then be possible to better define the complementary studies to be carried out. This map will also allow managers to better understand the confidence that can be placed in the analysis in a fairly transparent manner and help to define acceptable uncertainty thresholds.
The analysis of the RCE variability was carried out with 10,000 Monte-Carlo simulations using the 7 criteria considered in the calculation process as described in Table 3. The 10,000 simulations were carried-out on the IFREMER DATARMOR supercomputer. Figure 10 shows the distribution of the RCE variability. The red gradient refers to the top quartile of RCE values. 15% of the grid cells (11710 grid cells) were  top quartile of the RCE values between 50 to 75% of the simulations (g2: red meshes).
The green gradient refers to the bottom quartile of RCE values. Around 12% of the grid cells were in the bottom quartile in at least 75% of the simulations (g7: deepest green). This implies that even when varying the model assumptions, these areas always have the lowest RCE values. 16% of the grid cells were in the bottom quartile of RCE values between 50 to 75% of the simulations (g6: second deep green grid cells). In areas covered by g1 and g7 grid cells (deepest red and deepest green grid cells), the RCE values could be considered as relatively stable and well evaluated with regard to the data sets and criteria taken into account in the simulations.
Twenty seven percentage of the grid cells were in the top quartile of the RCE values in at most 50% of the simulations or were in medium RCE values (g3: light pink grid cells). 12% of the grid cells were very variable (g4: yellow grid cells), sometimes in the top quartile of the RCE values (up to 12% of simulations), sometimes in the bottom quartile (up to 12% of simulations) and sometimes in the mean RCE values. Around 7% of the grid cells were in the bottom quartile of the RCE values in at most 50% of the simulations or were in medium RCE values (g5: light green grid cells).
A brief description of the distribution of the group g1 (deepest red grid cells) allows highlighting the strongest and the most visible issues for the management of the risk of cumulative effects. This analysis also highlights that taking these issues into account could be beneficial for European public policies. The group g1 covers 30% of the grid cells in the French waters of the English Channel and North Sea MSFD marine sub-region and 23% of the grid cells in the French waters of the Western Mediterranean Sea MSFD marine sub-region. These two marine sub-regions are therefore particularly concerned by the risk of cumulative effects on benthic habitats.
Bottom trawling is the activity that contributes the most to the importance of the RCE in the group g1. This activity is present in 84% of g1 grid cells and contributes for 53% of the sum of the RCE over all the grid cells of the group. This result is very probably underestimated because only the fishing vessels monitored within the framework of the European Common Fisheries Policy (VMS survey) were taken into account in the analysis. Dredge and coastal artificialization are in second and third position and contributes, respectively, to 16 and 8% of the sum of the RCE over all the grid cells of the group. As human activities generate several pressures, the relative contribution of physical pressures is more nuanced. 6 physical pressures (surface abrasion, reworking of sediment, light deposition, change in suspended solids, heavy sub-surface abrasion, light sub-surface abrasion) are present in at least 90% of the grid cells of the group and contribute, respectively, between 12 and 22% of the sum of the RCE over all the grid cells of the group.
Analysis according to biozone seabed areas, as defined by Populus et al. (2017) reveals that coastal areas are particularly affected by the risk of cumulative effects of different pressures. The g1 group cover around 27% of the grid cells which are mainly in the intertidal biozone, 33% of the grid cells which are mainly in the infralittoral biozone and 38% of the grid cells which are mainly in the shallow circalittoral biozone.
The g1 group cover 27% of the grid cells included in the marine and coastal Special Areas of Conservation designated by France within the European Council Directive 92/43/EEC (Natura 2000 network). About 75% of the marine protected areas FIGURE 7 | Contribution of human activities to the RCE. "Activity presence" represent the percentage of grid cells where the activity is located; "impacted mesh" represents the percentage of grid cells were the activity contributes to the RCE; "contribution to RCE" represents the contribution of the activity to the total RCE in the studied area, in percentage.
(97 sites out of 130) have no grid cells in the g1 group, but some sites are particularly affected, such as FR2502021 (Baie de Seine orientale) of which 72% of the grid cells are in g1, FR2500079 (Chausey) of which 71% of the grid cells are in g1 and FR7200679 (Bassin d'Arcachon et cap Ferret) of which 63% of the grid cells are in g1. 24 sites have at least 10% of their grid cells in g1.
Analysis according to benthic habitats makes it possible to target more precisely the habitats most exposed to the risk of cumulative effects. 32% of the benthic habitats (67 EUNIS codes out of 210) have at least 50% of their total area located in the g1 group. Among them some are particularly affected by the risk of cumulative effects of different pressures, such as the benthic habitats listed in the Table 5. For example, all sublittoral seagrass beds habitats (EUNIS A5.53 and higher level codes), which are of high ecological and heritage importance, have at least 76% of their total area included in the grid cells of group g1.

DISCUSSION
We have developed a method and a tool contributing to the concrete implementation of the cumulative effect assessment concept on the scale of French European waters. These developments could contribute to lay the foundations of RCE assessment  and Decision support tools (Pınarbaşi et al., 2017) needed within the framework of the MSFD and MSP European policies. They also could contribute to feed the reflections on the major assumptions and challenges related to the RCE assessment (Halpern and Fujita, 2013;Korpinen and Andersen, 2016).
One of the main interests of the proposed method is to deal concretely with the question of uncertainties and the variability of the cumulative effects assessment. These major issues related to cumulative effects assessment (Stock and Micheli, 2016;Gissi et al., 2017) are integrated in the early stages of the analysis, by providing a confidence index for each data set and for each relationship of the activity-pressure matrix and the sensitivity matrix. Confidence indices are used to map confidence in the RCE (Figure 9) and, moreover, they are integrated in the Monte Carlo simulation to map the variability of the RCE according to 7 criteria (Figure 10). In the analysis presented here, 10,000 Monte Carlo simulations were performed taking into account 7 factors that can strongly influence the cumulative effects assessment. Of the 7 factors, 3 factors (X3: value of the sensitivity score, FIGURE 8 | Contribution of pressures to the RCE. "Pressure presence" represents the percentage of grid cells where the activity is located; "impacted mesh" represents the percentage of grid cells where the activity contributes to the RCE; "contribution to RCE" represents the contribution of the activity to the total RCE in the studied area, in percentage. X4: value of the activity-pressure relationship and X5: value of the activity intensity, Table 3) directly use the confidence indices associated with the data and matrices to guide random draws on these parameters. The other criteria take into account, in particular, 3 modes of combined effects of the pressures (X1: additive, antagonistic or synergistic effects), 4 methods of calculating the intensity of the pressures based on the intensity of the activities (X2) and the distance of the pressures' effect around the point source (X6). The demonstrator tool allows one to choose the criteria (among the 7) to use in the simulations. This option will make it possible to analyze the relative contribution of each criterion to the overall variability of the model, and to define precisely the methodological efforts to be undertaken on specific criteria. Finally, other criteria can be added, such as, for example, nonlinear responses of ecosystems to pressures (Halpern and Fujita, 2013).
The data describing the human activities are associated with a sufficiently precise typology to discriminate uses and associated pressures. The fishing activity is sub-divided into 8 activities based on fishing gear. Likewise, the aquaculture is sub-divided into 7 activities according to the main breeding techniques.
In addition, the sensitivity matrix also makes it possible to hierarchize and discriminate the effect of the different pressures on each habitat. These choices make it possible to assign different weights to the pressures in the RCE assessment and ensure that pressure layers will not have roughly equal importance in the analysis (Halpern and Fujita, 2013).
The spatial resolution of the regular square grid (1 × 1 ) is consistent with the resolution of the source datasets and the scale of the study area (whole French European waters). This resolution is much finer than the scale of MSFD evaluation and management units. Even if the spatial distribution of the benthic habitats, activities and pressures is unknown within each grid cell, this resolution could give enough precision to identify the most important and structuring areas concerning the RCE at this scale. To deal with the assumption that stressors are uniformly distributed within the grid cells (Halpern and Fujita, 2013) and if the descriptive data of the various parameters are available with higher spatial resolutions, it will be possible to perform more precise analyses on specific sites such as marine protected areas.
To ensure that habitats are not simply mapped as present or absent (Halpern and Fujita, 2013), the grid cells contain the   list and the relative area of each overlapping habitat defined with the EUNIS typology. Moreover, multi-source mapping of benthic habitats (Figure 4) integrates very recent benthic habitats data sets and more than half of the habitats are mapped at least at EUNIS level 4. To deal with the normalization of the pressure intensity between 0 and 1, we chose to use the maximum values observed over the period and the study area to set the value of 1 at each pressure, as is commonly done (Halpern and Fujita, 2013).
The tool and method have been developed to integrate existing and available data, they are also adaptable to integrate new or up-to-date data and new methods of analysis, which are essential to integrate these approaches into marine public policies and are becoming step by step more operational .
There are, however, still some challenges to improve the robustness, acceptability and deployment of this approach and contribute to answer some methodological and technical issues (Halpern and Fujita, 2013;Korpinen and Andersen, 2016).

Mapping Anthropogenic Pressures
Significant evolutions are still needed to improve the mapping of anthropogenic pressures generated by marine activities. This could be achieved through a more long-term work involving marine stakeholders and academic researchers with expertise on the human activities and their interactions with the environment. This could be a dedicated project based on the mobilization of scientific expertise, scientific publications and benchmarking of different knowledge. It could concern the following issues: Evaluating the Age, the Duration and the Frequency of Activities Taking into account historical data describing human activities and their temporal accumulation is a great challenge to improve the assessment of cumulative effects (Willsteed et al., 2017). As previously stated, the descriptive data on human activities currently used for the RCE analysis covers the period 2011-2016. However, human activities that contribute the most to the risk, notably the bottom trawl, have been practiced significantly for about a century (Joubin, 1922;Lorance et al., 2009;Boisson, 2012;Eigaard et al., 2017). This activity has transformed the original habitats, leaving a footprint for several decades (Eigaard et al., 2017). Mapped habitats are already a response to past anthropogenic pressures and could be considered as the result of a relatively old anthropization of the marine environment. Currently, the RCE analysis does not take into account the historical effects which have already modified the marine environment (Korpinen and Andersen, 2016). To advance this approach, a substantial work of collecting and interpreting older data sets should be conducted. It could delimit areas that have been used regularly by certain activities for several decades. Then in these areas, the resilience capabilities (La Rivière et al., 2016) of the habitats could be assessed locally by taking into account the age and temporal persistence of pressures as proposed by Knights et al. (2015). This method would ultimately vary the sensitivity of habitats according to the local and historical context of human activities and pressures. Likewise, the duration and frequency of activities have a great influence on the resilience capacities of habitats (La Rivière et al., 2016). As proposed by Knights et al. (2015) it is possible to evaluate the duration and frequency of activities and pressures based on expert judgment.

Evaluating the Geographical Area of Influence of the Pressures for Each Type of Activity
For each activity-pressure linkage, a specific spatial model could be defined to map the pressures more realistically. This modeling can implement a linear decrease in pressure intensity (Ban et al., 2010;Andersen et al., 2013) or an exponential decrease (Holon et al., 2015). Andersen et al. (2013) have automated the use of spatial models by defining some categories of maximum potential effect distances (local, 1, 5, 10, 20, and 50 km) and by asking a group of experts to evaluate the appropriate category for the different activity-pressure pairs. Sometimes, categories of pressure intensity are defined according to the intensity of the activity. For example, for each activity, Holon et al. (2015) defined maximum potential effect distances based on activity intensity categories. Spatial models are generally applied circularly around the pressure source. Other factors such as wave or tidal currents and river plumes will significantly modify these theoretical areas of pressure influence. Hydro-sedimentary modeling can be used to model the most likely area of influence, by type of pressure and by season, for example.

Improving the Consideration of Biological and Chemical Pressures Released by Land-Based Activities
Anthropogenic contaminants and nutrients reach the marine environment mostly directly from land-based sources. Most of the chemical and bacteriological contaminants produced by these activities are introduced into the marine environment by sea outfall (water treatment plants, industrial sites) and by the hydrographic network. Some of this pollution can also be released or disseminated in the marine environment by sea-based activities (Tornero and Hanke, 2016). A recent European Environment Agency (EEA) assessment concludes that contamination by synthetic chemicals and heavy metals is still a major environmental concern in European seas (Andersen et al., 2019).
The integration of these land-based pressures into cumulative effects assessment poses several methodological difficulties. In recent studies, contaminants and nutrient inputs are estimated from the location and intensity of the land-based activities concerned: for example, agricultural area, urban area, number of inhabitants, flow of the sea outfalls (Holon et al., 2015). These data, which are not associated with real levels of defined contaminants, and the use of the sea outfall location, make it possible to estimate inputs in general. The mapping of their area of influence at sea is carried out by applying exponential decreasing equations from point source and by taking into account a reduction as a function of depth (Holon et al., 2015). Other studies have produced input modeling by contaminant category (heavy metals, synthetic compounds) and by watershed using generally the same type of data, coupled with oceanographic models to map the dilution of these contaminants (Andersen et al., 2013). Clarke Murray et al. (2015) mapped the marine influence of the dominant stressor for each land-based activity using kernel density decay at the mouth of each estuary associated with a watershed. We propose to reconsider the problem of these land-based chemical pressures in the assessment of cumulative pressures.
First, considering that the hydro-geochemical cycles of contaminants and nutrients in the watersheds are complex and very variable over time and between watersheds (Desmit et al., 2018;Le Moal et al., 2019), it is unlikely that the commonly used proxies (number of inhabitants, agricultural area) give a realistic estimate of inputs into the marine environment. For this reason, it would be interesting to directly use the contaminant concentration measurements carried out in the mouths of rivers and estuaries and at sea, in particular under the Water Framework Directive (WFD). This option would allow the estimation to be based on real concentrations by year or by season and with specific and defined groups of pollutants. It would also have a strong link with the monitoring carried out under the WFD and the descriptors 8 and 5 of the MSFD.
Secondly, significant expertise is required to choose the molecules or group of molecules on which it would be relevant and possible to work. Which molecules exert pressure on benthic habitats and communities? Considering that pollutants have very variable chemical behaviors in the sea (sedimentation, transformation, dilution), which ones can be integrated into a modeling approach? What types of concentration measurements should be used, in biota, in sediment? Should concentration thresholds be used beyond which molecules concentration acts like a pressure on habitats and benthic communities? Or should we take into account the duration of exposure to these molecules that act over the long term (number of years)? Should we define sensitivity indices for each habitat such as the sensitivity matrix used for physical pressures, or should we develop sensitivities associated with water bodies by considering their oceanographic functioning and vulnerability (WFD water bodies for example)?
Thirdly, there is a significant challenge in modeling the likely distribution of these contaminants in benthic habitats and communities. Contaminant concentration data from WFD monitoring networks can be used in oceanographic models, such as MARS3D (Lazure and Dumas, 2008) to simulate the dispersion of chemical and bacteriological contaminants by integrating the process of transformation and transfer between the water/air/sediment compartments. The use of these tools throughout the coastal zone would allow average inter-annual or seasonal concentrations for targeted contaminants to be produced and then integrated into the RCE analysis.

Data Challenges
The quality and the spatial and temporal resolution of the data as well as their typology have a very important influence on cumulative effect analysis and on the interest of the results for management. As part of this project, a significant effort was made to collect and prepare the most accurate data possible. About 40% of the project cost was dedicated to data collection and preparation. No data was purchased, it is a cost associated with the time required to prepare and manage the data. This cost is partly a consequence of the large heterogeneity of data and data sources and producers which makes it difficult to construct consistent and usable data sets. The main difficulties encountered are the heterogeneity of the data, the regulations governing their acquisition and their right of diffusion and reuse, the relatively large number of institutions involved in the production of data, and sometimes the absence of data directly responding to our needs.
Overall, it can be noted that data acquisition was facilitated by the integration of the project into the national data collection program organized for the second evaluation of the MSFD. The possible improvements concern essentially four areas: the temporal resolution, the spatial resolution, the completeness of the data already collected, and the acquisition of missing data. Most of these issues are part of the confidence index assessment parameters, as presented in Supplementary Table 4.

Temporal Resolution
The data collected to describe human activities was generally available by year, aggregated by year. This enabled the calculation of average inter-annual intensities. However, the RCE estimate would be even more interesting for management if it were calculated at least per season. This issue of temporal resolution is only very rarely taken into account in the various cumulative effect assessment studies (Korpinen and Andersen, 2016). Indeed, this implies a real qualitative leap, a more efficient information production chain and greater data preparation capabilities. In France, most descriptive data on the intensity of human activities exist or can be calculated per month, but they are difficult to access and their preparation would require much more time. For benthic habitats data, we need to produce a synthesis of the knowledge to describe the seasonal dynamics of the habitats, in particular to have a sensitivity index per season.

Spatial Resolution
In this study we propose a spatial analysis in a regular square grid of 1 min of degree by 1 min of degree of resolution. This resolution was chosen because it is the lowest spatial resolution among our datasets (fishing effort). In this context, proposing a better spatial resolution would not provide more information and would give a misleading picture of the true resolution of the data. This resolution allows, as it is, a sufficiently precise mapping of the RCE for offshore areas. For the coastal zones and especially for the Provencal and Corsica coast of the Mediterranean Sea, this resolution is insufficient. In this area, the activities are concentrated on the coastline as there is no truly continental shelf (Figures 5-7). The method and tool developed here are able to work with a much more precise regular square grid resolution, like those used by Holon et al. (2015) in the French coastal area of the Mediterranean Sea. However, it is necessary to have data that are also mapped at high resolution.

Completeness of the Data
Several data sources were crossed to produce data as complete as possible. For the extraction of marine aggregates, we crossed spatial data from Ifremer (area of dredging concessions at sea) with declarative data on annual extracted quantity by concession from regional administrations. The same principle was applied to aquaculture activities, for buoy mooring areas and for artificial reefs. For artificial reefs, the data preparation stage did not lead to sufficiently complete data sets to allow their inclusion in the analysis (complete data in the Mediterranean and very incomplete data in the Atlantic).
Other particularly important data cover only a part of the territory, such as the modeling of the winter concentration of dissolved inorganic nitrogen (DIN), which is one of the indicators of eutrophication risk. The nutrient inputs at sea, causing eutrophication phenomena, still remain a strong issue for the management of European coastal waters (Ibisch et al., 2017;Desmit et al., 2018). Very recent modeling of the winter concentration of DIN  uses the ECO-MARS3D model (Lazure and Dumas, 2008) and can only be used from the Belgian border to the Loire estuary. For this reason, we did not integrate this data in the first analysis presented in this article. This mapping may be associated with MSFD thresholds defining eutrophication to identify areas for which the DIN concentration is a pressure. Work underway under the MSFD Descriptor 8 aims to achieve these models on all French waters, including the Mediterranean.
Other activities are monitored under European and longterm policies. This is the case of fishing activities. However, French fishing effort data are also intrinsically incomplete, as only vessels over 12 meters are monitored by the VMS. The fleets of small vessels represent a significant number of vessels that work mainly in the coastal area. In this area, the current RCE results are therefore clearly underestimated. The likely evolution of European policy (European Commission, 2018) including the obligation for Member States to monitor all fishing vessels in the VMS, even small ones, will provide a much more comprehensive view of coastal fishing activities.

Acquisition of Missing Data
Some data are not yet integrated in the database and used in the first analyses, either because they do not exist or because they are particularly difficult to collect. This is the case, for example, for beach replenishment and land reclaimed from the sea. In these two examples, the data exist but are not accessible at national level. These two examples illustrate the difficulty of obtaining current data describing the construction phases of coastal developments. Other activities would require considerable work to build consistent and harmonized data sets at the national level, such as all the recreational activities that could interact with benthic habitats. These activities are not monitored at national level and there is no coordinated survey or information structuring at this scale. There is a strong challenge in coordinating the different institutions and State services that produce or manage data.
Finally, the updating process for all the data is a particularly strong issue and despite the progress made within the context of the MSFD (easier access, setting up of web portal, etc.) is not yet resolved.

Technical Challenges
Several cumulative effects analysis tools have been developed in recent years to support the implementation of marine management public policies (Pınarbaşi et al., 2017;Menegon et al., 2018). Managers require easy-to-use, turnkey analysis tools and relatively simple and fast results and interpretation. The development of user-friendly tools, taking into account their needs and their skills is a strong issue to favor the appropriation and acceptability of these methods of analysis. For now, the demonstrator tool (based mainly on a suite of opensource software) presented here is not really suitable for nonexperienced users. In the current state, its use involves opening a configuration file and the analysis script in R language as presented in Figure 2. It is not necessary to have programming skills, but it is necessary to be extremely precise and focused to follow the setup and verification procedure before starting the analysis. Similarly, visualization and exploitation of results requires opening GIS software or using R to build new graphics. Geomatics specialists of the regional marine offices of the French Office for Biodiversity have been successfully trained to use the tool. This training also made it possible to draw up a list of improvements and developments that would be interesting to consider for their local needs, in particular for the management of marine protected areas. This includes the development of a single interface for the configuration of the analysis and its launch. The setting must be facilitated. The initialization must allow the configuration to be tested before launching the analysis to quickly identify any data entry errors or inconsistent data. The management of the calculation options and the Monte-Carlo simulations must also be rationalized and allow, among other things, an estimation of the calculation time. The management of the results and the backups must also offer more options and allow all the metadata of the performed analysis (back-up of configuration options, list of parameters and data used, log files, etc.) to be stored.
Finally, this interface must also offer an application for the visualization and exploration of the results, in the form of a map, graph and statistics, allowing one to request grid cells to obtain the different information for contextualizing the result. The final tool could be used as stand-alone software or as a webtool, providing user-friendly interfaces appropriate to decision-makers and regional authorities as proposed in a recent study (Menegon et al., 2018).

Validation and Analysis of the RCE Results
The CEA method developed here gives a general diagnosis of the risk of degradation of benthic habitats. However, the interpretation of the results must remain cautious. Mapping the RCE confidence index and the variability of RCE by performing Monte-Carlo simulation allows for better location and quantification of the confidence that can be placed in the outcomes. This answers the issues identified for the assessment of the confidence of these analyses (Korpinen and Andersen, 2016;Stelzenmüller et al., 2018).
The configuration of the Monte-Carlo simulations allows the study of the relative contribution of the 7 parameters (Table 3) involved in the variability of the RCE. The developed tool makes it possible to choose the parameters used in the simulations, and also to add parameters if necessary. These complementary analyses will make it possible to identify the parameters that contribute the most to the variability of the RCE and subsequently to target future developments on these parameters. In addition, it is possible to simulate the cessation or reduction of certain human activities and related pressures, to produce RCE scenarios potentially useful for management. Lastly, the data currently available enable year-by-year RCE analyses to be conducted. The tests performed so far are based on annual averages of activity intensity. Year-by-year analyses will provide evolutionary trends of the RCE over time and space and can help to identify a reference year.
Nevertheless, it is currently impossible to measure in the field the cumulative effects of multiple pressures and to measure the relative contribution of these effects and natural variability and climate change to the status of benthic habitats. The ecological status of the habitats is the result of all these conditions.
To make progress on this issue, it will be necessary to compare the RCE result with different ecological status indicators based on in situ monitoring of benthic communities as previously done in a recent study (Clark et al., 2016). At the French scale, the evaluation of the MSFD criterion D6C5 (condition of the benthic communities) is carried-out with the BenthoVal Biotic Index (Labrune et al., 2017;Bernard et al., 2018). This indicator was initially developed to determine the impact of different sources of disturbance on soft bottom benthic communities, through the analysis of losses of individuals within affected benthic communities, compared to non-impacted reference communities. It is based on site monitoring of the abundance of benthic fauna species. For the 2018 assessment of criterion D6C5, the indicator BenthoVal 2012-2016 quantifies the loss of species abundance between the 2 years sampled during the 2012-2016 period. A significant decline in the value of the indicator indicates habitat degradation due to disturbance. However, the interpretation of this indicator is also difficult due to the difficulty of locating reference sites without pressures and also due to the difficulty of qualifying and quantifying the anthropogenic pressures on the monitored sites. For these soft bottom benthic habitats, comparing the results of the RCE and BenthoVal index would be beneficial for the interpretation of both approaches. Overall, the cross-checking between these two approaches would be a concrete application of the concept of a feedback loop (Elliott et al., 2018) to better calibrate the model and the field monitoring and ultimately increase the confidence in our results.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FQ-A was responsible of the project, conceived the manuscript, figures and tables, and wrote the manuscript. FQ-A was responsible for benthic habitats data preparation, matrices preparation, database building, and RCE analysis. JB was responsible for methodological and demonstrator tool development and R programming and contributed to the writing of the draft article. MLR contributed to the writing of the draft article. GC was responsible for human activities data preparation and contributed to the writing of the draft article. DB supported the project. All authors contributed to the article and approved the submitted version.

FUNDING
This research is a result of the Carpediem project (2016)(2017)(2018) financed by the French Office for Biodiversity.