Impact of bottom trawling on sediment biogeochemistry: a modelling approach

Bottom trawling in shelf seas can occur more than 10 times per year for a given location. This affects the benthic metabolism, through a mortality of the macrofauna, resuspension of organic matter from the sediment, and alterations of the physical sediment structure. However, the trawling impacts on organic carbon mineralization and associated processes are not well known. Using a modelling approach, the effects of increasing trawling frequencies on early diagenesis were studied in five different sedimentary environments, simulating the effects of a deeper-penetrating gear (e.g. a tickler chain beam trawl) versus a shallower, more variable penetrating gear (e.g. an electric pulse trawl). Trawling events strongly increased oxygen and nitrate concentrations in surface sediment layers and led to significantly lower amounts of ammonium (43 %–99 % reduction) and organic carbon in the top 10 cm of the sediment (62 %–96 % reduction). As a result, total mineralization rates in the sediment were decreased by up to 28 %. The effect on different mineralization processes differed both between sediment types and between trawling frequencies. The shallow-penetrating gear had a slightly smaller effect on benthic denitrification than the deeper-penetrating gear, but there were no statistically different results between gear types for all other parameters. Denitrification was reduced by 69 % in a fine sandy sediment, whereas nitrogen removal nearly doubled in a highly eutrophic mud. This suggests that even relatively low penetration depths from bottom fishing gears generate significant biogeochemical alterations. Physical organic carbon removal through trawl-induced resuspension of sediments, exacerbated by a removal of bioturbating macrofauna, was identified as the main cause of the changes in the mineralization process.

1 Introduction

Bottom trawl fisheries provide for 23 % of global fish landings (Cashion et al., 2018), with the vast majority of this type of fishing taking place in productive coastal shelf seas (Amoroso et al., 2018). In bottom trawl fisheries, nets are dragged along the bottom with the help of weighted devices such as otter boards, shoes, or beams, while chains, ground ropes, and/or electrical stimuli are used to coerce fish into the net. Bottom trawl gears penetrate the seafloor, up to 35 cm deep for otter trawl boards and 10 cm deep for tickler chain rigged beam trawls, depending on the gear specifics and the sediment type (Paschen et al., 2000; Lucchetti and Sala, 2012; Depestele et al., 2016). Hence, during a trawling event, sediment is mixed down to a certain depth, and hydraulic drag introduced by the moving gear can cause the erosion of an additional sediment layer (Depestele et al., 2016, 2019; O'Neill and Summerbell, 2011; O'Neill and Ivanović, 2016). Sediment disturbances by bottom trawling occur on very large scales: 63 % of all North Sea sediments are trawled between 1 and more than 10 times per year (Eigaard et al., 2017).

Scientific literature is rich in studies showing the physical and ecological alterations to the benthic environment caused by bottom trawling. Acute impacts of bottom trawling include the homogenization of surface sediment (Depestele et al., 2019; Ferguson et al., 2020) and the removal of significant proportions of benthic fauna (Bergman and Hup, 1992; Bergman and Van Santbrink, 2000; Tiano et al., 2020). Consistent fishing pressure favours organisms with shorter life spans and/or increasing resistance to trawling, while communities become depleted of species with key functional roles (Kaiser et al., 2006; Hiddink et al., 2017; Sciberras et al., 2018). Both fining (Trimmer et al., 2005) and coarsening (Palanques et al., 2014; Mengual et al., 2016) of the sediment have been attributed to trawling, as well as chronic organic matter depletion (Pusceddu et al., 2014; Paradis et al., 2019). The effects of bottom trawling on biogeochemical dynamics, however, remain relatively understudied. Trawling has been linked with enhanced carbon mineralization rates due to organic matter priming (van de Velde et al., 2018) and/or trawl-induced increases in organic material (Polymenakou et al., 2005; Pusceddu et al., 2005; Palanques et al., 2014; Sciberras et al., 2016). These results seemingly contrast with findings of organic matter depletion (Mayer et al., 1991; Brylinsky et al., 1994; Watling et al., 2001) and reduced mineralization rates after acute trawling (Tiano et al., 2019), highlighting the lack of knowledge on this topic and the need for further investigation.

Geochemical alterations impact the capacity of the sediment to recycle organic matter back to bioavailable nutrients (i.e. the sediment biogeochemistry). These are important processes in shallow coastal seas where primary production is strongly dependent on nutrients regenerated in the sediment (Soetaert and Middelburg, 2009; Provoost et al., 2013). Observed biogeochemical changes caused by sediment resuspension can lead to the instantaneous release of nutrients from the sediment into the water column (Durrieu de Madron et al., 2005), temporarily enhanced oxygen consumption (Tiano et al., 2019) and increased nutrient concentrations in the bottom water (Riemann and Hoffmann, 1991; Almroth et al., 2009; Couceiro et al., 2013). Furthermore, trawling has been linked to an increase in the sediment oxygenated layer depth (Allen and Clarke, 2007; Tiano et al., 2019) and a reduction of the denitrification capacity of cohesive sediments (Ferguson et al., 2020). It has been argued, based on in situ measurements, that the sediment biogeochemistry in consistently disturbed sediments remains in a transient state; i.e. the sediments are permanently recovering from a disturbance event (Van De Velde et al., 2018). These effects can potentially be mitigated with alternative fishing gears or modified gear configurations; however, the effectiveness of which needs to be assessed.

To reduce fishing impacts, alternative bottom trawl gears such as pulse fishing gears are being investigated (van Marlen et al., 2014; McConnaughey et al., 2020). With pulse gears, the heavy tickler chains are replaced by electrodes, which emit electrical pulses that induce a cramping response in flatfish (Soetaert et al., 2015b). This causes fish to become temporarily immobilized, allowing their capture in a net which drags behind the electrodes. Pulse gears exhibit lower penetration depths (∼ 50 %) than conventional beam trawls (Depestele et al., 2019) and also erode less material into suspension through hydrodynamic drag due to a reduced towing speed (Rijnsdorp et al., 2020a). The lower penetration depth of the pulse gear compared to standard tickler chain methods has been shown to decrease the effects of bottom trawling on the sediment redox layer (Depestele et al., 2019) and on chlorophyll a reduction (Tiano et al., 2019), but it is yet unclear how specific mineralization processes might be affected on longer temporal scales.

