Offshore Windfarm Footprint of Sediment Organic Matter Mineralization Processes

Offshore windfarms (OWFs) offer part of the solution for the energy transition which is urgently needed to mitigate effects of climate change. Marine life has rapidly exploited the new habitat offered by windfarm structures, resulting in increased opportunities for filter- and suspension feeding organisms. In this study, we investigated the effects of organic matter (OM) deposition in the form of fecal pellets expelled by filtering epifauna in OWFs, on mineralization processes in the sediment. OM deposition fluxes produced in a 3D hydrodynamic model of the Southern Bight of the North Sea were used as input in a model of early diagenesis. Two scenarios of OWF development in the Belgian Part of the North Sea (BPNS) and its surrounding waters were calculated and compared to a no-OWF baseline simulation. The first including constructed OWFs as of 2021, the second containing additional planned OWFs by 2026. Our results show increased total mineralization rates within OWFs (27–30%) in correspondence with increased deposition of reactive organic carbon (OC) encapsulated in the OM. This leads to a buildup of OC in the upper sediment layers (increase by ∼10%) and an increase of anoxic mineralization processes. Similarly, denitrification rates within the OWFs increased, depending on the scenario, by 2–3%. Effects were not limited to the OWF itself: clear changes were noticed in sediments outside of the OWFs, which were mostly opposite to the “within-OWF” effects. This contrast generated relatively small changes when averaging values over the full modeling domain, however, certain changes, such as for example the increased storage of OC in sediments, may be of significant value for national / regional carbon management inventories. Our results add to expectations of ecosystem-wide effects of windfarms in the marine environments, which need to be researched further given the rapid rate of expansion of OWFs.

Offshore windfarms (OWFs) offer part of the solution for the energy transition which is urgently needed to mitigate effects of climate change. Marine life has rapidly exploited the new habitat offered by windfarm structures, resulting in increased opportunities for filter-and suspension feeding organisms. In this study, we investigated the effects of organic matter (OM) deposition in the form of fecal pellets expelled by filtering epifauna in OWFs, on mineralization processes in the sediment. OM deposition fluxes produced in a 3D hydrodynamic model of the Southern Bight of the North Sea were used as input in a model of early diagenesis. Two scenarios of OWF development in the Belgian Part of the North Sea (BPNS) and its surrounding waters were calculated and compared to a no-OWF baseline simulation. The first including constructed OWFs as of 2021, the second containing additional planned OWFs by 2026. Our results show increased total mineralization rates within OWFs (27-30%) in correspondence with increased deposition of reactive organic carbon (OC) encapsulated in the OM. This leads to a buildup of OC in

