Cumulative Human Impacts on Mediterranean and Black Sea Marine Ecosystems: Assessing Current Pressures and Opportunities

Management of marine ecosystems requires spatial information on current impacts. In several marine regions, including the Mediterranean and Black Sea, legal mandates and agreements to implement ecosystem-based management and spatial plans provide new opportunities to balance uses and protection of marine ecosystems. Analyses of the intensity and distribution of cumulative impacts of human activities directly connected to the ecological goals of these policy efforts are critically needed. Quantification and mapping of the cumulative impact of 22 drivers to 17 marine ecosystems reveals that 20% of the entire basin and 60–99% of the territorial waters of EU member states are heavily impacted, with high human impact occurring in all ecoregions and territorial waters. Less than 1% of these regions are relatively unaffected. This high impact results from multiple drivers, rather than one individual use or stressor, with climatic drivers (increasing temperature and UV, and acidification), demersal fishing, ship traffic, and, in coastal areas, pollution from land accounting for a majority of cumulative impacts. These results show that coordinated management of key areas and activities could significantly improve the condition of these marine ecosystems.


General Approach
We used a cumulative impact model that follows a 4-step process (Halpern et al. 2008;Selkoe et al. 2009). We first assembled the spatial data for each anthropogenic driver (D i ) and each ecosystem (E j ). Second, we log[X+1]-transformed and rescaled between 0-1 each driver layer to put them on a single, unitless scale that allows direct comparison, and converted ecosystem data into 1 km 2 presence/absence layers. Third, for each 1 km 2 cell of ocean we multiplied each driver layer with each ecosystem layer to create driver-by-ecosystem combinations, and then multiplied these combinations by the appropriate weighting variable (u ij ).
These weighting variables come from an expert survey that assessed the vulnerability of each ecosystem to each driver on the basis of 5 ecological traits (Halpern et al. 2007). The weighting values represent the relative impact of an anthropogenic driver on an ecosystem within a given cell when both exist in that cell, and do not represent the relative global impact of a driver or the overall status of an ecosystem. The sum of these weighted driver-by-ecosystem combinations then represents the relative cumulative impact of human activities on all ecosystems in a particular 1 km 2 cell. Finally, we used the same thresholds used in the global analysis to designate ecologically meaningful categories of the relative cumulative impact scores. In the global analysis, these were based on empirical data on the condition of coral reef ecosystems (Halpern et al. 2008). The sensitivity of our results to key steps in this process and further details on the groundtruthing method was analysed in previous articles (Halpern et al. 2008;).
In the present analysis we replaced some of the datalayers and included additional data to better reflect the specific pressures and ecosystems of the Mediterranean basin. A total of 22 spatial datasets of drivers (Table S1) and 17 ecosystem types (all the ecosystem layers used in Halpern et 2 al. 2008 except for those that do not exist in the Mediterranean Sea, i.e. the coral, mangrove and ice ecosystems) were assembled and used in the analyses and maps.
Some of the dalalayers measure impacts that are to some extent overlapping. In particular, nutrient input increases the risk of hypoxia, and these datalayers could be viewed as measuring the same impact on coastal marine ecosystem. However, each measures also distinct impacts.
While nutrient runoff is a major contributor to coastal hypoxia, other factors also determine the occurrence of hypoxia (e.g. hydrodynamics, primary productivity and water temperature, all affecting the formation and degradation of organic matter). These factors are included in the OxyRisk model that generated our risk of hypoxia datalayer (Djavidnia et al. 2005). Therefore, coastal areas receiving high nutrient loadings do not necessarily develop hypoxic waters, depending on the physical and biological setting. However, high nutrient loadings can have other impacts, beyond hypoxia, including altering primary productivity, shifting the composition of planktonic and macroalgal assemblages, and affecting benthic invertebrate assemblages and associated food webs. Thus, we included both datalayers (risk of hypoxia and nutrient input) in our analysis because they capture complementary impacts of increased nutrient availability.
Similarly, coastal population density and urbanization trends result in increased nutrient loading, but also in a suite of direct and indirect ecological impacts, through direct disturbance of people visitation and use of the shores (e.g., trampling, noise, collecting, littering) and habitat modification (conversion of natural habitat to hard structures, removal of coastal vegetation).
Moreover coastal engineering results in the construction of structures that modify coastal circulation and sediment dynamics, in addition to replacing natural habitats. So in this case as well, we included all of these datalayers because they capture different impacts of coastal ecosystems, in addition to overlapping ones.
3 Maps of individual data layers and of the full cumulative impact model can be viewed at the public website: http://globalmarine.nceas.ucsb.edu/mediterranean/ We conducted both a Mediterranean-wide analysis (displayed on the website above) and a more detailed analysis of cumulative impacts to the territorial waters (within 12 nm from the coastline) of Mediterranean and Black Sea countries that belong to the European Union (EU). This analysis used all of the data layers used in the first analysis and four additional layers (Table S1) that were available only for the EU member states.
A detailed description of data layers obtained from the previous global analysis can be found in the SOM of Halpern et al. 2008. We attempted to obtain or model measures of fishing efforts from commercial fleets and artisanal fisheries, but ended up using commercial fisheries catch statistics instead because actual effort data are available only for pelagic longline fisheries (data provided by one of the authors, Rebecca Lewison, UCSD, USA) but not for other fisheries. For pelagic longlines, we compared effort and catch, after aggregating datasets on the same coarse scale (5 degree cells) at which the effort data are reported (Fig. S1). Although visual inspection reveals a weak trend towards a positive correlation between effort and catch, this relationship is not statistically significant (r=0.33, df=14, ns), indicating that, for these fisheries and likely for others, catch data are not good proxies for fishing effort. Although we used catch as the only available measure of how fishing pressure may be distributed across the region, developing reliable estimates of fishing effort across the Mediterranean is a key research priority.
Below we describe the additional data layers used in these analyses. Halpern et al. (2008) modeled the incidence of invasive species as a function of the amount of cargo traffic at a port, on the basis of results from other studies showing a relationship between these two variables and in the absence of actual data for the global distribution of invasive species (Fig. S2). In this analysis, we replaced this layer with data on the actual distribution of a subset of invasive species in the Mediterranean (Fig. S3). We used maps compiled by The Mediterranean Science Commission (CIESM Atlas of Exotic Species, www.ciesm.org) on the distribution of exotic crustaceans, fish, molluscs and macroalgae in the Mediterranean. CIESM provided low resolution, per-species images showing the area impacted by each invasive species. From these images we derived spatial maps of species distribution by georeferencing the coordinates contained within the images with our study area, which converts the image coordinates into geographic coordinates via a transformation function. We used the GDAL software (http://trac.osgeo.org/gdal/wiki/FAQGeneral) with seven control points and a second order polynomial fit. These georeferenced maps were then digitized to capture the areas of frequent, infrequent and local occurrence. Finally, these polygons were rasterized to the resolution and extent of the overall model.

5
We overlayed the distributions of each of these species to produce a datalayer of the number of harmful marine invasives reported within each 1 km 2 of intertidal and shallow subtidal (down to 60 m depths) coastal habitat.

Oil spills
Pollution at sea, originating from shipping (e.g., oil spills, waste discards, oils and fuels leakage) was modeled in Halpern et al (2008) based on the data on shipping traffic (see above, and details provided in Halpern et al. 2008) (Fig. S4). Here we replaced this data layer with data on the actual occurrence and magnitude of oil spills within the Mediterranean. We obtained data on the occurrence and magnitude of ship accidents resulting in oil spills between 1977-2009, compiled by the Regional Marine Pollution Emergency Response Centre for Mediterranean Sea (REMPEC) (Fig. S5). Using the included attribute of oil released quantity, we performed the same plume modeling techniques used for our land-based threats (see above, 'land based drivers') to distribute the quantities into the surrounding ocean waters. For those accidents where no release quantity was provided, we used a default value of 1 ton released. A total of 93,600 km 2 was affected by one or more accidents. Halpern et al. (2008) quantified temperature warming as the difference in the frequency of positive temperature anomalies between the first and last 5 years of a twenty-year time series  of sea surface temperature (SST) data from remote sensing (Fig. S6). In the present analysis we replaced this data layer with estimates of rates of change in SST as these represent a more direct measure of temperature warming. The detailed spatial analysis of rates of change in SST, utilizing remote sensing data for 1985-2006, is taken from Nykjaer (2009). As shown in Nykjaer (2009) the change in temperature is not constant throughout the year but occur primarily 6 during summer months. In order to capture the seasonal influence of SST we repeated all analyses using either the annual or the summer (calculated for the months of June, July and August) SST rate of change. However, the resulting cumulative impact maps using either the mean annual or summer SST rate of change were very similar. The map of the resulting annual SST anomalies for the entire area and selected sub areas is shown in figure S7.

Risk of hypoxia
The semi-enclosed setting of the Mediterranean and the high coastal population density make coastal areas of the Mediterranean Sea prone to anthropogenic eutrophication and associated formation of hypoxic and anoxic waters. To capture this pressure on Mediterranean coastal marine ecosystems, we included in our analysis an additional data layers representing the risk of development of hypoxic conditions. The risk of development of low-oxygen conditions was predicted based on a biophysical model (OxyRisk) developed by Djavidnia et al. (2005). We used the modeled output to identify areas susceptible to hypoxia (Fig. S8).

Coastal erosion, aggradation, and armoring.
Data on coastal erosion, aggradation (Fig. S9), and armoring ( Fig. S10) were available, only for the Mediterranean countries that belong to the EU, from the EU project CORINE managed by the European Environmental Agency. This program was set up with the aim to gather and integrate information on the state of the environment in the European Community, as a basic instrument for the European Environmental policy. The concrete goal of a component of this project, CORINEcoastal erosion (CCEr), was to provide a description of European coastal systems' morphology and highlight the causes of its modification. Such causes fall under two nonexclusive categories: (i) major natural events (subsidence and rising of the marine level; storms and marine overhigh tides, seisms, mass movements); (ii) human actions causing movements of sediments (reduction of river contributions; development of estuaries, engineering of the coast; construction on the coastal dunes; harbor work and coastal defense construction; sediment, water, gas or oil extraction).
Data layers are downloadable free of charge from the EEA web site: http://dataservice.eea.eu.int/dataservice/metadetails.asp?table=COASTEROS&i=1 The data is provided as a series of line segments of the EU coastline, and includes attributes on each segment indicating the level of erosion, aggradation and coastal defenses present within that line segment. For our purposes, we classified the erosion and aggradation classes based on the certainties included in the attribute data, with confirmed reports of erosion and aggradation taking priority over undocumented reports. These classes were then used to estimate the extent and impact of these drivers, which resulted in: 33,807 km 2 aggrading, 80,259 km 2 eroding, and 71950 km 2 undergoing defense.
Datalayers on the presence of harbors or coastal defense structures, and on the erosion or aggradation of the coastline are delivered as vector format with a scale 1:100,000.

Urbanization trends
Coastal urbanization trends (Fig. S11) data were also produced by the EU project CORINE. This datalayer was downloaded from http://www.eea.europa.eu/data-and-maps/data/urban- 122 (Road and rail networks) and 511 (Water courses), when neighbors to the enlarged core classes, cut by 300m buffer. Forests & scrub (311,312,313,322,323,324), when they are completely within the core classes.
An UMZ layer was obtained, which was rasterized (100 m resolution).

Normalization of driver data layers
Nearly all of the anthropogenic driver data had extreme left-skewed distributions with very small numbers of extremely high outlier data. The extreme outlier data may or may not be precisely accurate, and so we sought to minimize the impact of this potential source of error on our model results. To do so we log[X+1]-transformed each data layer, except benthic structures. Benthic structures were treated as binary data since presence of an oil rig is an all-or-nothing event. All data layers were then rescaled between 0-1, with the highest log-transformed value for each driver set = 1. The transformation of data appropriately reduces the effect of extreme outliers when rescaling the data to assign the relative impact of different levels of the anthropogenic drivers considered here (Malczewski 2000), while rescaling allows for direct comparison among drivers with dramatically different native scales and units of impact. We also rescaled all driver layers without log-transformation of the data for comparison. This method preserves the true relative magnitude of the data; we found no qualitative differences in our results (not shown) and 9 so focus on those that derive from log-transformed data. Our approach assumed a linear relationship between driver magnitude and impact on ecosystems. This assumption ignores thresholds that likely exist but are known for very few driver-by-ecosystem combinations, but it allows for direct comparison between and among drivers.

Data representation and projection
To create a uniform coastline (land-sea interface) against which we compared all driver and ecosystem data, we used the SRTM30-PLUS data (http://topex.ucsd.edu/WWW_html/srtm30_plus.html) that represents merged SRTM30 and ETOPO2 data. Land-based data that occurred within the ocean or ocean-based data that occurred on land, in reference to this data mask, were clipped and removed. This issue emerges solely for coastal cells, where a 1km 2 region has both land and ocean in it but must be assigned a binary value (land or ocean) or where the two layers do not meet (i.e. a gap exists). Our approach of applying a uniform land mask ensured consistency across all data used in our analyses, but removed or added area from coastal marine ecosystems. Shallow soft bottom and rocky reef ecosystems gained area as a result of this process due to the need to fill gaps between land and ocean datalayers -we expanded neighboring ocean substrate data into these gaps. In contrast, seagrass ecosystems lost area since these habitats are in very shallow water and so were more prone to being classified as land.
All data were represented at 1km 2 resolution, even though several layers had native resolutions at coarser scales. In doing so, we assumed the coarse-scale value was evenly distributed across all 1km 2 cells within that region. For the climate change drivers (UV anomalies, SST change, and changes in ocean acidification), this assumption is reasonable given the scale at which those drivers act. The land-based drivers, human population data, and oil/gas development data were all at our native 1 km 2 resolution, and spread of these data into the ocean at the same resolution is reasonable. Regardless, when coarse-scale data are distributed equally to all 1 km 2 cells within the larger area, the coarser scale pattern is essentially recreated while the finer resolution information is preserved where and when it is appropriate. Finally, prior to all analyses, we converted all data to the WGS84 Mollweide projection as it is an accurate single global projection that preserves area and allows data transfer and analysis among operating systems and software. Table S1. List and characteristics of spatial datalayers representing drivers of change of marine ecosystems used to assess and map cumulative human impacts to the Mediterranean and Black Seas. Data used were obtained from the previous global analysis conducted by Halpern et al. (2008) Figure S1. Correlation between fishing effort (in No. Hooks/year, for the years 1999-2000) and catch (in metric tons/year, normalized to local primary productivity, see Halpern et al. 2008) for pelagic longlines. Figure S2. Modeled distribution of marine invasive species. This layer was used in the global analysis by Halpern et al. (2008), and was replaced here with a data layer on the actual distribution of 20 harmful invasive species (Fig. S3). Figure S3. New data layer used in these analyses on the distribution of harmful invasive species in Mediterranean coastal waters. Figure S4. Modeled pollution at sea layer used in the global analysis by Halpern et al. (2008).
This layer was replaced with actual data on the occurrence and magnitude of oil spills in the Mediterranean (Fig. S5). Figure S5. Occurrence and magnitude of accidents resulting in oil release (data from REMPEC).
This data layer was used in the present analyses. Figure S6. Frequency of positive temperature anomalies used in the global analysis by Halpern et al. (2008). In the present analysis, this data layer was replaced with rates of sea surface temperature increase (Fig. S7). Figure S7. Annual average rate of increase in SST (between 1985SST (between -2006).
18 Figure S8. Risk of hypoxia. Figure S9. Coastal erosion (red) and aggradation (green) along the coastlines of Mediterranean EU countries. Figure S10. Coastal armoring (black) along the coastlines of Mediterranean EU countries. Figure S11. Urbanization trends within Mediterranean EU countries.