The aim of this study was to explore the possible impacts of bottom trawling on the sediment carbon and nitrogen cycling for two gears with different penetration depth distributions and with increasing trawling frequency. We use a dynamic diagenetic model, to which trawling disturbances were added. We parametrized the model for five locations in the North Sea, with sediments ranging from coarse sands to fine mud. Our hypotheses were (1) that the effects of bottom trawling would differ depending on the sedimentary environment and (2) that fishing gear with reduced sediment penetration would incur less changes in biogeochemical cycling.

2 Materials and methods

2.1 Model setup

2.1.1 Model description

To model the effects of bottom trawling on sediment biogeochemistry, disturbance events were added to a dynamic implementation of the early diagenesis model OMEXDIA (Soetaert et al., 1996a, b). This model describes the concentrations of organic matter, oxygen, nitrate, ammonium, dissolved inorganic carbon (DIC), and oxygen demand units (ODUs, reduced reaction products of anoxic mineralization). These are calculated on a 1D grid, with 100 layers increasing in thickness, starting from 0.01 cm at the sediment water interface (SWI) and extending up to a sediment depth of 100 cm. The incoming flux of organic matter (detritus) consists of a labile, fast-decaying fraction (FDET) and a semi-labile, slow-decaying fraction (SDET) and is mineralized in either oxic mineralization, denitrification, or anoxic mineralization (Table 1). With oxic mineralization and denitrification, the consumption of oxygen and nitrate as terminal acceptors is explicitly modelled (Table 1; Reactions R1, R2). Anoxic mineralization processes with alternative oxidants such as manganese oxides, iron oxides, sulfate, and organic matter are collected into one process that produces oxygen demand units (ODUs) as reaction products (Table 1; Reaction R3). ODU reoxidation and nitrification, the biological oxidation of ammonia to nitrate, are two additional processes that consume oxygen (Table 1; Reactions R4, R5). Mineralization rates are dependent on carbon availability (first-order kinetics) and oxidant availability (Michaelis–Menten type kinetics) and are inhibited by concentrations of inhibiting solutes (e.g. oxygen inhibits denitrification and anoxic mineralization). FDET, SDET, O₂, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, ODUs, and DIC are the seven state variables from which the concentrations are modelled in every layer.

Table 1 Diagenetic reactions used in OMEXDIA. x denotes the molar C : P ratio and y the molar N : P ratio in organic matter per mole of phosphorus (for Redfield stoichiometry, x=106, y=16).

Exchange of state variables between the different layers is caused by advection (sediment accretion, v) or molecular diffusion (for solutes), as well as bioturbation (for solids). The solute flux J_D due to molecular diffusion and advection is described by Fick's first law (Fick, 1855),

$$\begin{array}{}\text{(1)}& J_D=\mathit{\phi }{D}_i{\displaystyle \frac{\partial C}{\partial z}}+\mathit{\phi }vC,\end{array}$$

where the effective diffusion coefficient is estimated as $D_i=D_{\mathrm{0}}/{\mathit{\theta }}^{\mathrm{2}}$, with D₀ the molecular diffusivity of the solute; θ²=1–2ln(φ) the factor correcting for sediment tortuosity (Boudreau, 1996); φ the sediment porosity, which was kept constant with depth z; and C the concentration of the state variable. Molecular diffusion coefficients were calculated using the R package marelac (Soetaert and Petzoldt, 2018). Bioturbation is depth-dependent and assumes a constant biodiffusivity value Db₀ in a layer with thickness L_mix. Below this depth, bioturbation decreases rapidly to zero, determined by the attenuation coefficient for bioturbation (Db_coeff).

$$\begin{array}{}\text{(2)}& D{b}_z=D{b}_{\mathrm{0}}{e}^{-\left(z-L_{\mathrm{mix}}\right)}/D{b}_{\mathrm{coeff}}\end{array}$$

2.1.2 Model parametrization

The model was parametrized for five different sedimentary settings in the Southern North Sea (Fig. 1, Table 2): a coarse sandy sediment (hereafter denoted as “Coarse”) with a median grain size of 433 µm located on a sandbank in the Belgian part of the North Sea (BPNS); two sediments with intermediate grain sizes (216–220 µm) – one with a low nutrient load situated on the Dogger Bank (FineL) and the other nearshore (BPNS) with a comparatively high nutrient load (FineH); and finally two muddy sediments with a high silt content (74 %–88 %) – one with a low nutrient load situated offshore, on the Fladen Grounds (MudL), and one with a comparatively high nutrient load, situated nearshore in the BPNS (MudH).

Table 2 Characteristics of the selected sites. Low and high nutrient classification is based on relative differences in nutrient build-up for the same sediment type (see Fig. S1 in Supplement). MGS, SD.1, and SD.9 are median grain size and the boundaries of the 10th and 90th percentile of the grain size respectively (in µm). The percentages of sand (63–1000 µm) and mud (<63 µm) are weight percentages of a dried sediment sample sieved over a 1 mm sieve.

^* a: Toussaint et al. (2021); b: De Borger et al. (2021).

Figure 1 Sampling locations in the North Sea. Offshore stations MudL and FineL sampled by De Borger et al. (2021), and nearshore stations (inset) in the Belgian part of the North Sea sampled by Toussaint et al. (2021). Low (L) nutrient content relative to high (H) nutrient content. Basemap: © Esri, depth raster by GEBCO Compilation Group (2020).

Biogeochemical data from these stations were collected from box core samples (30 cm i.d., 25–30 cm sediment height) in two separate sampling campaigns on the North Sea, one in September 2017 (Coarse, FineH, MudH; Toussaint et al., 2021) and the other in May–June 2018 (FineL, MudL: De Borger et al., 2021). From these box core samples, sediment characteristics and the distribution of nutrients in the sediment were subsampled, and nutrient exchange rates across the sediment–water interface were determined in incubation experiments (see Sect. S1 in the Supplement for an extensive description of the methodology used to derive this information).