INTRODUCTION
Rising concern for global climate change has increased the urgency to lower carbon emissions from individual and industrial energy consumers (United Nations, 2015). The transition from fossil fuel driven energy production toward renewable energy sources from solar, hydrodynamic, and wind sources, is a booming industry in alignment with the goal of achieving "carbon-neutral" anthropogenic development (European Commission, 2020). Offshore windfarms (OWFs), groupings of wind turbines on submerged sediments, offer part of the solution in the energy transition. Compared to terrestrial windfarms, the advantages of a flat, wind-swept marine area, and the possibility to build larger structures, overrule increased construction and maintenance costs (Inger et al., 2009;Bergström et al., 2014).
With respect to marine life, the presence of the turbine substructure and the scour protection layer (a layer of rocks around the foundation base that prevents erosion) on the seabed represents the addition of a new type of habitat, which spans the full water column (Krone et al., 2013). Many organisms, mostly filter-and suspension feeders, thrive on this new habitat. For wind turbine substructures in the North Sea, a clear vertical succession of fouling fauna can be observed, from the barnacle Semibalanus balanoides in the intertidal zone, to the blue mussel Mytilus edulis in the upper subtidal, the amphipod Jassa herdmani below that, and the anemone Metridium senile in the lower regions of the foundation (De Mesel et al., 2015;Mavraki et al., 2020). The attraction of organisms extends beyond the surface of the turbine foundations. Fish species such as cod (Gadus morhua) and pouting (Trisopterus lusculus) are known to aggregate near turbine foundations (Reubens et al., 2011;Langhamer, 2012), and benthic communities of macrofauna have been observed to change, alongside a fining of the sediment and an enrichment with organic matter (OM) (Bergström et al., 2012;Coates et al., 2014;Leewis et al., 2018). Deposition of fecal pellets by the fouling fauna, as well as biomass falling from the structures is a likely source of this enrichment (Krone et al., 2013;Lefaible et al., 2019).
For sediment biogeochemistry, this increased carbon deposition may have far-reaching effects. In sediments, OM is mineralized to free inorganic nutrients available for water column processes by the respiration activities of microbes and fauna (Soetaert and Middelburg, 2009;Provoost et al., 2013). Sediments in shelf seas represent a considerable sink of carbon, either at the short-term by means of temporary storage before mineralization takes place, or over longer timescales as partially degraded carbon is buried in deep sediment layers where its reactivity decreases (Legge et al., 2020). Moreover, shelf sea sediments are of vast importance for the removal of excess nitrogen from the ecosystem through denitrification, the mineralization of OM through nitrate reduction, which produces nitrous oxide and dinitrogen gas (Middelburg et al., 1996;Galloway et al., 2004;Seitzinger et al., 2006).
Recently, Slavik et al. (2019) used a model that included the filtration by blue mussels attached to wind turbine foundations to show how this process can decrease primary production up to 8% in the OWF footprint, but also have effects noticeable 50 km away from these concentrations of filter feeders. Similarly, by upscaling carbon tracer experiments, Mavraki et al. (2020) calculated an uptake of 1.3% of the primary production standing stock by M. edulis and J. herdmani in the Belgian Part of the North Sea (BPNS), an area characterized by rapid OWF expansion. This highly localized removal of primary produced carbon will likely cause a significant transfer of reactive organic carbon (OC) in the form of fecal pellets and dislodged biomass from the water column to sediments in the vicinity of OWFs.
Observational evidence of mineralization processes in OWF sediments, let alone of changes thereof, is rare (Toussaint et al., 2021) due to the difficulties of sampling sediments there. However, consequences of organic enrichment on sediment biogeochemistry are well described for example in fish and mussel aquaculture. Here, increased deposition of organic material (i.e., feces, fish feed, dead biomass) causes an increase in mineralization activity in the sediment, which increases sedimentary oxygen consumption, and an increase in anoxic mineralization processes. This leads to higher CO 2 release from sediments, increased exchange of nutrients, and a higher total OC content in the sediment (Nizzoli et al., 2005;Kalantzi and Karakassis, 2006;Valdemarsen et al., 2010;Rampazzo et al., 2013;Bannister et al., 2014).
In this research, we model the potential effects of increased carbon deposition related to filtration by fouling fauna in OWFs on benthic mineralization processes. To do so, we used output generated by Ivanov et al. (2021) who implemented a model describing filtration by blue mussels on turbine substructures in a hydrodynamic model of the Southern North Sea, including the full BPNS and parts of the bordering French and Dutch coastal zones (Figure 1; Ivanov et al., 2020). In biomass, M. edulis is by far the dominant species present on turbine foundations (Joschko et al., 2008;Krone et al., 2013). Whereas alterations to sediment and OM transport to sediments from the water column under different scenarios of turbine placement in a new concession area are evaluated in an accompanying paper (Ivanov et al., 2021), we used the OM deposition fluxes as input for a model of sediment biogeochemistry (OMEXDIA, Soetaert et al., 1996a), to investigate potential changes to mineralization pathways in the sediment. This was done for three cases: the first containing the current OWFs present in the model domain, the second a likely future scenario starting form 2026 when additional OWFs have been constructed and fully colonized in planned concession zones in the BPNS, and the Dunkirk OWF in the French coastal zone (Figure 1), and a third case containing no OWFs as a point of reference.
Seeing the analogy with aquaculture developments, we expected to find increased carbon concentrations in the OWF sediments through increased OM input, and increased anoxic mineralization, and nitrogen removal from the sediment through denitrification.

MATERIALS AND METHODS
The response of the sediment metabolism to alterations in carbon deposition induced by OWFs was represented by a dynamic model of early diagenesis (OMEXDIA, Figure 2; Soetaert et al., 1996a,b). The diagenetic model was forced by the carbon deposition estimated by Ivanov et al. (2021). This carbon deposition was estimated under different scenarios of OWF development (and without, see section "Study Area and Scenarios"), by implementing a representation of the filtration and feces production processes associated to monopile fouling FIGURE 1 | Median grain size (µm) of sediments in the modeled domain, with the Belgian Part of the North Sea (BPNS) delineated in black. Current OWF developments (2021) outlined in purple, and future developments (operational ± 2026) in orange, with letters corresponding to concession zone characteristics in Table 3. Black dots correspond to stations sampled for sediment biogeochemistry by Toussaint et al. (2021), used as reference material. mussels in a three-dimensional hydrodynamic-wave-sediment transport model (COAWST; Warner et al., 2010). In essence, COAWST models the distribution of OM and sediment in the water column, which are subject to hydrodynamics, depositionerosion processes, and simple degradation processes. A part of the OM is deposited and incorporated on the sediment. This part is what is used as input for the OMEXDIA model, which calculates how OM in the sediment is remineralized to free nutrients. The coupling between the 3D hydrodynamic model and the 1D model of sediment biogeochemistry was performed only in one direction (i.e., the water column provides the carbon flux to the sediment but no feedback from the sediment to the water column is considered). The spatial grid consisted of 91 by 101 cells (see Figure 1) covering the BPNS and surrounding waters, resulting in 9,191 dynamic implementations of OMEXDIA, each involving 50 vertical layers. Lateral exchanges between the diagenetic models are considered to be irrelevant.

Description of Pelagic Organic Matter
The annual cycle of OM deposition was modeled using the COAWST model framework (Warner et al., 2010), consisting of: a three-dimensional hydrodynamic model (ROMS; Regional Ocean Modeling System), a module describing OM and mineral particle dynamics through the water column and sediment bed (CCSTM; Community Coastal Sediment Transport Module), and a module for wave activity (SWAN; Simulating Waves Nearshore) (Supplementary Figure 2). The COAWST model was solved with a horizontal resolution of 5 km, downscaled to 1 km over the BPNS using a two-way nesting procedure (Ivanov et al., 2020). The ROMS module was configured and validated for the modeling domain by Ivanov et al. (2020), whereas implementation and validation and the CCSTM and SWAN modules were performed by Ivanov et al. (2021). For clarity, fundamentals of the CCSTM and SWAN modules are noted below, while for details we refer to Ivanov et al., in press.

