Sediment data
North Sea
A large dataset of surface sediment OC and mud content in the North Sea (Supplementary Fig. 2) based on measurement of over 30,000 sediment samples from various German and European surveys conducted between 1960 and 2014 has been compiled and interpolated to a 0.03°×0.03° grid by Bockelmann et al.33-35, andwas used to calculate OC/mud. In addition, field data from other exiting literature and recently-collected sediment samples in both frequently trawled and untrawled areas were analyzed (Supplementary Table 1). The recently collected sediment samples include 12 stations in the Skagerrak and 23 stations in the southern North Sea (Supplementary Fig.6). Sediment and macrobenthos samples were collected from these stations by means of a Multicorer (Oktopus Kiel). Vertical distributions of TOC, TOC/TN ratio and macrobenthos biomass in the upper-most 20 cm sediments were measured following the procedure in Zhang et al.21 In the southern North Sea stations, bioturbation rates were mostly estimated by evaluation of Chlorophyll profiles. Chlorophyll was extracted from freeze-dried sediment with 90 % acetone and measured subsequently. The depth-averaged chlorophyll gradient was then evaluated with the Fick’s first law of diffusion. In the Skagerrak stations, bioturbation rates were estimated by simulating the depth distribution of 210Pb using a steady-state transport-reaction model20.
Global shelf seas
Global surface mud and OC content values (Supplementary Fig.7a, b) were compiled in the sediment database dbSEABED (https://instaar.colorado.edu/~jenkinsc/dbseabed/), which is based on over 5 million described samples at 3 million sites worldwide and mapped to a 0.1°×0.1° grid using a 3D-inverse-distance-weighted interpolation, including water depth as a vertical coordinate36,37. This ensures highly local representation of actual seabed properties.
In addition, the global surface sediment OC product by Atwood et al.16, which are based on sampling from 11,578 global sediment cores and gridded at 1×1 km2 resolution using a random forest regression model (Supplementary Fig.7c), was also included in the analysis. In order to compare the two OC datasets, the carbon stock data Cstock (in Mg C km-2) in Atwood et al.16 were re-gridded to 0.1°×0.1° and transformed to carbon content Cperc (in percentage) by inverting their applied relationship between stock and dry bulk density, yielding:
Both datasets have particularly high point densities of measurement in the nearshore and on the continental shelf, which overlap strongly with areas of trawling activity. A general consistency (r = 0.57) is shown in both OC datasets for shelf seas, and supports the analysis of OC/mud and its relationship with trawling intensity (Supplementary Figs.2-3).
Bottom Trawling Data
North Sea trawling effort
Annual aggregated spatial data on bottom trawling in the North Sea including gear-type information38 were combined with daily fishing effort data10 in order to generate a daily time-series of fishing effort on a 0.03°×0.03° grid. The annual aggregated data were used for analysis of the linkage between trawling density and OC/mud, while the daily time series served as input for numerical modeling.
A daily time series of bottom trawling impact was generated based on the dataset of Kroodsma et al.10, which contains daily vessel locations and fishing effort, as well as estimates of vessel power, length and class, but does not distinguish between specific gear types. For each vessel designated as a “trawler”, we assigned a gear type according the dominant "métier", following the definition by Eigaard et al.39 The relevant metier at each vessel’s average position was set according to the information in ICES38. A "métier" groups fishing trips by gear type and target species, such that vessels operating in the same métier are expected to have similar impacts on the seabed per unit area contacted. A total of 14 métiers have been defined, of which eight are otter trawlers, three are beam trawlers, two are demersal seines and one is dredges. Among all métiers in the North Sea, otter and beam trawlers are dominant and account for more than 90% of total annual fishing effort38. The seine ropes of demersal seines are expected to have no subsurface impacts. Dredges are expected to have a high impact per surface area contacted40, but they are utilized mainly off the British coast and do not have a large spatial coverage38.
Eigaard et al.39 provided empirical formula to estimate gear widths from vessel length or vessel power, as well as average towing speeds and the length proportions of gear components for each métier. To avoid extrapolating to an unrealistic distance too far from the data points in Eigaard et al.39, we limited the gear width to a maximum of 300 m, thus preventing excessively large gear widths which would otherwise occur in less than 5% of all vessels.
For each grid cell, the daily surface swept area ratio (SAR) is calculated as:
where W is the width of the gear or gear component, ν is the vessel speed, T is the time trawled on that day, and A is the total grid cell area. We estimated the daily SAR in three depth intervals within the seabed: 0-2 cm, 2-5 cm, and 5-10 cm. We distinguished between muddy and sandy areas when computing the SAR at different penetration depths. The depth to which the individual gear components penetrate in sandy and muddy sediment were determined using information given in Eigaard et al.39(Supplementary Table 2). The silt fraction of the sediment is assigned according to the surface sediment map provided in Bockelmann et al35.
A comparison between the annual aggregated bottom trawling effort from ICES, which are based on the Vessel Monitoring System (VMS) 38 and logbook data, and daily fishing effort data, which are based on the self-reported Automatic Identification System (AIS)10, shows a consistency for recent years (2015-2020, Supplementary Fig.8). The annual aggregated trawling effort from the AIS-based data are 17%, 16% and 8% lower than the ICES data for 2015-2017, respectively. Starting from 2018 both datasets yield comparable results within ±5% difference. However, the AIS-based data underestimate the bottom trawling effort in early years (e.g. 424%, 46% and 25% lower than the ICES data for the year 2012 to 2014, respectively). This deficit impedes the use of the AIS-based data as input for numerical modeling for long-term (decadal) historical hindcast. In order to provide a consistent daily-resolved long-term time series of trawling data prior to 2015 for numerical modeling, we averaged the daily SAR fields of 2015-2020 and scaled them with reported historical data. For the period of 1985-2014, a historical reconstruction of annual aggregated bottom trawling effort has been compiled by Couce et al.41 The outcomes of Couce et al.41 suggest a fairly stable distribution of fishing effort among countries, seasons and gear types for the 30-years period. Such stable distribution is also seen in the ICES 2015-2020 data38 (see also Supplementary Figs.9-11). Based on this stability, the spatial and daily distribution of fishing effort in the historical period 1950-2014 is assumed to follow the same spatial and seasonal pattern as in 2015-2020, but is scaled by the ratio of the annual aggregated trawling effort for that specific year to the average of the years 2015-2020. For the period of 1950-1984, in which annual aggregated trawling data are missing, we employed three scaling methods to generate the spatio-temporal distribution of daily effort. Method-1 is based on the relationship between trawling effort and total annual landings. The time series of bottom trawling effort for 1985-2020 shows a high correlation (Pearson's r = 0.86) with the annual landing data of demersal and benthic fish (ICES42; Supplementary Fig.12). Based on this correlation, we extrapolated the time series of annual aggregated bottom trawling effort for 1950-1984 through a linear regression between the trawling and landing data. Method-2 is based on the relationship between the total trawling effort in the North Sea and the UK’s portion. Among all North Sea countries, trawling data from the UK have the longest temporal coverage dating back to the 1910s43. The trawling effort of UK vessels shows an almost linear (r = 0.98) relationship with the total effort in the North Sea in 2003-2020. Based on this relationship, the total trawling effort for 1950-1984 was assumed to also scale linearly with the UK trawling effort. In the Method-3, we considered an increase of technological efficiency, or “creep factor,” which is used in fisheries science to adjust for the gradual increase in the effectiveness of fishing gear resulting from the successive introduction of technological improvement to fishing gear and vessels45. A technological creep factor of 2.4% per year45 was adopted in scaling the trawling effort with the total annual landings. The resulting time series from the three scaling methods (Supplementary Fig.12) help to evaluate the sensitivity of sedimentary OC stock to long-term temporal variability of trawling effort.
Global trawling effort
Global annual surface SAR (Fig. 1a) is estimated using the dataset of Kroodsma et al.10 For the estimate of gear width and vessel speed using the relationships in Eigaard et al.39, it is assumed that all bottom trawlers are beam trawlers targeting demersal fish (métier "TBB_DMF"). Only areas with depths less than 300 m are considered. The annual fields are averaged for the years 2015-2020. The overall pattern is similar to previous global trawling intensity estimates 8,46, with highest values in the Northwestern European shelf (North and Baltic Seas), parts of the Mediterranean (e.g. Adriatic Sea) and the Chinese shelf seas. One notable discrepancy of our SAR estimate compared with the qualitative trawling intensity map of Oberle et al.46 is that our data show little trawling activity in Southeast Asia (e.g. Sunda shelf), which is due to the low coverage of AIS data in this region47.
Numerical Modeling
We applied a 3-dimensional physical-biogeochemical model to simulate the impact of bottom trawling on macrobenthos and OC cycling in surface sediments. The model consists of three major components, with HAMSOM (HAMburg Shelf Ocean Model48) for hydrodynamics, ECOSMO (ECOSystem Model49) for biogeochemistry and TOCMAIM (Total Organic Carbon-MAcrobenthos Interaction Model20) for benthic fauna. OC is divided into three pools depending on the degradability, namely labile (i.e., of high nutritional quality for macrobenthos with first-order oxic remineralization rate k = 20 yr-1), semi-labile (i.e., of intermediate nutritional quality with k = 2 yr-1), and refractory (i.e., of low nutritional quality with k = 0.01 yr-1). With a Neumann boundary condition of OC fluxes at the sediment-water interface calculated by the pelagic model component, sedimentary OC content is solved by a mass balance equation including the impacts of deposition/erosion, oxygen-dependent first-order degradation, macrobenthic uptake, respiration and bioturbation21,50. Temporal change of macrobenthic biomass is calculated based on food availability in the form of OC, temperature, oxygen, and mortality caused by predation and bottom trawling. Bioturbation scales with macrobenthic biomass and is inversely related with local OC resource20.
The model has been applied to the North Sea to investigate the spatio-temporal variability of macrobenthic biomass20,21, sediment OC storage21 and benthic oxygen consumption50. Details of mathematical descriptions, sensitivity analysis of model parameters and application to station data are provided in Zhang & Wirtz20. Model coupling between benthic and pelagic components as well as its application to the North Sea are introduced in Zhang et al.21,50.
The long-term hindcast simulation encompasses the period 1950 – 2016. The map of surface sediment OC content (Supplementary Fig.1b) is used as initial condition for modeling assuming that all OC is refractory. Input of labile and semi-labile OC to sediment is provided by the pelagic model component which simulates the primary production and associated OC fluxes during the period 1950 – 2016 21. Because macrobenthic biomass is dependent on both quality and quantity of sedimentary OC, simulation of the first three decades is featured by an increase of biomass due to gradual accumulation of labile and semi-labile OC in sediments. Afterwards the simulated macrobenthic biomass and OC fluxes across the sediment-water interface approach a relatively stable level21. Simulation results for the period 1985-2016 were used for evaluation of the benthic OC fluxes (Fig.2c & Supplementary Fig.4).
Numerical implementation of trawling impacts
The impact of bottom trawling on OC storage in sediments is implemented through enhanced OC resuspension and subsequent transport, physical mixing in the penetrated sediments and enhanced mortality of macrobenthos, which subsequently affects macrobenthic uptake of OC, respiration and bioturbation intensity.
In order to estimate the resuspension rate of each trawler in the daily resolved data, we follow the study of O'Neill & Ivanović51, which demonstrated that sediment entrained behind towed fishing gear is directly related to the hydrodynamic drag of the fishing gear components and to the seafloor sediment type. Their empirical formula for the amount of sediment mobilized per contacted area (in kg m-²) is:
Maximum resuspension per area contacted is limited to 6 kg m-² in order to avoid excessive resuspension beyond reported values46, which pertains to about 5% of all vessels. For the year 2017, the average resuspension per area contacted calculated in this way is 2.1 kg/m² with a standard deviation of 1.6 kg/m², which is well within the range of reported values. The associated average resuspension rate is 213 kg s-1, with 48% of vessels having an erosion rate below 100 kg s-1 and 88% below 500 kg s-1.
The resuspension is divided among mud and sand according to their fraction in the seabed (Supplementary Fig.1a). Accordingly, OC in surface sediment is also resuspended according to their weight proportion in the sediment. At each model time step, the trawling-induced resuspension flux is added to the hydrodynamic resuspension flux.
For implementation of trawling impact on benthic fauna, three impact levels in mortality, namely low (10th percentile, corresponding to depletion rate of 0.11 for SAR=1), medium (50th percentile / depletion rate of 0.2) and high (90th percentile / depletion rate of 0.3), were included according to field assessments7, 44. Trawling-induced daily mortality is then calculated by:
where I = 0.11/ 0.2/ 0.3 refers to the low/medium/high mortality rate, respectively. M is divided into three depth layers (0-2 cm, 2-5 cm and 5-10 cm) according to the corresponding SAR values in these zones (Supplementary Table 2).
In addition, the impact of physical mixing was implemented as a mixing coefficient13 scaled by SAR in the impacted sediment depth layers. Due to lack of data support for determining the physical mixing coefficient in field, we adopted two constant values, namely 2.4 and 0.24 cm2 day1 representing high and low mixing efficiency, respectively, in the model experiments. The mixing coefficient is multiplied by the daily SAR to calculate daily matter fluxes (OC and oxygen).
The three generated time series of trawling, combined with the different impact levels of trawling-induced faunal mortality and physical mixing, produce in total 18 different scenarios (Supplementary Table 3). These scenarios were simulated and compared to a non-trawling scenario (No-trawling).
Validation of numerical modeling
Our simulation results have been confirmed by field data of OC, oxygen consumption rate and macrobenthic biomass in seafloor surface sediments of the southern North Sea21,50. In this study, the model domain was extended to the entire North Sea and additional field data from the northern part (Supplementary Table 1) were compiled. The measured OC data, macrobenthos biomass and estimated bioturbation rates at these field stations were used to assess the simulation results.
A comparison between the simulated mean values of the 18 experiments with trawling impact and field data suggests a generally satisfactory model performance in capturing spatial distribution of OC content (Fig.3b), macrobenthic biomass and bioturbation rates (Supplementary Fig.13). The correlation coefficient (Pearson’s r) between the measured and simulated values for macrobenthic biomass (in ash free dry weight) and bioturbation rate is 0.91, and 0.71, respectively. The Root Mean Square Error (RMSE) for biomass is 2.23 g m-2, equivalent to 40% of the spatially-averaged biomass in the North Sea. The simulated spatial distribution of biomass along latitude (Supplementary Fig.13b) also shows consistency with data from two large-scale surveys56. The Normalized RMSE for bioturbation rate is 0.176. It is worth noting that the estimated bioturbation rates from fitting tracer profiles (chlorophyll and 210Pb) vary over one order of magnitude within the samples from each field station, confirming a large small-scale heterogeneity of bioturbation driven by variations in food supply20, metabolism rate of benthic fauna57 as well as sediment composition and community structure21,58. Despite that the model grid (0.03°×0.03°) is too coarse to resolve such heterogeneity at small spatial scales, the mean bioturbation rate from the samples at each station was reasonably captured by the simulations.
- Bockelmann, F. D. Mud content of Noarth Sea surface sediments. World Data Center for Climate (WDCC) at DKRZ (2017). https://doi.org/10.1594/WDCC/coastMap_Substrate_Mud
- Bockelmann, F.D. Total organic carbon content of North Sea surface sediments. World Data Center for Climate (WDCC) at DKRZ (2017). https://doi.org/10.1594/WDCC/coastMap_Substrate_TOC
- Bockelmann, F.D., Puls, W., Kleeberg, U., Müller, D., & Emeis, K.C. Mapping mud content and median grain-size of North Sea sediments – A geostatistical approach. Mar. Geol. 397, 60–71 (2018).
- Jenkins, C. Building offshore soils databases. Sea Technol. 38, 25–28 (1997).
- Bostock, H. et al. Distribution of surficial sediments in the ocean around New Zealand/Aotearoa. Part B: continental shelf. N. Z. J. Geol. Geophys. 62, 24–45 (2019).
- ICES. OSPAR request 2018 for spatial data layers of fishing intensity/pressure. Data Outputs (2019). https://doi.org/10.17895/ices.data.4686.
- Eigaard, O.R. et al. Estimating seabed pressure from demersal trawls, seines, and dredges based on gear design and dimensions. ICES J. Mar. Sci. 73, i27-i43 (2016).
- O'Neill, F.G., Robertson, M., Summerbell, K., Breen, M., & Robinson, L.A. The mobilisation of sediment and benthic infauna by scallop dredges. Mar. Environ. Res. 90, 104–112 (2013).
- Couce, E., Schratzberger, M., & Engelhard, G. H. Reconstructing three decades of total international trawling effort in the North Sea. Earth Syst. Sci. Data 12, 373–386 (2020).
- ICES. Greater North Sea Ecoregion Fisheries Overview - Data Output file. Data Outputs (2021). https://doi.org/10.17895/ices.data.9158.
- Engelhard, G. H. North Sea Historical Effort by Rectangle by Month: 1913 -1980. Cefas, UK. V1 (2016). https://doi.org/10.14466/CefasDataHub.24
- Bergman, M.J.N., & Van Santbrink, J. W. Mortality in Megafaunal Benthic Populations Caused by Trawl Fisheries on the Dutch Continental Shelf in the North Sea in 1994. ICES. J. Mar. Sci. 57 (5), 1321–1331(2000).
- Palomares, M. L. D., & Pauly, D. On the creeping increase of vessels’ fishing power. Ecol. Soc. 24(3), 31 (2019).
- Oberle, F.K.J., Puig, P., Martín, J. Fishing Activities, In: Micallef, A., Krastel, S., Savini, A. (Eds.), Submarine Geomorphology. Springer International Publishing, Cham, pp. 503–534 (2018). https://doi.org/10.1007/978-3-319-57852-1_25.
- Taconet, M., Kroodsma, D., Fernandes, J. Global Atlas of AIS-based fishing activity. FAO, Rome (2019).
- Schrum, C. Thermohaline stratification and instabilities at tidal mixing fronts - Results of an eddy resolving model for the German Bight. Cont. Shelf Res. 17, 689-716 (1997).
- Daewel, U., & Schrum, C. Simulating long‐term dynamics of the coupled North Sea and Baltic Sea ecosystem with ECOSMO II: Model description and validation. J. Mar. Syst. 119–120, 30–49 (2013).
- Zhang, W. et al. Quantifying importance of macrobenthos for benthic-pelagic coupling in a temperate coastal shelf sea. J. Geophys. Res. Oceans 126, e2020JC016995 (2021). https://doi.org/10.1029/2020JC016995
- O'Neill, F.G., & Ivanović, A. The physical impact of towed demersal fishing gears on soft sediments. ICES. J. Mar. Sci. 73 (1), i5–i14 (2016). https://doi.org/10.1093/icesjms/fsv125
- O'Neill, F.G., Noack, T. The geometry and dynamics of Danish anchor seine ropes on the seabed. ICES. J. Mar. Sci. 78, 125–133 (2021).
- O’Neill, F.G. & Summerbell, K. The hydrodynamic drag and the mobilisation of sediment into the water column of towed fishing gear components. J. Mar. Syst. 164, 76 – 84 (2016).
- O'Neill, F.G., Noack, T. The geometry and dynamics of Danish anchor seine ropes on the seabed. ICES. J. Mar. Sci. 78, 125–133 (2021).
- Solan, M. et al. Worldwide measurements of bioturbation intensity, ventilation rate, and the mixing depth of marine sediments. Sci Data 6, 58 (2019).
- ICES. Structure and dynamics of the North Sea benthos. ICES Cooperative Research Reports (CRR). Report (2007). https://doi.org/10.17895/ices.pub.5451
- Teal, L. R., Bulling, M. T., Parker, E. R., & Solan, M. Global patterns of bioturbation intensity and mixed depth of marine soft sediments. Aquat. Biol. 2(3), 207–218 (2008).
- Morys, C., Forster, F., & Graf, G. Variability of bioturbation in various sediment types and on different spatial scales in the Southwestern Baltic Sea. Mar. Ecol. Prog. Ser. 557, 31–49 (2016).
- Ståhl, H. et al. Factors influencing organic carbon recycling and burial in Skagerrak sediments. J. Mar. Res. 62, 867-907 (2004).
- Rosenberg, R., Hellman, B., & Lundberg, A. Benthic macrofaunal community structure in the Norwegian Trench, deep skaggerrak. J. Sea Res. 35 (1-3), 181-188 (1996).