Model parameters included both measured concentrations in the bottom water, as well as process rate parameters that were derived following a two-step steady-state fitting procedure (Table 3). Using the measured DIC flux as the upper boundary organic carbon input flux, the O₂ flux and porewater profiles of O₂, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ were first fitted manually by tweaking a limited set of model parameters. The degradation rate of low degrading material (r_Slow) and the biodiffusivity constant (Db) were constrained by fitting NH${}_{\mathrm{4}}^{+}$ and O₂ profiles. Mechanistically, decreasing the bioturbation rate Db reduces the build-up of NH${}_{\mathrm{4}}^{+}$ with depth, increases the oxygen penetration depth, and changes the shape of the NO${}_{\mathrm{3}}^{-}$ profile (deepening the NO${}_{\mathrm{3}}^{-}$ peak). Whereas the degradation rate of the semi-labile organic matter also impacts deep NH${}_{\mathrm{4}}^{+}$ concentrations, it has a larger effect on the shape of the NH${}_{\mathrm{4}}^{+}$ profile, with lower degradation rates causing a more gradual build-up with sediment depth. Subsequently, parameters affecting the NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ profiles were tuned (the nitrification rate r_nit, as well as denitrification constants ksNO3denit and kinO2denit). Higher nitrification rates increase the build-up of NH${}_{\mathrm{4}}^{+}$ and increase concentrations of NO${}_{\mathrm{3}}^{-}$, typically producing a nitrate concentration peak within the oxic zone. The shape of the oxygen profiles further constrained the oxidation rate of oxygen demanding units (ODUs) and inhibition constants for anoxic mineralization (kinO2anox, kinNO3anox). Then followed an automated constrained parameter fitting step using an optimization algorithm. In this second step, the fitted parameters were allowed to vary in a range ± 10 % around the manually fitted parameter values. Also the DIC fluxes were refitted within a narrow range (0.98–1.02 of measured value) to allow freedom to the fitting algorithm. A random-based minimization algorithm (Price, 1977) implemented in the R package FME (Soetaert and Petzoldt, 2010) was used. This algorithm pseudo-randomly sampled the parameter space until the parameter set was found, which returned the minimal model cost, defined as the sum of variable costs (modelled – measured values), scaled using the mean–standard deviation relation determined for each nutrient.

Table 3 Parameters used for modelling the different sediment types.

Using the steady-state condition as the initial condition, several time-variable boundary conditions were imposed for the dynamic model simulations. A sinusoidally varying detrital carbon deposition flux, with the model-derived carbon flux (Cflux, Table 3) as the annual average and imposing an amplitude of 1, was used as the upper boundary organic carbon flux (Fig. 2a). The uniform amplitude of 1 for all sites was chosen to simplify temporal variations between sites. This resulted in differing organic carbon fluxes for each location. Additional time-variable boundary conditions (daily bottom water concentrations of O₂, NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$, as well as bottom water temperature) were extracted from the Copernicus Marine Environmental Monitoring Service implementation of the ERSEM model (European Regional Seas Ecosystem Model, Butenschön et al., 2016; Copernicus Marine Service Information, 2020) for each location.

Figure 2 (a) Example of sinusoidal varying organic carbon deposition flux, based on average DIC efflux measured by Toussaint et al. (2021) for the Coarse sediment. (b) Simulated depletion of bioturbation, relative to the maximum (y axis), with increasing trawling frequency (yr⁻¹, x axis) in fine sands. (c) Imposed mixing and erosion depths (cm), for fine sandy and muddy sediments, for the two gear types (tickler chain and pulse gear).

2.1.3 Disturbance modelling

Trawling disturbances were modelled as events causing the instantaneous removal of the surface layer due to hydraulic erosion (Depestele et al., 2016), followed by the mixing of a layer below that due to the actual gear penetration (Fig. 3). The hydraulic erosion was implemented as a reset of the sediment water interface (SWI) to the depth of the eroded layer. The mixing was implemented as a homogenization of solids (FDET, SDET) over the mixing depth (Fig. 3a), whereas solutes in the mixing depth (O₂, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, DIC) were set equal to the bottom water concentration of the respective solute to represent a complete flushing of the mixed layer with bottom water (Fig. 3b).

Figure 3 Implementation of a trawling event, on sediment concentrations of solids and liquids. Dotted grey lines: initial concentration profile; black line: profile after the event. (a) Hydraulic erosion removes a layer of sediment, moving the sediment–water interface (SWI) downward; the effect is implemented similarly for solids and liquids. (b) Subsequent mixing of the sediment homogenizes the solid concentration over the mixed layer depth, whereas liquids are set to the overlying bottom water concentration (BWconc). Depths of both impacts are defined in the text.

Modelled trawling events also caused an immediate reduction in bioturbation rates, due to the mortality of benthic fauna after a trawl pass (Fig. 2b). Benthic mortality is mostly dependent on the total penetration depth of the gear (Hiddink et al., 2017) but also varies with habitat (Pitcher et al., 2017). The instantaneous reduction in bioturbation was included as a proportional depletion (d), dependent on the sediment type (lowest in coarse sand, highest in mud) and the penetration depth (Fig. 2c). It was calculated based on the total gear penetration depth (TPD, i.e. the sum of the eroded layer depth and the penetration depth, in cm) and the mud content (% mud, particles <63 µm) of the sediment, as described in Eq. (3) derived from Sciberras et al. (2018).

$$\begin{array}{}\text{(3)}& d=\left(\mathrm{TPD}\cdot \mathrm{3}+\mathrm{mud}\cdot \mathrm{0.3}\right)/\mathrm{100}\end{array}$$

The subsequent recovery of the bioturbation was modelled based on the logistic growth equation (Eq. 4), with the recovery rate (r) the inverse of the longevity of the species community, kept constant at 0.04 yr⁻¹ (Rijnsdorp et al., 2016; Hiddink et al., 2019).

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}}=r\cdot N\left(\mathrm{1}-{\displaystyle \frac{N}{K}}\right)\end{array}$$