Inclusion of OM Transport
The sediment module (CCSTM) was set up to describe the dynamics of three size-classes of mineral particles (i.e., silt, medium sand, coarse sand, with median grain sizes of 4, 22, and 750 µm respectively), and of two classes of OC of differing degradability (Fast and Slow degrading OC, resp. FDOC and SDOC, the inputs to the sediment diagenesis model). FDOC and SDOC originate from primary production at the water surface, which is forced by the daily vertically-integrated primary production product delivered by the Copernicus Marine Environment and Monitoring Service (Copernicus Marine Service Information, E. U, 2020).

Implementation of Filtration by Mytilus edulis
Filtration of OM by the blue mussel M. edulis on turbine foundations was implemented in grid cells containing OWFs as a conversion of silt, and the two classes of OC (FDOC, SDOC), to fecal silt, fecal fast, and fecal slow degrading matter (Ivanov et al., 2021). These fecal substances have exactly the same characteristics as the "non-fecal" forms, except for the sinking speed, which is 1,000 times higher (18 mm s −1 ) to account for the increased size and mass of the fecal pellets (Callier et al., 2006).
The filtration rate was determined by the M. edulis biomass in each grid cell, directly dependent on the available substructure surface in the cell. For this, all turbine foundations were assumed to be monopiles (including both C-Power concession zones; Table 3). All mussels were parametrized to be individuals of an average size [300 mg dry mass, Bayne et al. (1993)], occurring at a density of 6,468 ind. m −2 of substructure surface over a depth interval of 6 m starting from the water surface (Kerkhof F., personal communication). The M. edulis life cycle and seasonal biomass dynamics were not taken into account. The daily OC flux to the sediment from fecal pellets was added to the daily carbon flux stemming from sinking of primary production, to achieve the net carbon flux to the sediment for each grid cell, for each day of the year (Figure 2).

Sediment Early Diagenesis Model Description
Mineralization processes in the sediment layer were modeled using a dynamic implementation of the early diagenesis model OMEXDIA (Soetaert et al., 1996a,b). In this model, the concentrations of two classes of reactive OC present in OM (fast decaying detritus FDET, and slow decaying detritus SDET), oxygen (O 2 ), nitrate (NO 3 − ), ammonium (NH 4 + ), dissolved inorganic carbon (DIC), and oxygen demand units (ODUs, reduced substances generated in anoxic mineralization) are described on a 1D grid. This vertical grid has 50 layers of increasing thickness, starting from 0.01 cm at the sediment-water interface (SWI), extending down to a sediment depth of 200 cm.
In each sediment layer, the fractions of organic detritus (carbon) are mineralized to DIC and other reaction products through oxic or anoxic mineralization, or in denitrification (resp. reaction 1-3 in Table 1), with these processes mediated by concentrations of O 2 , and NO 3 − in each layer through first order or Monod reaction kinetics (Soetaert et al., 1996a). Oxygen is additionally consumed in nitrification (the oxidation of ammonium), and the reoxidation of reduced substances produced in anoxic mineralization (reaction 4-5 in Table 1). To investigate changes in OM reactivity, the average reactivity was calculated as follows (Eq. 1): where rFast and rSlow are the degradation rates (d −1 ) of fast and slow decaying detritus. Transport of modeled substances between layers is caused by advection (sediment accretion), molecular diffusion (for solutes), and bioturbation (for solids). Bioturbation was implemented as a depth-dependent biodiffusion coefficient (Db, cm 2 d −1 , Eq. 2), starting from the surface biodiffusion coefficient (Db 0 ), remaining constant for a certain mixed depth L mix , and decreasing below that layer to 0 according to a constant attenuation coefficient (Db coeff ).
The magnitude of the surface biodiffusion coefficient Db 0 , was modified in each grid cell i as dependent on three temporally and/or spatially variable factors: bottom water temperature, sediment median grain size, and the incoming carbon flux (Eq. 3).
The temperature effect assumed a doubling of organismal metabolism for every increase in temperature of 10 degrees (Q 10 of 2, Eq. 3, Davis and McIntire, 1983;Wrede et al., 2018). Median grain size (MGS, retrieved from the COAWST model, Ivanov et al., 2021) was used to scale the surface biodiffusion coefficient (Db 0 ), consistent with its relation to the calculated bioturbation potential of the species community [BP c , Solan et al. (2004)]. The BP c -MGS relation was based on measurements of Braeckman et al. (2014) and Toussaint et al. (2021). BP c values were highest for median grain size values of 100-175 µm, which represent fine sandy sediments in the BPNS, and were low where sediments consisted mostly of silt (MGS < 63 µm) and were too anoxic to allow for deep living species (Pearson and Rosenberg, 1978), and in coarse sands starting from MGS > 400 µm (see Supplementary Figure 1). This relation was parametrized as a skewed normal distribution (dsn, see Supplementary Information). Benthic activity was also assumed to be stimulated by high incoming OM flux (Cflux). This was represented by scaling Db 0,i according to the ratio between the actual deposition flux and the maximal deposition flux at the same location (maxFlux i , Eq. 3, Dauwe et al., 1998;Zhang and Wirtz, 2017). The mixed layer depth (in cm) was estimated following Zhang and Wirtz (2017), using the ratio of fast decaying OC relative to total OC deposited (pFast) to represent the OC quality or palatability (Eq. 4).