Here N is the bioturbation rate (cm⁻² d⁻¹), with K the full bioturbation rate for a species community at carrying capacity. The maximum reduction of the bioturbation was set to 90 % to account for quasi-immediate recolonization of trawled sediment by scavengers (Sciberras et al., 2018) and deeper living species that can survive intense trawling activity (Rijnsdorp et al., 2018). Electrical pulses were assumed to not affect benthos mortality in addition to the physical effects, since current available research shows very limited to no increased mortality by electrical pulses when compared to control situations (ICES, 2020; van Marlen et al., 2009; Murray et al., 2016; Soetaert et al., 2015a, 2016).

2.2 Simulations

The model description was implemented in R (R Core Team, 2020), the concentration changes of simulated species due to transport were calculated using the R package ReacTran (Soetaert and Meysman, 2012), and the resulting system of differential equations was solved using the deSolve package (Soetaert et al., 2010). Dynamic model simulations were initialized with a steady-state solution calculated with annually averaged boundary conditions as input parameters, using the R package rootSolve (Soetaert, 2009). This is necessary to build up an organic carbon, ammonium, ODU, and DIC stock in the sediment (Soetaert et al., 1996b). Dynamic simulations were run for 15 years, with daily output, to generate sufficient independence from starting conditions. Reported modelling results stem from the last simulated year.

The frequency of the trawling events imposed ranged from 0 (the baseline) to 5 yr⁻¹, based on realistic values of bottom trawling intensities in the North Sea (Rijnsdorp, 1998; Eigaard et al., 2017). Events were distributed randomly throughout the year, given the absence of a clear seasonal pattern in trawling intensities (Rijnsdorp et al., 2020a). For each trawling frequency and site, 30 model simulations were performed, each with a different pair of penetration and erosion depths, generated from a log-normal distribution of penetration depths. For the tickler chain trawl, a log-normal distribution of penetration depths was generated given the average values (95 % confidence limits) for sand and mud summarized by Pitcher et al. (2017): 3.2 (1.5, 6.7) cm and 1.9 (1.0, 3.7) cm for mud and sand respectively; the erosion depth was set to 22% of the penetration depth. The penetration depths for the pulse gear were set to 50 % of the tickler gear: 1.6 (0.75, 3.38) cm and 0.95 (0.49, 1.83) cm for mud and sand. The erosion depth for the pulse gear was set to 70 % of the tickler gear (Fig. 2c; Depestele et al., 2016, 2019; Rijnsdorp et al., 2020b).

2.3 Statistical analysis

Linear models were constructed with selected model output variables to analyse the effects of the different gear types, trawling frequency, and the sedimentary context on the rates of the different mineralization processes and on the total mineralization (the sum of the separate mineralization processes). A normal distribution was adopted for the process responses. To deal with heterogeneity of variances of the residuals (for all models) in the linear models, a generalized least squares (GLS) structure was added (Pinheiro and Bates, 2000; Zuur et al., 2009; West et al., 2014), which allowed for unequal variances among treatment combinations to be included as a variance structure (Pinheiro and Bates, 2000; West et al., 2014). To find the most suitable variance structure, models with different variance structures were compared using Akaike information criterion (AIC) scores (Akaike, 1974) and plots of fitted values and individual model terms versus the residuals (Zuur et al., 2009). For all models, a variance structure was selected that allowed for variances conditional on the station and trawling frequency. This variance structure was of the form sediment × frequency, using the varIdent function of the R package nlme (Pinheiro et al., 2019). Subsequently the fixed model component was optimized by manual stepwise selection, using the likelihood ratio test and associated p values as validation for removing excess terms (Zuur et al., 2009). During this step, the philosophy was adopted to not include significant interaction terms containing a certain variable when said variable was not significant by itself. The minimal adequate model was represented using restricted maximum likelihood estimation (REML, West et al., 2014). GLS models were implemented using the R package nlme (Pinheiro et al., 2019).

3 Results

3.1 Baseline model simulations

Baseline model simulations (undisturbed) show the differences in organic matter cycling between the chosen locations. In the coarse-sand station (Coarse), the average total mineralization rate was 13.6 mmol C m⁻² d⁻¹, with 89 % of this due to oxic mineralization and 6 % due to anoxic mineralization, and 5 % was denitrified (see Table S2, Supplement). The two muddy stations had either very high or very low total mineralization rates (MudH: 82 mmol C m⁻² d⁻¹, MudL: 8.5 mmol C m⁻² d⁻¹) and similar for the two fine sandy stations (FineH: 30 mmol C m⁻² d⁻¹, FineL: 8.1 mmol C m⁻² d⁻¹), related to the difference in organic matter deposition between nearshore and offshore locations (Table 3). Oxic mineralization dominated in FineL (oxic: 81 %, anoxic: 12 %, denitrification: 6 %), FineH (oxic: 78 %, anoxic: 22 %, denitrification: 0 %), and MudL (oxic: 72 %, anoxic: 18 %, denitrification: 10 %), whereas the mineralization in the nearshore muddy station (MudH) was dominated by anoxic processes (oxic: 27 %, anoxic: 68 %, denitrification: 5 %).

3.2 Impact on biota

Trawling-induced depletion of fauna substantially decreased average annual bioturbation rates. Bioturbation decreased with increasing penetration depth (Fig. 2b, c), resulting in the strongest decreases in muddy sediment (MudL, MudH) and larger decreases in the deep-penetrating gear versus the shallow-penetrating gear (Table 4). In the Coarse sediment, the annually averaged bioturbation decreased gradually, from 81 % of its original value at one trawl per year to 19 % at five trawls per year in the tickler gear and 49 % at five trawls per year for the pulse gear. For the fine (FineL, FineH) and muddy (MudL, MudH) sediments the maximum depletion was reached after four (five for the shallow gear) and two trawling events respectively.

Table 4 Percentage of bioturbation (average ± SD) remaining after sustained trawling activity at a given trawling frequency (yr⁻¹), for the different sediment types. T: deeply penetrating tickler gear, P: shallow-penetrating pulse gear.

3.3 Nutrient and organic carbon distribution