Process Reaction
Oxic mineralization x denotes the molar C:P ratio, y the molar N:P ratio in organic matter per mole of phosphorus (for Redfield Stoichiometry, x = 106, y = 16).

Parametrization and Implementation
Porosity and mean grain size were extracted from the COAWST model for each grid cell. Daily fluxes to the sediment of fast degrading OC, and slow degrading OC ). An advection rate w of 0.1 cm y −1 , an average value for a dynamic coastal sea (Boudreau, 1997;Mouret et al., 2009) was used for the entire domain.
The biogeochemical model 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). Molecular diffusion coefficients were calculated using R-package marelac (Soetaert and Petzoldt, 2018). The initial conditions for a dynamic run were obtained by the steady-state solution, using the R-package rootSolve (Soetaert, 2009).
Several runs of COAWST and OMEXDIA were done to finetune model parameters. During these calibration experiments, the SCOC and DIC fluxes produced by OMEXDIA forced by the COAWST carbon deposition were compared with field observations collected at six reference stations (Provoost et al., 2013;Toussaint et al., 2021). Parameters that were calibrated were degradation rates for OC in the water column and in the sediment, the sinking rates of OC and fecal pellets, and the base bioturbation value Db 0 . Details of the validation are provided in the Supplementary Table 1, and in Ivanov et al. (2021).

Study Area and Scenarios
Three simulations were performed, corresponding to different states of OWF development in and around the BPNS (Figure 1). Simulation "baseline" represents the situation before 2008, when no OWFs were present. Simulation "current" is the current situation, in which all turbines in the eastern concession zone in the BPNS are constructed (399 turbines, density of 1-5 turbines km −2 ), as well as the Borsele OWFs in the Dutch coastal zone (173 turbines, density of 0.75-1 km −2 ) bordering the BPNS ( Table 3). The number of turbines in the concession zones was used to estimate their density in the model (Table 3). Simulation "future" is a likely scenario after 2026 (Figure 1), in which the western concession zone in the BPNS (210 turbines) is developed with a turbine density of 2 km −2 , with as little as possible intrusion in the Natura 2000 area "Vlaamse Banken" (Figure 1). Intrusion into Natura 2000 areas may not be allowed, and requires compensation measures (EC, 2001). In this "future" scenario, the planned OWF in the French coastal zone (46 turbines, two turbines km −2 ) bordering the BPNS is also included. Planned capacity  Turbine foundation types correspond to monopile (M), gravity based foundation (G), jacket foundation (J).
and number of turbines for the "future" scenario were retrieved from 4C offshore (2020). Modeling results were aggregated over different scales: over the OWF area (OWF locations, grid cells containing OWFs-an area totaling 442.5 and 636.6 km 2 in the "current" and "future" scenario, respectively), over the full modeling domain (9,140.4 km 2 ), and for the BPNS (3,366 km 2 ). Effects on sediments outside of the OWFs were also quantified, by considering only grid cells without OWFs in either the BPNS, or the full modeling domain (Outside OWFs). Changes to biogeochemical processes in these different spatial scales were quantified relative to corresponding cells in the baseline scenario. The duration of the dynamic biogeochemical model simulations was set to 20 years to evaluate the effects of OWF placement near the end of their expected lifetime of 20-25 years. Output of the last year of the model simulations is presented in the "Results and Discussion" sections, averaged throughout the year and represented as a daily rate in line with the constant mussel biomass on the turbine substructures, and to focus on spatial differences across the modeling domain.

The Baseline Scenario
The baseline simulation (no OWFs) showed a gradient of decreasing OM processing away from the shore, with averaged total mineralization rates in excess of 40 mmol C m −2 d −1 nearest to shore, and lowest rates (5-10 mmol C m −2 d −1 ) in the northwestern part of the model domain ( Figure 3A). This pattern largely followed the grain size distribution (Figure 1). The bulk of the OC mineralization occurred through anoxic mineralization, which followed the same patterns as total mineralization ( Figure 3B). Rates of oxic mineralization related inversely to the grain size (compare Figure 3C with Figure 1), with higher total oxic mineralization rates in coarser sediments on the sandbanks in the center of the BPNS and in the northwest of the model domain, and lowest rates in finest sediments nearshore, and in the northeast of the modeling domain ( Figure 3C). Denitrification rates were highest near the outflow of the Scheldt estuary (0.4-1.1 mmol C m −2 d −1 ), decreasing to below 0.4 mmol m −2 d −1 20-30 km away from this source ( Figure 3D).

Within OWF
Rates of all mineralization processes increased within the OWF areas both in the "current" and "future" scenarios, resulting in a higher total mineralization of OC within these areas. On average, total mineralization rates within OWFs were 26.8 and 29.5% (values reported as "current" and "future" scenario) higher compared to baseline conditions, with maximal increases in individual cells of up to 69% (Figures 4A,B). This was in line with the similar increase of the OC flux (Table 4).
Increases in total mineralization rates were predominantly supported by increased anoxic mineralization (Figures 4C,D), followed by oxic mineralization (Figures 4E,F). Denitrification increased by 2.3 and 2.5% on average (Figures 4G,H and Table 5), resulting in additional removal of 7 and 10.7 kmol y −1 N respectively. These higher denitrification rates within OWFs originated from the increased OM influx. When considering the relative contribution of denitrification to OM respiration, instead the proportion of denitrified OM decreased by 17.8 and 19.6% ( Table 5).
Less of the incoming OC flux was returned to the water column as DIC, with this conversion decreasing by −3.6 and −4.5% in the "current" and "future" scenario, respectively ( Table 4). Consequently the amount of reactive carbon stored in sediments of OWFs increased. The amount of OC in the upper 10 cm increased by 9.7 and 11% on average (Table 4), and up to 37% in individual cells (Figure 5). The constant deposition of fresh OM additionally increased the reactivity of the OC within OWFs by similar values (Figure 5).

Outside OWF
Offshore windfarms also affected biogeochemistry in sediments beyond their area ("Outside OWF full domain, " Tables 4, 5), with effects mostly in the opposite direction as observed for the within-OWF changes (Figure 4). Outside of the OWF area, the incoming OC flux decreased for both scenarios (0.2 and 0.3% vs. baseline). Likewise, the total mineralization rates were lower (0.1 and 0.1% vs. baseline), with associated decreases in individual mineralization processes (Figures 4A,E,G,D). In the "future" scenario, oxic mineralization increased also across most of the modeling domain outside of the OWFs, while denitrification increased only for sediments closest to shore (Figures 4F,H). Overall, this decreased the contribution of denitrification to total mineralization (part denitrified) in the future "scenario" by 0.2%, with proportionally less N removed as a result, whereas in the "current" scenario this proportion increased by 0.1%. The average storage of reactive carbon outside of the OWFs decreased by 0.4% in the "current" scenario (Table 4), and decreased more northeast of the OWFs (Figures 4, 5A,C). For the "future" scenario (Figures 4, 5B,D), carbon storage increased by 0.8% (Table 4).
Considering only the BPNS, patterns are reversed ("Outside OWF BPNS, " Tables 4, 5). OC deposition fluxes increased slightly outside of the OWF area in both scenarios, by 0.3 and 0.4%, resulting in increased mineralization rates of 0.3 and 0.6% ("current" and "future"). The conversion efficiency of incoming OC to DIC on the other hand, decreased by 0.2 and 0.6%, relative to the baseline scenario (Table 4). This was associated to an increase in the amount of carbon stored in the sediments of the BPNS outside of OWF grid cells of 0.3 and 2%.

Full Domain Effects
The effect of the presence of OWFs in the entire modeling domain was an average increase in the total benthic mineralization of 1.4% in the "current" scenario, increasing to 1.9% in the "future" scenario, and an increased reactive carbon content of the sediment by 0.2, and 0.5% (resp. "current" and "future, " Table 4). Total denitrification also increased over the full domain by 0.03 and 0.2%, whereas the proportion of denitrified organic material decreased by 0.6 and 1.4% ( Table 5). As a result of this balance between enhanced benthic respiration and decreased contribution of denitrification, the total annual N removal increased by 0.1 and 0.2%. Within the BPNS only, the direction of changes to individual processes was similar to patterns seen in the full 4 | Averaged values and relative changes between tested scenarios and the baseline of carbon mineralization processes for the current (C) and future (F) scenarios, subdivided in values over locations with OWFs, the BPNS excluding OWFs (far-field BPNS), the domain excluding OWFs (Far-field domain), the BPNS, and the full domain.

Total mineralization
Org. C flux DIC flux Conversion Carbon storage Values for total mineralization, organic carbon flux, and DIC flux in mmol m −2 d −1 , the conversion represents the proportion of the incoming organic carbon flux released as DIC, and carbon storage is the total reactive carbon in the top 10 cm of the sediment (mol m −2 ).

Scope of the Results
The work shown here represents a first investigation of the large-scale effects of blue mussels on offshore wind substructures on OM mineralization in the sediment. With a mechanistic implementation of the filtration activities of M. edulis in a coupled hydrodynamic-biogeochemical sediment model, we estimated potential changes to mineralization intensities and pathways for two scenarios of OWF development in a coastal sea. That said, the limited availability of validation material for the baseline state, and some particular assumptions made in the model represent difficulties in assessing the robustness of the results. Top down, the mineralization of OM in sediments depends on the amount of OM that is deposited from the water column. In our approach, this amount is regulated by primary production in near the water surface, redistribution of OM by blue mussels when there are OWFs present, and the hydrodynamic regime of the area. In COAWST, the CMEMS primary production product was used to force primary productivity on the water column, whose deposition on the sediment surface is only as realistic as the modeled hydrodynamics forces that affect it in the water column. Overall, the COAWST model reproduced observed SPM distributions quite well (Ivanov et al., 2021), with exception of some nearshore, muddy areas where there appears to be a mismatch between modeled and actual bottom shear stress (Ivanov et al., 2021). As a result, modeled DIC and O 2 fluxes compared quite well with reference material (Supplementary Table 1), except in the previously mentioned muddy area, and in a location where the seabed is subject to periodic disturbances in the form of dumping of dredging material (station 130, Supplementary Table 1; van de Velde et al., 2018). However, these types of sediments do not occur in the domain where OWFs are placed, increasing our confidence in the model performance in those areas.
The lack of extensive validation data makes it impossible to numerically estimate the model error. Traditional methods for the validation of spatial models rely either on satellite products, or datasets with temporally coherent measurements (e.g., buoys) of variables with which the model output may be compared, and some type of error can be calculated. For sediment biogeochemistry, these types of measurements do not yet exist in the modeling domain. For now, available information derived from point measurements such as sediment grain size and porosity, OM sinking and degradation rates, bioturbation 5 | Averaged values and relative changes between tested scenarios and the baseline of carbon mineralization processes for the current (C) and future (F) scenarios, subdivided in values over locations with OWFs, the BPNS excluding OWFs (far-field BPNS), the domain excluding OWFs (Far-field domain), the BPNS, and the full domain.