With higher trawling frequencies, concentrations of oxygen and nitrate in the sediment generally increased, whereas ammonium and organic carbon contents were always reduced (Figs. 4, 5). The magnitude of concentration changes was similar for both gear types (Table S1, Supplement shows the mean percentage change relative to baseline concentrations reported). Increases in oxygen (Fig. 4a–e) and nitrate (Fig. 4f–j) concentrations were largest in the oligotrophic stations FineL and MudL, where concentrations of oxygen in the upper 5 cm increased 15–16 fold (respectively 1604 % and 1516 %), while nitrate concentrations increased 9–19 fold (respectively 909 % and 1911 %) at the highest trawling intensities (see Fig. S2 in Supplement for the range of nutrient concentrations throughout the year). In contrast, O₂ and NO${}_{\mathrm{3}}^{-}$ concentrations initially decreased at MudH by 25 % and −50 % maximally for one to two trawls per year, before increasing by 52 % to 81 % (O₂) and then increased by 123 % to 188 % (NO${}_{\mathrm{3}}^{-}$) at five trawls per year. Ammonium (NH${}_{\mathrm{4}}^{+}$) concentrations decreased strongly in all sediments, with a decrease of up to 69 % in Coarse; 68 % in MudH; and >90 % in FineL, FineH, and MudL (Fig. 4k–o, Table S1).

Figure 4 Annually averaged modelled concentrations (mmol m⁻³) of oxygen (a–e), nitrate (f–j), and ammonium (k–o) in the different types of sediment used as the basis for the disturbance simulations (columns). The black dotted line is the 0 trawl default, and the full and dotted coloured lines are tickler and pulse gear respectively, with increasing trawling frequencies as different colours.

Figure 5 Total organic carbon concentrations in the upper 10 cm of the sediment (mmol m⁻³) for the different stations (x axis on log scale). Profiles represent the average profile per trawling intensity (n=30) for the deep-penetrating gear (full) and shallow-penetrating gear (dotted line). For station MudH a cutout of the first 0.5 cm is inflated to show the displacement of OC to the surface.

Increasing the trawling frequency reduced the total amount of reactive organic carbon (labile + semi-labile, OC) in the sediment and reduced the penetration depth of the OC (Fig. 5, Table S1). Trawling frequencies of 3–5 yr⁻¹ led to near-total depletion of reactive OC in all sediments in the upper 10 cm (>95 % removed for Coarse and FineH and >90 % for FineL, MudL, and MudH, Fig. 5). The mean OC profiles for both gears at a given frequency were often visually different (dotted vs. full lines on Fig. 5), but the average concentrations over 10 cm did not differ significantly. A redistribution of organic carbon was visible in the upper centimetre of the sediment, where organic carbon concentrations were higher in the impacted than in the baseline simulation (example in the cutout of the top 5 mm shown for MudH, Fig. 5). In FineL, MudL, and MudH the ratio of labile organic carbon (FDET) as a proportion of the carbon pool (FDET + SDET) increased between 25 % and 34 % (Fig. S3). This effect was only noticeable in the upper 0.2–0.5 cm; below this, depth values of this ratio in all trawling frequencies converged to 0 due to the depletion of labile organic carbon.

3.4 Total mineralization rates

The trawling frequency had a significant negative impact on all studied mineralization process rates (oxic, anoxic, denitrification) and on the total organic carbon mineralization, as confirmed by the negative coefficients in the GLS models (“Freq”, Table 5). Changes in process rates also differed between the studied sediments, as seen by the inclusion of an interaction term between the sediment type and the trawling frequency (Freq : Sed). The sediment biogeochemical response to increasing trawling frequency was often non-linear, warranting the inclusion of a squared frequency term (Freq²). The gear type was only included as a significant explanatory variable in the model for denitrification, where the deeper-penetrating gear (tickler) decreased denitrification rates more than the shallow-penetrating gear (Table 5).

Table 5 Generalized least squares (GLS) models for the total mineralization, oxic mineralization, anoxic mineralization, and denitrification as a function of increasing trawling frequency (Freq and Freq²), the fishing gear type (Gear), and the sediment context (Sed) and interactions between model terms Freq : Sed and Freq² : Sed.

^* S1: MudH; ${}^{**}$ S2: FineL; ${}^{***}$ S3: Coarse; ${}^{****}$ S4: MudL.

The total mineralization rate was impacted negatively in all cases, and decreases ranged from 5 % for one trawl per year for MudL to −28.9 % for five trawls per year, for FineL (Fig. 6a–e, values summarized in Supplement Table S2).

Figure 6 Rates of total mineralization (a–e), and the three main mineralization processes (f–j: oxic; k–o: anoxic; p–t: denitrification) (y axis, mmol m⁻² d⁻¹) for each gear type (blue boxes: tickler gear; red boxes: pulse gear) and for increasing trawling frequency (x axis, yr⁻¹).

The change in oxic mineralization rates (base: 12.0, 23.3, 6.5, 6.1, and 23.1 mmol m⁻² d⁻¹ for Coarse, FineL, FineH; MudL, and MudH respectively) showed different patterns depending on the station (Fig. 6f–j). For Coarse and FineH there was a consistent decrease in oxic mineralization rates with increasing trawling frequency, with maximum decreases at five trawls per year of 21 % and 23 % for the tickler gear and 21 % and 25 % for the pulse gear. In contrast, for FineL and MudL oxic mineralization remained the same or increased at a frequency of one trawl per year (8 % and 11 %), followed by a decrease at higher trawling frequencies with maximal decreases of 15 % and 10 % at five trawls per year for the tickler gear; the values for the pulse gear were similar (Fig. 6f–j). For MudH, oxic mineralization increased by a maximum of 56 % at five trawls per year for both gears.

Anoxic mineralization rates (base: 0.9, 1.0, 6.6, 1.5, and 57.7 mmol m⁻² d⁻¹ for Coarse, FineL, FineH; MudL, and MudH respectively) were also affected similarly by the two gear types and decreased for all stations, though with differing magnitudes (Fig. 6k–l). The lowest decrease was in the Coarse sediment, where the decrease in the anoxic mineralization rate was similar for all trawling frequencies (range of 18 % to 25 %), and the highest decrease was modelled at MudL, where anoxic rates decreased between 73 % and 83 %.