Total denitrification
Part N removed Part denitrified N removed Values for total denitrification are in mmol m −2 d −1 , part N removed is the proportion of incoming N removed from the sediment, and part denitrified is the proportion of total mineralization occurring through denitrification.
activity (related to benthic biomass), and nutrient exchange rates were used as input values or as reference material (e.g., Verfaillie et al., 2006;Braeckman et al., 2014;Provoost et al., 2013;Zhang and Wirtz, 2017;Toussaint et al., 2021). The challenge of validating benthic-pelagic coupling models is well known (e.g., Luff and Moll, 2004;Griffiths et al., 2017). Promising technical developments aimed at facilitating the comparison between insitu observations and model outputs include eddy covariance derivation of benthic nutrient and oxygen fluxes, which allows for an increase of the spatial footprint of measurements (Berg et al., 2003), and regular benthic sampling through automated benthic stations, which can better constrain the temporal variability of in-situ flux acquisition (e.g., Toussaint et al., 2014;Moriarty et al., 2017). The implementation of the filtration effect of M. edulis and its effect on modeled benthic mineralization processes was similarly based on available reference material. The mussel biomass in the model was constant throughout the year, and mussels were assumed to be adults with an average dry biomass of 300 mg DW ind. −1 . Reported biomass of M. edulis on offshore constructions varies strongly in time and space, but is mainly constrained between 110 and 768 mg DW ind. −1 (e.g., Bouma and Lengkeek, 2012;Krone et al., 2013, converted values). The chosen biomass of 300 mg DW ind. −1 thus represents slightly smaller sized individuals compared to the known average, and corresponds to the average size of individuals used in the experiments of Bayne et al. (1993) in which M. edulis filtration behavior was characterized. Since predictions of filtration effects, including those presented here, hinge on this important value, it would serve future impact assessments well to collect more data about the community composition and seasonal dynamics of M. edulis and other organisms on offshore structures. For example, the potentially large contribution of the amphipod J. herdmanii to OM assimilation of turbine substructures has recently come to light (Mavraki et al., 2020). Similar to Slavik et al. (2019), seasonal variation in mussel biomass was not implemented despite a tendency of biomass to increase toward summer (Krone et al., 2013). This is, however, not a consistent feature [see Bouma and Lengkeek (2012)], and would have introduced the need to dynamically model the biomass of mussels on the turbine substructures. A dynamic description of the mussel biomass would only have been possible with a feedback loop between the sediment and the water-column. Without such a feedback loop, the one-way mass transfer to the sediment in the shape of dead M. edulis biomass would introduce an erroneous distribution of carbon in the model. This does not mean that there was no seasonality in the production of fecal pellets; phytoplankton concentrations in the water column varied according to the CMEMS PP product, and the filtration rate is metabolically scaled according to the ambient temperature (see Ivanov et al., in press). Both vary throughout the year.
Specifically, feedback from OMEXDIA to COAWST was omitted to manage the computational load, while still representing biogeochemical processes in the sediment in high resolution. Thus, we implicitly assumed that primary production in the water column would not be affected by potential changes in ambient nutrient concentrations resulting from the sediment, whereas in reality this will definitely be the case (Floeter et al., 2017). Specifically, a dynamic coupling with a water column biogeochemistry model (Luff and Moll, 2004;Slavik et al., 2019), would make it possible to investigate whether water column nutrient concentrations in OWFs can be expected to decrease, or increase as a result of the transfer or OM to the sediment, and whether this in turn affects phytoplankton concentrations. Especially in shallow shelf seas, nutrients regenerated in the sediment co-regulate primary production events in the water column (Ruardij and Van Raaphorst, 1995;Soetaert and Middelburg, 2009). A negative feedback on phytoplankton blooms for example, would decrease fecal pellet production and subsequent storage of OC in the sediment, thus introducing self-regulation in the system. Which way this feedback may go is a pathway of future investigation, as this chain of events will depend on the complex interplay between water depth, hydrodynamics, nutrient regeneration, and filtration by fouling fauna. So in summary, the presented model results are based on available literature data, and a pragmatic methodology allowing for a first investigation of these complex processes, but absolute values should be interpreted with the previously mentioned remarks in mind.