Denitrification rates (base: 0.6, 0.5, 0.1, 0.8, and 4.2 mmol m⁻² d⁻¹ for Coarse, FineL, FineH, MudL, and MudH respectively) decreased with increasing trawling frequencies at all stations except MudH, with a maximum reduction of 74 % (tickler) and 68 % (pulse) at FineL (Fig. 6p–t). Trawling frequencies of 1–2 yr⁻¹ did not consistently alter denitrification for the Coarse sediment and FineH and MudL. For the MudH trawls, denitrification rates increased by 50 % and 49 % towards five trawls per year for tickler and pulse gear respectively.

As a result of the changes to denitrification, the removal of reactive N from the sediment changed. The sediments where denitrification decreased most (FineL, MudL) had 35 % and 51 % of N produced by mineralization removed by denitrification when undisturbed, and this reduced to 11 % and 45 % respectively for five trawls per year. For the Coarse sediment the fraction of N removed increased from 26 % when undisturbed to 30 % for one trawl per year and then decreased again to 25 %. In MudH more N was removed as well, with a near doubling as a peak at five trawls per year (48 %, up from 26 % as the base).

All previous results represent average changes throughout the year, but trawling also showed instantaneous effects, as illustrated by the decrease in denitrification rates (to nearly 0 mmol m⁻² d⁻¹) immediately after a trawl events (Fig. 7).

Figure 7 Nitrate concentrations (mmol m⁻³) in the sediment for a no-trawling simulation (a, d, g) and for three trawls per year (d, e, f), and associated annual denitrification rates (mmol m⁻² d⁻¹, black lines indicate default (zero trawls per year) and red dotted lines indicate the trawling simulation). Represented stations are Coarse (a–c), FineL (d–f), and MudH (g–i). Black vertical lines indicate when trawling events take place. ”Doy” indicates day of year.

3.5 Relative changes

The relative contribution of the mineralization processes to the total mineralization changed markedly between trawling frequencies and stations (Fig. 8). In general, the proportion of oxic mineralization increased (Fig. 8a–e) at the expense of anoxic mineralization (Fig. 8f–j). The largest changes occurred when switching from zero to one trawling event yr⁻¹, and values remained stable from two events per year onwards. The proportion of oxic mineralization increased most at MudH (tickler: 116 %, pulse: 112 % for five trawls per year), and the smallest changes occurred at Coarse (<1 % for both gears). The proportion of anoxic mineralization, on the other hand, decreased in all simulations. The largest changes were modelled at FineL and MudL (69 % and 78 % respectively) and the smallest for Coarse and FineH (14 % and 18 % respectively). The proportion of mineralization performed by denitrification decreased in FineL and MudL (71 % and 10 %, five events per year), doubled at MudH (100 %), and remained practically the same for Coarse and FineH (Fig. 8k–o).

Figure 8 Relative contributions (y axis, as proportion) of the three main mineralization processes (a–e: oxic; f–j: anoxic; f–o: denitrification) to the total mineralization for each gear type (blue boxes: tickler gear; red boxes: pulse gear) and for increasing trawling frequency (x axis, yr⁻¹).

4 Discussion

4.1 Organic carbon depletion

Simulated trawling of the seafloor impacted the sediment biogeochemistry in all environments and for all trawling frequencies. The amount of total mineralizable carbon in the sediment consistently decreased with higher trawl frequencies, but the changes in mineralization pathways differed from case to case. The main drivers of the biogeochemical changes were found to be the depletion of organic carbon (OC) in the sediment (i.e. the substrate for mineralization itself), the redistribution of this OC nearer to the SWI (Fig. 5), and the increasing oxygenation of the sediment. With each trawl pass, a part of the organic-carbon-rich top layer is removed (Fig. 5). This is associated with an injection of oxidized reactants from the water column (O₂, NO${}_{\mathrm{3}}^{-}$) deeper into the sediment and a homogenization of OC concentrations in the mixed layer during the mixing phase. Simultaneously, part of the benthos in the sediment is removed, often strongly decreasing the bioturbation rate, affecting the rate at which organic matter is distributed in the sediment (especially after multiple trawling events yr⁻¹, Table 4). Sediment mixing alone could potentially increase OC contents at the bottom of the mixing zone, but successive trawling events, and the removal of bioturbators that can transport OC far below the mixing zone, resulted in a redistribution of OC closer to the sediment–water interface in all simulated sediments.

The fine sandy station with low organic matter content (FineL), as well as both muddy stations (MudL, MudH), showed smaller decreases in surface organic carbon concentrations compared to the eutrophic fine sandy (FineH) and Coarse sediments (Table S1), mainly because baseline bioturbation rates in the former were 3 orders of magnitude larger (Table 3). As such, bioturbation seems to cause an increased resistance to carbon loss by facilitating transport to deeper layers, making it less vulnerable to surface disturbances. While the physical OC depletion caused by the penetrating gear is aggravated by the loss of bioturbating fauna in the sediment, this effect is context-dependent as bioturbators show variable levels of resistance to trawling (Hale et al., 2017; Tiano et al., 2020). Our modelled results provide further evidence that surviving fauna help buffer and mitigate the biogeochemical effects of trawling (Duplisea et al., 2001). The comparison with the work of Duplisea et al. (2001) is in fact remarkable, as Duplisea et al. (2001) used a food-web-based model to investigate changes to carbon cycling, whereas a diagenetic model was used in this work. Both are very different approaches, which highlight both a shift to more oxic mineralization and the importance of benthic fauna as a stabilizing factor. Tiano et al. (2019) observed decreases in sedimentary Chl a in the upper 1 cm immediately after trawling of 41 % and 83 % for pulse and tickler gears respectively. Also, organic matter (OM) depletion as a result of long-term fishing has been reported, even at water depths beyond 500 m (Martín et al., 2014; Pusceddu et al., 2014; Paradis et al., 2019), where comparisons between trawled and untrawled sites yielded a difference in OC between 20 % and 60 % (Paradis et al., 2019), or 60 %–100 % of the daily input flux of organic carbon was removed from sediments by trawling (Pusceddu et al., 2014). In fact, these deep-water sediments are particularly sensitive to trawling disturbances, a concerning feature given the steady expansion of fishing practices into deeper waters in recent decades (Morato et al., 2006; Puig et al., 2012; Watson and Morato, 2013). Deep-water species communities are slow growing and thus recover slowly, organic matter deposition rates are low, and the generally finer grained sediments found in the deep are easily resuspended following a trawl passage (Norse et al., 2012; Mengual et al., 2016). All three of these factors increase the impacts of trawling events on organic matter cycling in the model presented in this work, and further modelling work could be useful to investigate the potentially large impact that deep-sea habitat experience. There are also studies reporting enhanced OC concentrations in trawled areas, in contrast with our results (Palanques et al., 2014; Pusceddu et al., 2005; Sciberras et al., 2016), possibly due to differing hydrodynamic and morphological conditions of the North Sea compared to other areas.