Strong Local Effects and Consequences for Surrounding Sediments
Modeled impacts of OWFs extended over the full modeling domain, with strongest effects on sediment biogeochemistry within the OWF itself, and weaker but clear effects on surrounding sediments. Deposition of highly reactive OC increased as a result of fecal pellet production by M. edulis. Total benthic mineralization within the OWFs thus increased by 26.8-29.5% in two scenarios that either considered the "current" situation, where OWFs comprised 443 km 2 , or a "future" situation where the domain of OWFs was increased to 636 km 2 . At the same time, OM deposition and mineralization decreased in regions surrounding the OWFs. Rates of oxic mineralization, denitrification, and anoxic mineralization as well as the carbon storage all decreased in sediments surrounding the OWFs (Figures 4, 5), with effects extending to the boundaries of the modeling domain, suggesting impacts beyond the model boundaries. In all likelihood, impacts extend beyond the artificial modeling domain, which was chosen to represent the BPNS, until net changes are evened out. Indeed, in this model run, the OWFs merely redistribute the same amount of produced OM, as there was no feedback provided to the water column and subsequently the rates of primary productivity.
For both scenarios of OWF development, extended halo's of altered benthic processing were mostly situated to the northeast of the OWFs (e.g., see Figure 4A), in line with the main direction of the residual current present in the region (Supplementary Figure 3: current rosette; Ivanov et al., 2020). Since OWFs act as filters of OM that locally shunt OM from the water column to the sediment, water transported further NE became relatively depleted of OM compared to the baseline (Ivanov et al., 2021). As a result, OM deposition and subsequent mineralization processes decreased in the wake of the OWFs.
The "future" placement of additional OWFs in the west of the BPNS and in the French coastal zone toward 2030 changed mineralization patterns considerably compared to the "current" scenario. With increasing OWF area, the effect halo's of different OWFs begin to overlap, with more pronounced effects throughout the domain (Figure 4). Most striking was the increase of oxic mineralization in the majority of sediments in the modeling domain ( Figure 4F). Within the OWFs, OM is concentrated, resulting in higher deposition, which is partly mineralized by oxic mineralization, but mostly by anoxic mineralization. This combined effect of OM trapping in multiple OWFs decreased the OC load on the sediments outside of the OWF area to a point where most of it is mineralized in surficial sediment layers, by consuming oxygen.
Denitrification was also affected by the expansion of the OWF area in the model domain, showing on one hand an increase of denitrification in nearshore sediments, and on the other hand a decreasing denitrification further away from the coastline (Figure 4H). This interesting pattern can be explained by a shift from high to more intermediate OM loadings discussed previously. Denitrification is known to peak at intermediate carbon loadings (Soetaert et al., 2000), and at high nitrate concentrations, conditions that are met nearshore due to high NO 3 − loadings of the Scheldt outflow. While the zones of denitrification shifted spatially, the overall denitrification throughout the entire modeling domain barely increased (0.03-0.17% for "current" and "future" scenario, Table 5). From this, we derive that the amount of N removed from the aquatic environment is only marginally increased by this scenario, with total increases of 3.2 and 13.1 kmol N on an annual basis, which represents about 0.1-0.4% of the annual N input of the Scheldt estuary (45 kT N y-1; Brion et al., 2006).
Whereas the physical effects of wind turbine foundations on the water column have so far only been noticed within 1-5 km away from OWF areas (Rivier et al., 2016;Floeter et al., 2017), the scale over which e.g., water filtration, feces production, deposition and resuspension interact to alter OM dynamics is still relatively unknown. Our results point to a region of influence that is of similar magnitude as the effect on water column primary productivity (NPP). Slavik et al. (2019) described a decrease of NPP of up to 8% on an annual basis within the OWF and further reductions up to 20 km away. Conversely, NPP increased up to 50 km away. It is thus increasingly likely that the placement of OWFs has far-reaching effects on benthic-pelagic ecosystem functioning.
These spatial effects raise the issue of transboundary impacts. Biogeochemical alterations induced by OWFs installed in one country can have an impact on the functioning of the benthic ecosystem in another country. The region on the western boundary of the BPNS, including part of the newly foreseen OWF expansion area, is a Natura 2000 area (Figure 1), an area containing threatened habitats and species protected under European law. Among others, this area contains rare gravel bed habitats, which are highly sensitive to anthropogenic disturbances, but are of great importance for species that require hard substrates for reproduction such as the common whelk, the Atlantic bobtail, and the spotted catshark [respectively Buccinum undatum, Sepiola atlantica, Scyliorhinus canicula; Degraer et al. (2009)]. For example, increased particle concentrations originating from OWF can clog the filtering apparatus of filter feeders, or limit visibility for visual predators in this valuable Natura 2000 area (Essink, 1999). This region, already directly impacted by the placement of OWFs in the North, and by domain-wide effects caused by the compounding effects of multiple OWFs (Figures 4, 5), additionally may be impacted by the placement of the Dunkirk OWF, whose increased particle production is carried into the BPNS by the residual current.
Whereas the EU is at the forefront of transboundary marine spatial planning (Li and Jay, 2020), the transboundary biogeochemical effects of OWFs are not yet being considered in projects tasked with investigating this complex problem (European MSP platform, 2020). By indicating the potential transboundary effects that the Dunkirk OWF may have on the Belgian Natura 2000 area, our results indicate that this will be necessary with future expansions of OWFs.