The decrease in total mineralization rates may partly be offset by redeposition of organic matter, which was not considered in our model. Not all eroded organic matter stays in the water column, but a part resettles on the sediments. How this redistribution occurs depends on the sediment type and the local hydrodynamics, which determine the distance over which eroded sediment particles are transported (Le Bot et al., 2010; Robinson et al., 2005). It can be expected that for coarser, heavier sediments a fraction will be redeposited in the trawling track but that for muddy sediments, lighter and rich in organic matter (Mayer, 1994), eroded material remains in suspension long enough to be transported elsewhere (Palanques et al., 2014). In the North Sea, suspended material is transported from the Southern Bight northward by anticlockwise residual currents. Ultimately, the partially degraded materials are deposited in the Skagerrak (Dauwe et al., 1998). So given the intensity, and the persistence with which vast areas of the southern North Sea are trawled (e.g. total annual sediment mobilization by the Dutch trawling fleet varied between 8 and 17×10¹⁴ kg of sediment between 2010 and 2016, Rijnsdorp et al., 2020a), we expect that trawling-induced sediment resuspension plays a significant role in the northward transport and actively contributes to organic matter depletion in southern areas.

4.2 Changes to mineralization pathways

Trawling activities generally caused strong increases in sedimentary oxygen and nitrate availability and decreases in the ammonium content (Figs. 4, 7, Table 3). Sediment oxygenation increased both because of a direct injection of oxygen-rich bottom water in deeper sediment layers during a trawling event and because oxygen consumption by mineralization processes decreased as a result of strong decreases in OC and ammonium. As a result of increased oxygen availability, the importance of oxic mineralization generally increased with trawling, whereas anoxic mineralization decreased (Fig. 8).

The strongest increases in the proportion of oxic mineralization were modelled for the sediment characterized by a high silt percentage and organic matter load (MudH). These types of sediments also have a low permeability, high mineralization rates, and a low oxygen penetration depth (Braeckman et al., 2014), with a lesser importance of oxic relative to anoxic mineralization in undisturbed conditions. Fishing gears penetrate deepest in these muds and as such provide oxygen to deeper layers, although this is consumed rapidly. Van De Velde et al. (2018) found an increase in mineralization rates of over 200 % after a disturbance event in muddy sediments (from the same origin as MudH). This was attributed to multiple possible factors, such as self-priming by mixing refractory with labile organic matter, burial of phytoplankton in settling sediment, and the introduction of oxygen to redox shuttle mechanisms. Our results show an enhancement of total oxic mineralization at MudH but no increase in the total mineralization rate, perhaps because the aforementioned processes were not included in the used model.

The higher oxygen concentrations also had a clear inhibiting effect on denitrification rates (Fig. 7). Denitrification in coastal shelf sediments accounts for an estimated third of all nitrogen loss in Earth's marine surfaces (Middelburg et al., 1996), making these regions crucial to counteract nitrogen eutrophication (Galloway et al., 2004; Seitzinger et al., 2006). For all stations, trawling events caused an instantaneous dip in denitrification rates, because of the injection of O₂, and removal of the electron acceptor (NO${}_{\mathrm{3}}^{-}$) from the sediment, similar to the results of van der Molen et al. (2013) (Fig. 7, black vertical lines). However, a discrepancy was noted between the biogeochemical impacts of trawling in cohesive sediments with high organic matter concentrations versus sandier and comparatively low nutrient sediments, consistent with literature findings (Polymenakou et al., 2005; Van De Velde et al., 2018; Tiano et al., 2019). The coarsest sandy sediments (Coarse) were by default deeply oxygenated Fig. 4a–b), with denitrification maximally inhibited by oxygen concentrations. Larger pore spaces in these sediments allow for bottom water to penetrate more deeply into the sediment matrix, bringing oxygen and other reactants into deeper sediment layers (Huettel and Gust, 1992). Cohesive sediments mostly lack such advective transport. As a result, oxygen often fuels rapid mineralization in coarse-grained sediments (Huettel and Rusch, 2000; Ehrenhauss et al., 2004). The increasing trawling frequency in coarse sediments thus had little effect on oxygen penetration, and nitrate concentrations only marginally increased (Fig. 7a–b), resulting in minor changes to mineralization pathways such as denitrification on average, although instantaneous effects could be prominent (Fig. 7c). In oligotrophic finer sediments (FineL, MudL), there was a massive increase in both oxygen and nitrate concentrations as a result of trawling (Fig. 7d–e). Whereas increasing NO${}_{\mathrm{3}}^{-}$ concentrations would stimulate denitrification on their own, the rise of oxygen concentrations strongly inhibited denitrification, leading to a drop in denitrification rates throughout the year (Fig. 7f). The more eutrophic fine sediment FineH displayed a similar pattern of increased oxygenation mineralization as the other fine sandy sites, but the already low baseline denitrification rates did not decrease further. In the unperturbed simulation of MudH, mineralization was predominantly anoxic (68 %), with denitrification limited by NO${}_{\mathrm{3}}^{-}$ availability. Increasing the oxygenation in this type of sediment caused the denitrification to double, by increasing the nitrate availability (Fig. 7g–i).

Ferguson et al. (2020) found that denitrification rates in Moreton Bay, Australia, were reduced between 11 % and 50 % within 3 h after a trawling event, and this rate decreased after successive trawling events during the studied period. These decreases were attributed to homogenization of the sediment, which removes oxic microniches created by fauna and thus zones of intense coupled nitrification–denitrification (Ferguson et al., 2020). Though the key role of redox microniches is not directly investigated here, we acquired decreases in denitrification rates in a similar range in all sediments apart from MudH, especially when trawling frequencies increased.

Within the marine environment, shelf sediments are sites characterized by high nutrient concentrations, therefore offering resilience against reductions in nutrient loadings. Soetaert and Middelburg (2009) showed that storage of ammonium in sediments significantly delays the response of shallow systems to oligotrophication, as the efflux of nitrogen from the sediment will compensate for part of the losses in the water column. The increased reduction of the ammonium concentrations with fishing intensity will affect this buffering capacity of the sediment. The nitrogen buffering capacity of the investigated sediments, representative for a large fraction of North Sea sediments, was affected similarly by both gears. Firstly, the stock of nitrogen (as NH${}_{\mathrm{4}}^{+}$) in the sediment was directly affected by porewater flushing during trawling (with decreases >99 % in some cases, Fig. 7f). Secondly, lower availability of reaction substrate (OC, NH${}_{\mathrm{4}}^{+}$) decreased denitrification rates, reducing N removal to the atmosphere.

4.3 Reducing gear disturbance of the seafloor

In our work the differences with respect to organic matter mineralization dynamics between gear types with differing penetration depths (e.g. 3.2 ± 1.2 cm vs. 1.6 ± 0.6 cm in mud) mostly remained suggestive rather than statistically conclusive. The fishing gear type was only included as a significant predictor for the denitrification rates, with a small coefficient (−0.00038; Table 5). This is because the freshly deposited stocks of organic carbon are present near the sediment surface, and any gear that penetrates the sediment impacts this layer, especially for multiple trawling events per year (Fig. 5, Table 5). This indicates that only a thin layer of surface sediment needs to be disturbed to generate significant biogeochemical changes (Dounas et al., 2005). Many biogeochemical processes are mediated by the dynamics of oxygen near the sediment–water interface, which itself is influenced by the composition and permeability of the sediment. A shift towards fining (an increased proportion of finer grain size classes) has been described in certain trawled areas, with expected consequences for sediment biogeochemistry, such as an increased rate of sulfate reduction (Trimmer et al., 2005). But the opposite occurs just as well (Depestele et al., 2019; Mengual et al., 2019; Tiano et al., 2019). In these cases resuspended fine-grained material is exported away from the trawling site, leaving a coarsened trawling track, with the results subtly different between gear types (Depestele et al., 2019; Tiano et al., 2019). This means that the effects of fishing gears on grain size sorting should be better characterized for various sediment types to constrain the uncertainty around predictions of gear impacts on sediment functioning.

This does not imply that the penetration depth is irrelevant. Other studies have reported clear positive effects of reducing the penetrations depths of fishing gears, such as decreased sediment mobilization and homogenization, and reduced organic matter depletion (Depestele et al., 2019; Tiano et al., 2019). Conversely, assuming that the eroded layer and the mixing depth scale with more deeply penetrating gears than those tested here, the depletion effects should become more pronounced. In a single trawling event more of the (reactive-carbon-rich) top sediment layer would be removed, there would be a higher mortality of organisms, and more of the nutrient build-up would be removed.

Aside from the penetration depth, the largest impacts occurred when increasing the trawling frequency from zero to one trawling event per year, and the response of mineralization processes to increased trawling frequency was often non-linear, making them more difficult to predict. This would imply that management strategies aimed at maintaining the ecosystem functions provided by shelf sediments should be focused on spatial controls, bottom impact quotas, and effort control of trawling gears that per definition require contact with the bottom to catch commercially viable target species (McConnaughey et al., 2020). An effective strategy limits the impacted surface area and allows carbon stocks and faunal communities in the sediment to recover from a disturbance, resulting in the recovery of vital biogeochemical functions such as denitrification and carbon burial. This includes technical adaptations to improve catch efficiency. Whereas our study only focused on direct head-to-head comparisons between the two gear types, pulse trawls are associated with lower spatial footprints due to their relatively higher catch efficiencies compared to beam trawls (Poos et al., 2020; Rijnsdorp et al., 2020a, b; ICES, 2020).

Shifting the fishing effort from peripheral areas to core fishing grounds would also reduce the area where the top sediment layer is removed on a regular basis, along with the associated reduction in mineralization of organic matter. This can be achieved through time-restricted bottom trawling, in which fishing grounds are closed off temporarily. Similarly, areas with high denitrification rates, crucial for eutrophication mitigation, can be closed to trawling completely, as done in other regions of the world (Ferguson et al., 2020). Site-specific conditions such as rates of biogeochemical recovery and sedimentation rates need to be known to determine the resilience of ecosystems to trawling and to fine-tune management plans (Paradis et al., 2021).

5 Conclusions

With the addition of perturbation events to a model of early diagenesis, and a description of faunal mortality and recovery, we simulated the effects of increasing bottom trawling frequencies on sediment biogeochemistry. The results showed that bottom trawl fisheries strongly impacted the sediment biogeochemistry, and the magnitudes of the changes were dependent on the sedimentary context and trawling frequency. Two types of fishing gears were investigated. The exposed top sediment layer rich in organic matter was targeted similarly by both fishing gears, resulting in a similar loss of organic carbon, which was further exacerbated by the loss of bioturbating fauna. A shift towards increasingly oxic mineralization at the cost of anoxic mineralization was observed, driven by an often strongly increased oxygen availability in the sediment. The removal of fixed nitrogen by denitrification was not affected similarly in all sediments. Denitrification increased in nearshore cohesive mud and decreased elsewhere, with highest decreases in offshore sediments with lower carbon loads. Our modelling results corroborate multiple patterns found in other studies and can serve to interpret research and search for mitigation strategies. Trawling impacts are hard to mitigate by only reducing the penetration depth of the gear, so additional management strategies are needed to allow for partial or full recovery of biogeochemical functions in between trawling events.