Increased Carbon Storage
In the model, carbon cycling in the sediment was altered by the presence of OWFs. Twenty years after the installation of OWFs, the deposited OC was returned less efficiently to the water column within the OWF areas, resulting in an increased buildup of OC in the sediment ( Table 4). The additional deposition primarily translated to an increase in anoxic mineralization rates (by max. 75% locally), with oxic mineralization increasing less (30-40%, Figure 4). With higher deposition, concentrations of OC in the upper 10 cm of sediment increased by 10-11% within OWFs, while this is 0.2-0.5% when extrapolated over the entire area covered by our model. As such, OWFs installed on the permeable sediments of the Southern North Sea became local sinks for carbon relative to the baseline scenario. As such, the current and future windfarm concession zones in the BPNS increased the total amount of reactive carbon trapped in the upper 10 cm of the sediment by 28,715-48,406 tons of carbon respectively, coinciding with 0.014-0.025% of Belgium's annual greenhouse gas emissions (118.5 million tons in CO 2 equivalents in 2018, VMM et al., 2020). In that sense, the carbon potentially stored in OWFs represents a small but significant carbon offset in carbon accounting, which can be added to the intrinsic carbon economy of OWFs. In Belgium, an energy production of 3391 GWh was attributed to OWFs in 2018 (at 0 tCO 2 GWh −1 ; operational cost only; Belgian Debt Agency, 2018). Compared to a gas turbine for example (380 tCO 2 GWh −1 ), this corresponds to 1288580 tCO 2 not emitted by using wind-generated power. When including the life cycle (material sourcing, construction, . . .), this factor ranges between 14 and 111 (depending on the reference chosen; Dolan and Heath, 2012). This would reduce CO 2 emissions between 1041037 and 2858613 tCO 2 . To this, our estimated quantities in the sediment contribute an additional 1, to 4.6 %.
Naturally, this storage is of limited duration. Firstly, sediment disturbing activities such as bottom trawling may or may not be prohibited within OWFs (e.g., the United Kingdom and France allow trawling outside a 50 m radius around individual turbines, whereas in Belgium, Netherlands, and Germany it is prohibited completely). Bottom trawling will resuspend sediments, likely lowering the build-up of OM as presented here . Secondly, the expected lifespan of a wind turbine is 20-25 years (Nghiem and Pineda, 2017) after which, in theory, the concession zone needs to be restored to its original condition by the concession holder (Kruse et al., 2019). In case of such full decommissioning, the increased OM input will cease while the decommissioning activities themselves can trigger the release of accumulated carbon to the water column as they will most likely include sediment disturbing activities. However, when alternative partial decommissioning scenarios would be considered whereby part of the subtidal structure would either remain in place, be repurposed, or translocated (Fowler et al., 2020); OM filtration by fouling fauna and subsequent local carbon storage in sediments should be considered an important decision criterion.

Growth Urges Research
The currently installed wind energy capacity in the North Sea is 18.5 GW, and is expected to rise to 70 GW by 2030 in order to achieve 30% of renewable energy targeted by the EU (Nghiem and Pineda, 2017;Selot et al., 2019). Thus under the impulse of the ongoing "blue acceleration, " the rapid development of the ocean economy (Jouffray et al., 2020), new offshore energy provisions are expected to contribute heavily to the expansion of man-made structures at sea collected under the term "ocean sprawl" (Duarte et al., 2013;Bishop et al., 2017;Jouffray et al., 2020). Concurrently, decommissioning of OWF developments, as well as of oil and gas installations reaching their end-of-life is expected to increase strongly in the coming decades (Kruse et al., 2019;Fowler et al., 2020). This means that despite recently intensified research efforts, a full cycle of OWF development to decommissioning will take place before the ecological effects are fully known. Recently, several knowledge gaps on the ecological effects of decommissioning of offshore structures have been identified (Fowler et al., 2020), but possible effects on benthic ecology extending beyond the OWF areas are still overlooked. Even then, additional factors such as the ongoing expansion of OWF's to deeper waters (Ramirez et al., 2020), where biotic and abiotic factors are different still, will again lead to a novel set of feedbacks to the ecosystem.
Despite these current shortcomings, modeling studies such as ours generate new hypotheses to validate or falsify in future research on the effects of OWFs on biogeochemistry in the sediment (e.g., do changes in the benthic community around offshore structures affect denitrification rates differently than those modeled as a result of only increased OM deposition and solute availability), and indicate where research should be performed to improve on models (e.g., OM deposition rates and quality). Furthermore, this approach is not limited to the chosen model domain. Both COAWST and OMEXDIA are flexible models that have been applied to a variety of marine settings [OMEXDIA: (e.g., Soetaert et al., 1996a;Laurent et al., 2016;Khalil et al., 2018;De Borger et al., 2021), COAWST: Warner et al., 2010;Fennel et al., 2013;Bolaños et al., 2014;Zang et al., 2020], including other regions where OWF development may take place. So for now, we have provided a simplification of a complex ecological chain of events, that may be used as a basis in future model developments.

CONCLUSION
Water-column filtration induced by the presence of large densities of filter feeding fouling fauna on offshore wind turbines caused the displacement and aggregation of fresh OC from the water column to sediments. We show that sediments in OWFs become sites of intense OC mineralization, with a halo of reduced mineralization outside of the OWF areas. With this, the body of evidence is growing that OWFs have both local and larger scale ecosystem effects. Field measurements of sediment functioning near wind turbines are needed to give weight to our results and guide future model developments. Given the current proliferation of OWFs in North Sea waters and other shelf seas around the globe, more studies toward the positive and negative ecosystem effects of these structures are urgently needed.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
ED implemented and performed biogeochemical model simulations, analyzed results, and prepared the manuscript. EI implemented and performed COAWST model simulations, and contributed to the manuscript. AC, UB, JV, MG, and KS conceptualized the study and contributed to the manuscript. All authors contributed to the article and approved the submitted version.