A spatial model for nutrient mitigation potential of blue mussel farms in the western Baltic Sea

consistingofaspa-tial and


Introduction
Both European and international coastal and marine policies have the objectives to reach healthy ocean ecosystems, sustained marine ecosystem services, and integrated utilization of marine areas by minimizing conflicts in the marine environment (Ehler and Douvere, 2009;EU, 2014;EU, 2008;EU, 2000;IOC, 2014).In this respect, eutrophication is a global challenge, particularly impairing near-coastal marine water bodies (Smith, 2003).
Within the European Union, Environmental targets for reaching good ecological status in coastal water bodies, including target nutrient (nitrogen and phosphorus) and chlorophyll-a (ChlA) concentrations, were defined within the Water Framework Directive (WFD) (EU, 2000) and the Marine Strategy Framework Directive (MSFD) (EU, 2008).Marine aquaculture, on the other hand, is an essential element of implementing the Maritime Spatial Planning Directive (MSPD) (EU, 2014;Gimpel et al., 2018).Considerable reductions of nutrient contents in the western Baltic Sea have been achieved over the last decades, mainly attributed to land-based measures on reducing inputs from both point and diffuse sources (Carstensen et al., 2006;EEA, 2019;Riemann et al., 2016).However, large pools of internal nutrient loadings in the Baltic Sea delay recovery from eutrophication (Carstensen et al., 2006;Riemann et al., 2016), and most water bodies in the western Baltic Sea are still far from reaching the environmental targets (HELCOM, 2018).The Danish River Basin Management Plans for the years 2015-2021(EPA, 2016)), for example, include targets of nitrogen load reduction for coastal water bodies.These total up to 13.000 t N /year for all of Denmark, and locally, this target reduction can reach up to 33 t N / (km 2 *year) with respect to water body area (only water bodies larger than 5 km 2 considered).Consequently, in situ mitigation measures could effectively complement land-based actions to reduce nutrient contents in coastal water bodies (Petersen et al., 2014).This would contribute to minimizing threats of eutrophication to valuable coastal ecosystems, such as the loss of benthic macrophytes (Timmermann et al., 2019), or frequent occurrences of harmful algal blooms (Delegrange et al., 2015).However, potential conflicts about marine space with other economic sectors (e.g.fishing and marine transportation) may arise when implementing marine mitigation measures (Maar et al., 2020b).
Mussel mitigation farming, a form of marine aquaculture optimized for nutrient removal, has been proposed as a tool to improve coastal water quality (Petersen et al., 2014).It is designed for acting as effective in situ mitigation measure in eutrophic coastal water bodies or as compensation measure for marine point sources (i.e.finfish culture) (Maar et al., 2020a;von Thenen et al., 2020).In contrast to farming mussels for human consumption, where the mussels should be large and of similar size, the target of mitigation farms is to fix the highest amounts of nutrients possible in a short time.Such floating farm systems provide a large suspended substrate area for mussel larvae to settle and grow.Once settled, the mussels can filter large volumes of water, immobilizing phytoplankton biomass and other types of suspended particles.In eutrophic systems, nutrients are transformed into phytoplankton biomass, which is then transformed into mussel biomass after filtration, and extracted from the marine environment by harvesting after certain growth periods (Petersen et al., 2019a).Depending on the productivity of longline mitigation farms, respective prices for nitrogen and phosphorus reduction have been estimated between 13-53 €/ kg N and 180-880 €/kg P for areas with high and low productivity, respectively (Maar et al., 2020b).Blue mussels, or more specifically mussels of the Mytilus edulis/trossulus complex, are widely present in the shallow areas of the western Baltic Sea (Kotta et al., 2020).There, primary spawning occurs from May until June (Riisgård et al., 2015), when planktonic larvae are released into the water and widely dispersed over a period of up to six weeks (Larsson et al., 2017).Provided that supply of wild mussel seeds is sufficient, mussel mitigation farms are currently considered a potentially effective part of measures to improve water quality in these coastal areas (Kotta et al., 2020).The western Baltic Sea is a highly heterogenic marine water body, isolated from the open Atlantic Ocean by the narrow entrance through Skagerrak between Denmark, Norway, and Sweden.It is characterized by relatively shallow depths of mostly b50 m, a pronounced salinity gradient from north-west to south-east, and a fine-structured coastline with numerous islands, fjords, and bays.
Estimates from several studies in the Limfjord system reveal that longline mitigation farms are able to extract 0.6-1.4t N /ha within one growth season from July till March (Petersen et al., 2014;Nielsen et al., 2016;Taylor et al., 2019) and new farm constructions, such as with nets, appear even more efficient for mitigation culture.Currently, there are only a few Geographic Information System (GIS) based studies and applications available with respect to site-selection of mussel mitigation farms (e.g.Bagdanavičiūtė et al., 2018;Gimpel et al., 2018;Kotta et al., 2020).These studies evaluated the spatial variability of mussel growth conditions, either based on fixed thresholds or seasonal average values of habitat factors, but do not consider their temporal dynamics.Mussel growth, however, is a temporal process and relies on the simultaneous interaction of growth-controlling parameters, such as temperature (Temp), salinity (Sal), and food availability.A derivation of specific growth estimates based on e.g.fixed thresholds or yearly averages is, therefore, highly critical in seasonally and spatially dynamic ecosystems, such as the western Baltic Sea.Sainz et al. (2019) have recently published a MSP study on mussel farms in southern California, where they point out the importance to consider temporal variability in environmental conditions for spatial planning of marine aquaculture.Furthermore, large numbers of blue mussels are capable of filtering large volumes of water of the majority of phytoplankton content (Cranford, 2019;Petersen et al., 2019b), which can lead to within-farm food limitation and reduced growth performance of mussels (Petersen et al., 2019b;Taylor et al., 2019).Food limitation is strongly related to the supply of suitable food into a farm.ChlA is a usual proxy for the available food for mussels (Petersen et al., 2019b), even though the species composition of phytoplankton can affect ingestion rates.Consequently, food limitation can be described as a function of hydrodynamics, ChlA concentration, and the actual ChlA ingestion capacity of the blue mussels present in a farm (Petersen et al., 2019b).The required hydrodynamics to avoid food limitation for specific farm setups do not find consideration in the recently available GIS-based site-selection approaches (e.g.Bagdanavičiūtė et al., 2018;Gimpel et al., 2018) and the spatial model for the Baltic Sea (Kotta et al., 2020).Although the existing assessments can predict the large-scale pattern of mitigation potential in the Baltic Sea region, these models do not take into account shortterm dynamics of growth conditions (e.g.seasonality in Temp or phytoplankton blooms) and mussel densities at farm-scale (e.g.substrate density in farms interacting with natural food availability).Therefore, the actual mitigation potential of mussel farms may largely deviate from their modelled full potential, and a suitable spatial model for site-specific implementation of coastal and marine policies in the western Baltic Sea is still lacking.Eco-physiological models, such as the Dynamic Energy Budget Model (DEB-model) (Maar et al., 2015;Sainz et al., 2019), or Bio-Energetic Growth models (Larsen and Riisgård, 2016) are capable to predict the growth potential of mussels in the field, given that in situ experiments for calibration/validation and input data of ambient water quality conditions are available in suitable temporal resolution.Nevertheless, such models are currently limited in their direct implementation into large-scale GIS applications due to complex setups and high demands of computing capacity.
In this context, both the BONUS OPTIMUS (bonus-optimus.eu)and the MuMiPro (mumipro.dk)projects aim at exploring and optimizing the potential of mussel mitigation farming in the western Baltic Sea from different perspectives, including site-selection, farming methods, and economic feasibility.
In this paper, we develop a modular modelling approach to upscale local point-estimates of growth and biomass of individual blue mussels to hypothetical mussel mitigation farms on a transnational scale involving Danish, German, and Swedish water bodies in the western Baltic Sea.Therein, we modularly linked growth experiments in the field, eco-physiological model simulations, large-scale monitoring data, and spatial models by means of empiric statistical model functions.The results, illustrate an evidence-based and comprehensive large-scale assessment on the spatial distribution of nutrient reduction potentials of blue mussel mitigation farming across the western Baltic Sea.

Overview of the modular modelling approach
The developed model contains four inter-linked modules, which use data from monitoring, field experiments, eco-physiological modelling, and bathymetry to estimate harvest and nutrient reduction potentials, as well as locally required hydrodynamic conditions for different setups of Mussel Mitigation Farms in the western Baltic Sea (Fig. 1).The required input and output parameters to run the final model are listed for each module in Table 1.All variables used throughout the model are summarized in the suppl.Information.
• Module 1 provides spatial (1 km 2 resolution) and temporal (monthly resolution) models of the long-term average pelagic habitat factors temperature (Temp), salinity (Sal), ChlA concentration along with estimates of their natural variability.This module is based on data from the operationally performed monitoring programs of Denmark (NOVANA), Germany (LUNG, LLUR), and Sweden (SharkWEB).All spatial model calculations operated on a model raster with 1000*1000 m 2 spatial resolution in the UTM32N coordinate system.The extent of the model raster was aligned with the spatial boundaries of monitoring data plus a buffer of at least 20 km (X: 400,000-964,000 m; Y: 5,948,000-6,570,000 m).For spatial evaluation, we used the European coastline shapefile (EEA, 2015) for defining landwater boundaries and the Exclusive Economic Zones from the Maritime Boundaries Geodatabase (Flanders Marine Institute, 2018).Both shapefiles were projected to the UTM32N coordinate system.

In situ monitoring data
We used data from the operationally performed monitoring programs of Denmark (DCE, 2018), Germany (LUNG, 2018;LLUR, 2018), and Sweden (SMHI, 2018) throughout the western Baltic Sea from 2007 to 2017 to cover ten mussel growth seasons.All position data were transformed to the UTM32N coordinate system.The spatial boundary of monitoring data was set to X-coordinates between 500,000-950,000 m and Y-coordinates between 5,950,000-6,550,000 m.Additionally, the Limfjord system was fully included ranging down to X-coordinates of around 450,000 m.In total, data was evaluated from 970 monitoring locations (Denmark: 338; Sweden: 541; Germany: 91) using the parameters Temp (°C), Sal (psu), and ChlA (μg/L).Monitoring data was filtered to the upper 10 m water column to represent the depth range relevant for suspended mussel mitigation farms.For spatial modelling purposes, the study area was classified into open and separated water bodies (Fig. 2), such as fjords and bays with only narrow connections to the open sea.For the whole study area, the number of stations monitored per month of the year ranged between 284 (ChlA in November) and 642 (Temp in August).In the open water bodies, the number of stations monitored per month with at least five years coverage ranged between 80 (Sal in December) and 229 (Temp in August).An overview of all monitoring data used in the model is presented spatially in Fig. 2 and quantitatively in Table_S 1.

Blue mussel growth experiments
From 2017 to 2019, in situ blue mussel growth experiments were performed on longline systems at 10 sites in the study area over one or two mussel growth seasons (July-May): Limfjord (Løgstør, Sallingsund, and Skive), Mariager Fjord (inner and mid), Vejle Fjord, As Vig, Flensburg Fjord, Roskilde Fjord, and Greifswald Bay (Fig. 2).All experiments relied on natural recruitment of mussel spat.Growth data were acquired at differing time intervals by sampling triplicates of 30 cm sections of the spat collector material.Biomass dry-weight (m bio [g DW ]) and shell length of the mussels (l mus [mm]) were measured as described in Taylor et al., 2019.Average values of m bio at the different time-steps were used to calibrate (2017-2018) and validate (2018-2019) both the eco-physiological DEB-model, as well as to validate the empirical mussel growth model (Module 2 in Fig. 1) used for spatial modelling.The longline systems were partly individual test-lines (Mariager Fjord, Vejle Fjord, Flensburg Fjord, Roskilde Fjord, Greifswald Bay), or test-lines integrated within commercially run mussel farms (Limfjord, As Vig).Spat collector material (5 cm polypropylene belts) was applied at 1 m depth and set up as 0.4 m wide loops with 0.3 m or 0.6 m spacing between and 2 m or 3 m depth-range.

DEB-model
The DEB-model is a mechanistic eco-physiological model describing the energy flow through the mussel in response to Temp, Sal, and ChlA (proxy for food) (van der Veer et al., 2006;Saraiva et al., 2011;Maar et al., 2015).ChlA (including planktonic and resuspended fractions) was used as a quantifier for food, because phytoplankton biomass is considered as the main food source for mussels in these micro-tidal, eutrophic areas with high ChlA concentrations and relatively low amount of resuspended seston concentrations (Filgueira et al., 2018;CEFAS, 2016).Due to the low expected concentrations, adverse effects of inorganic resuspended material on ChlA assimilation have not been considered in the model.Further, the food quality of phytoplankton expressed at the assimilation efficiency is higher compared to other food sources (Filgueira et al., 2018;Navarro et al., 1996).In the model, mussel filtration ceases at ChlA concentrations ≤0.5 μg/L due to valve closure (Dolmer, 2000).Filtration was kept unaffected by Sal assuming local adaptation (Riisgård et al., 2013a), although abrupt changes in Sal have been observed to induce short-term valve-closure (Almada-Villela, 1984).The ingested food is assimilated by a constant assimilation efficiency (KA) except at high ChlA concentrations N20 μg/L, where the assimilation efficiency, and hence growth, is reduced exponentially due to food over-saturation.The exponential function (see Eq. ( 5)) is based on net growth efficiency data from Riisgård et al. (2013b) with ChlA concentrations varying from 4.1 to 67.8 μg/L (n = 7, R 2 = 0.85, their Fig. 5, ignoring one value at lower salinity).The assimilated food goes to the reserve density (energy storage).A high energy density implies a well-fed mussel, whereas a low energy density is a sign of starvation.
Energy is allocated from the reserve density with a fixed fraction of 0.7 to somatic growth and somatic maintenance, whereas the other fraction is allocated to maturity or reproduction (gonads) and the related maintenance.The growth response to low Sal is described as an extra maintenance cost due to osmoregulation below a Sal threshold of 16.2 psu (Almada-Villela, 1984;Maar et al., 2015).Osmoregulation takes place over the surface of an organism and the maintenance cost is assumed proportional to the surface of the mussel, Temp and Sal (Kooijman, 2000;Maar et al., 2015).The extra maintenance costs were distributed to both somatic-and reproductive maintenance and the reserve utilization was increased to cover this expense.Spawning takes place above a certain threshold for Temp (9.6 °C) and gonad content.During severe starvation, the somatic and reproductive tissue can pay for somatic maintenance.The DEB-model results were calibrated with observations from four experimental test-line sites in Skive Fjord, Sallingsund, Løgstør, and Mariager Fjord (Fig. 2) for the growth season 2017-2018.Validation was performed against observations from experimental test-line sites in eight different areas (Fig. 2) for the growth season 2018-2019 (Taylor et al., 2020).Hydrodynamics were not directly part of the DEB-model calculations.However, as the mussels for calibration and validation have been growing on farm-like longline setups, the effect of mussel clusters, forming on the collector substrate, on competition for food were indirectly covered.
The validated DEB-model was used to predict the potential growth of mussels at 50 monitoring stations in Denmark (Fig. 2) using data from the above-mentioned Danish national monitoring program NOVANA.Each site had at least 12 samples per year from two successive years of Temp, Sal, and ChlA.Furthermore, the minimum depth for monitoring sites was set to 5 m.The modelled period covered growth seasons (July to May) for the years 2008-2017 while the number of seasons modelled varied from 1 to 9 depending on available data for each site.The DEB-model parameterization, as well as the predictions of biomass for 50 sites from 2008 to 2017, were used as input to build, respectively calibrate the empirical mussel growth model (Module 2) for the western Baltic Sea.

Farm setup -the 'standard mitigation farm'
A standard mitigation farm in our model is based on longline farms, already applied in the field (Petersen et al., 2014;Nielsen et al., 2016;Taylor et al., 2019).It consists of three sections with 30 longlines of 200 m length, each (Figure_S 1).Space between longlines is set to 8 m.Each section covers an area of A sec = 250 m * 250 m = 6.25 ha; the whole farm covers 18.8 ha.Along each longline, spat-collector material (5 cm polypropylene belts) is continuously deployed at a depth of 1 m below water surface with loops reaching down another 2-8 m (r dep ), dependent on a specified maximum value max(r dep ) and local bathymetry b.To ensure the collector loops do not touch the sediment, at least 1 m space is maintained from the bottom.The width of each loop is set to w loop = 0.5 m and its length is defined to be l loop = 2 * r dep + w loop long.The distance between two loops is kept flexible between 0.2 and 1.5 m and the loop interval on the longlines i loop , therefore, varies between 0.7 and 2.0 m.A standard farm in the final model is, consequently, setup flexible with respect to max(r dep ), b and i loop .A whole model farm contains a total collector length l col [m] of: where Consequently, at max(r dep ) of 8 m and i loop of 1 m, l col of the standard mitigation farm is in between 90,000-306,000 m for 4 m, respectively ≥10 m water depth.We did not limit potential farm areas by water depth, as there are examples, where longline systems are described feasible for large depths of up to 100 m (Goseberg et al., 2017;Langan and Horton, 2003;Mizuta and Wikfors, 2019).To our knowledge, however, there is no full-scale commercial or mitigation farming for blue mussels in offshore environments, yet.
In this paper, we present results for different setups of the model mitigation farm.In order to specify easily the respective farm layout we use a nomenclature combining the potential range of r dep from min to max ('2-max(r dep )') and i loop : e.g.'2-2_0.7'specifies a farm with r dep = 2 m and i loop = 0.7 m and '2-8_2.0'a farm with 2 m ≤ r dep ≤ 8 m and i loop = 2.0 m.

Bathymetry
We merged two bathymetry datasets: (1) A Danish bathymetry model (DMSA, 2012) and ( 2) the EMODnet Bathymetry DTM (EMODnet, 2018), which were both projected to the target UTM32N coordinate system in ArcGIS and had spatial resolutions of ~50*50 m 2 and ~100*100 m 2 , respectively.Both datasets were transferred to the model raster using the 'rasterize' function of the 'raster'-package in R. For each raster cell, we calculated the (1) mean bathymetry, (2) a selection layer for suitability for mussel farms, where depths ≥4 m covered at least the area needed for a standard mitigation farm (18.8 ha), and (3) the mean bathymetry of all depths ≥4 m to derive respective r dep for potential mussel farms.In this step, both raster layers were combined using the 'merge' function of the 'raster'-package in R, prioritizing the layer based on the higher resolved Danish bathymetry data.

Hydrodynamics
Hydrodynamics are important in this model, as food-flux into a mussel farm is assumed as a function of ChlA concentration and horizontal flow velocity (Petersen et al., 2019b).Here, monthly-mean absolute horizontal flow velocities serve as a reference for the required hydrodynamics to avoid food limitation (output of Module 4).Unfortunately, there is no such dataset covering the whole spatial extent of this study.We used monthly-mean flow velocities of the Baltic Sea Physics Analysis and Forecast (BSPAF) (EU Copernicus, 2019) for the upper 10 m water column over three available years from March 2016 to February 2019.The Limfjord system is not part of the BSPAF's model domain.There, we calculated monthly-mean absolute horizontal flow velocities from model runs with FlexSem (Larsen et al., 2020) for the years 2009-2012, and 2017 as reference (Larsen, 2020).

Module 1: Spatial and temporal habitat factor model
We evaluated the monitoring data of pelagic habitat factors with respect to spatial (1 km 2 resolution) and temporal (monthly) distribution on a long-term basis.Due to the specific value distribution of ChlA concentrations, we used ln(ChlA) for spatial data processing.First, the data was aggregated by 'Station', 'Year', and 'Month', and we calculated respective median values for all available parameters at each monitoring location.We chose the median as a robust measure towards individual outlier values.These values were further aggregated by 'Station' and 'Month' to calculate monthly long-term averages for each monitoring location over the whole time period from 2007 to 2017.Standard deviations and the number of years covered were calculated, respectively, for each month and monitoring location.Standard deviations were used for further assessment of effects of natural variability in pelagic habitat factors on expected mussel growth.
We assumed that the deviation of values from individual years from their respective long-term average is normally distributed.To justify this assumption, we ran Shapiro-Wilk test for normality (shapiro.testfunction, 'stats' R-package) for each parameter at each monitoring location.Therein, only values from months with at least five years coverage were considered.Rejection rates of normality (p-value b0.05) for Temp, Sal, and ln(ChlA) were 11%, 21%, and 32%, respectively, and sample sizes ranged in between 5 and 132.Due to the large majority of monitoring locations passing normality tests, we decided to accept this assumption for the whole model domain.

Open water bodies
For the open water bodies, continuous geostatistical modelling was applied to estimate long-term monthly average Temp, Sal, and ChlA conditions along with their standard deviations, across the whole open water body spatial selection.To account for uncertainty related to long-term average estimates based on only few data values, only those locations with at least five years coverage for a single month were processed (Table_S 1).We applied ordinary block-kriging interpolation.Therefore, we used respective methods of the 'gstat'-package in R (Pebesma, 2004;Gräler et al., 2016).For ordinary kriging, data-based variogram models are required to describe the structure and degree of spatial dependence of a spatial random variable.In our case, the simplest approach of spatially modelling Temp, Sal, and ln(ChlA) data and their standard deviations would be based on one variogram model for each parameter, valid across all months of the year.To be able to calculate respective empirical variograms, we modified the Ycoordinates of monitoring locations by adding 1500 km for each additional month after January, in order to exclude inter-monthly interference of the values.By visual judgement, we selected the type of underlying variogram function and fit respective models.The received variogram models were validated by leave-one-out cross-validation (Figure_S 2) and also by modelling estimates for those monitoring locations represented by less than five years coverage for individual months (Figure_S 5).
For the long-term means of monthly pelagic habitat factors, we derived predictive powers of R Temp 2 = 0.99, R Sal 2 = 0.97, and R ln(ChlA) 2 -= 0.80.The absolute deviations of predicted vs. observed long-term mean values decrease markedly with higher yearly coverage and reach relatively stable distributions after coverage by values from 4 to 5 years (Figure_S 5).Consequently, it is plausible to limit the spatial modelling of open water conditions to monitoring locations with at least five years coverage.Modelling of respective standard deviations has generally lower predictive powers (R σ(Temp) It is obvious from the respective variograms (Figure_S 3) that there is less spatial dependence related to the data and consequently a stronger random component involved.However, we still observe distinct spatial and temporal patterns of monthly standard deviations in the open water bodies (e.g.Figure_S 6).Therefore, we decided to use the spatially modelled results instead of generalized constant standard deviations over the whole model domain.The final kriging estimation was limited to both the open water selection of the model raster dataset and a convex hull around all monitoring locations of each individual parameter and month with a 20 km buffer.

Separated water bodies
The individual separated water bodies (Fig. 2) were inadequately covered by monitoring stations for spatially resolved modelling.Therefore, we calculated a weighted average of data values from all monitoring locations within these water bodies for each individual month.The respective number of years covered was applied as a weight variable.For each one of the separated water bodies, we calculated individual monthly standard deviations based on the assumption that 'the variance of the pooled set is the mean of the variances plus the variance of the means' (Rudmin, 2010).Thus, data of all monitoring locations within a separated water body was used and weighted with respect to the number of years covered according to the following equation (adapted after Rudmin, 2010): σ: estimated standard deviation for a separated water body; σ 1−n : standard deviation of individual monitoring locations 1 to n; f 1−n : number of years covered for individual monitoring locations 1 to n; ø: weighted average of data values from all monitoring locations within the water body; ø 1−n : average of data values of individual monitoring locations 1 to n.
For some separated water bodies, individual months are only covered by monitoring values from one year.For these months, we estimated respective standard deviations by averaging the standard deviations from all other months.There is one separated water body that was only covered by monthly values from the year 2017 (Kolding Fjord: NOVANA Station No. VEJ0003749).Therefore, no estimate of respective natural variability could be calculated.

Module 2: Blue mussel growth model
This module contains an empirical blue mussel growth model that builds upon the parameterization of the DEB-model.The mussel growth model is based on monthly resolved inputs of Temp, Sal, and ChlA values from Module 1, and it is calibrated against the DEB-model's predictions of biomass for 50 sites from 2008 to 2017.As for the DEBmodel, hydrodynamics were not directly considered here, but competition for food in mussel clusters was indirectly covered.From the parameterization of the calibrated DEB-model, we derived the following ecophysiological response functions to transfer the pelagic habitat factors Temp, Sal, and ChlA to normalized growth performance factors fT, fS, and fC in the range of 0-1 (Figure_S 4): (Arrhenius temperature T A : 5800 K; Lower temperature boundary T L : 275 K; Upper temperature boundary T H : 296 K; Arrhenius temperature for rate of decrease at lower boundary T AL : 45430 K; Arrhenius temperature for rate of decrease at upper boundary T AH : 31376 K; Reference temperature T 1 : 289 K) (Maar et al., 2018;van der Veer et al., 2006.) (Section 2.4; Almada-Villela, 1984;Maar et al., 2015) fC (Section 2.4; Dolmer, 2000;Riisgård et al., 2013b) We further assumed that the actual growth of blue mussels is equally dependent on each of the three growth performance factors.Therefore, we derived an integrated growth performance factor fTSC: From the DEB-Model runs for 50 monitoring stations (Fig. 2) throughout the study area from 2008 to 2017, we extracted average biomasses of individual mussels for mon harv November, to represent potential harvest in autumn, and mon harv March/April, to represent potential harvest in spring; following the abbreviated production period of mitigation culture (Taylor et al., 2019).Due to the overlapped growth period at harvest in spring, we also decided to average the biomasses of the DEB-Model results over a two-month period.Referring to the calculation of monthly Temp, Sal, and ChlA values within the spatial modelling of pelagic habitat factors, monitoring data of all these stations were aggregated as monthly median values of Temp, Sal, and ln(ChlA) within individual years.We transformed these monthly median values to the growth performance factors fT, fS, fC, and finally the integrated factor fTSC by applying Eqs. ( 3)-(6), as described above.Furthermore, we assumed that mussel biomass growth (as dry-weight) is directly related to the integrated monthly fTSC, which an individual mussel has experienced throughout its lifetime after settlement.
After statistically fitting a respective empirical model function, this function was validated against mussel growth data from seven longline test sites in Danish water bodies for growth seasons from 2017 to 2019, where in situ data for Temp, Sal, and ChlA (Taylor et al., 2020) was sufficiently monitored (see above).Validation calculations were performed in the following way: For each growth experiment time-series, we derived a theoretical starting dosage as fTSC 0 from the mussel biomass of the first sampled point by inverting the function in Eq. ( 7).For all further points in the time-series, the time range from start to sample were divided into monthly (30 days) intervals.For each interval we calculated median values of Temp, Sal, and ln(ChlA) and derived corresponding monthly fTSC dosages.To account for incomplete months at the end of each series, we further multiplied the fTSC dosages with the days/30 fraction for each interval.Those monthly fTSC dosages were summed up with each corresponding fTSC 0 and converted to final model biomass (Eq.( 7)).Modelled biomass was evaluated against monitored biomass and deviations compared with the 95% prediction interval of the model function.
Based on data from similar mussel growth experiments published in Taylor et al., 2019 andNielsen et al. (2016)

Module 3: Mussel farm model
Each meter of spat-collector material offers a certain settling area for mussels.When we assume a mussel shape with constant proportions throughout the modelled lifetime, each mussel occupies a part of this settling area (A set ), which is related to its length (l mus ): A set ~lmus 2 .When we further assume that the biomass (m bio ) of a mussel is related to its size, then m bio ~lmus 3 applies.Consequently, the density of individual mussels of the same size on the settling material (ρ mus ½ ind m ) could be described by: The biomass density on the settling material (ρ m bio ½ g DW m ) would then follow the relationship: The biomass of a whole mussel mitigation farm (Harvest [g DW ]) would then become the product of Eqs. ( 1) and ( 9): We used field data on mussel length, biomass, mussel density, and biomass density from longline setups similar to the intended final farm setups (; Taylor et al., 2019;Nielsen et al., 2016) to verify the relationships, assumed above, and to fit respective model functions.In doing so, we assume that settlement of the spat-collector material with mussel larvae is happening under 'ideal' conditions and that no spat limitation is present.This assumption is based on the wide presence of blue mussels in the shallow areas of the western Baltic Sea (Riisgård et al., 2015) and the potential wide dispersion of larvae over a period of up to six weeks (Larsson et al., 2017).We further assume, that individual mussels grow homogeneously under the same pelagic habitat conditions and that these conditions remain similar throughout a mitigation farm.Other external factors that could affect the settling density of mussels on spatcollector material, such as predation by e.g.eider ducks and sea stars, or mussels sloughed by strong currents and physical exposure, are not considered here.
The resulting fit model function for biomass density ρ mbio on the settling material can be used to model biomass contents of a whole farm, based on the total collector length l col (Eq.( 1)) within an individually specified 'standard mitigation farm'.Because of the defined minimum r dep of 2 m, the spatial output of the mussel farm model is limited to those areas with water depths ≥4 m.

Module 4: Avoidance of food limitation model
Significant food limitation within a farm would limit the validity of the above-described mussel growth (Module 2) and mussel farm (Module 3) models, and should consequently be avoided.We assume that within-farm food depletion is dependent on the ChlA-ingestion capacity of mussels within a farm and on the food-flux (ChlA-flux) into the farm (Petersen et al., 2019b).This food-flux is linearly related to hydrodynamics, here treated as the horizontal flow velocity v h , and relates to the residence time of water within the farm t res .Within the farm the water passes the suspended canopy covered in blue mussels.Temp and Sal, are assumed to remain constant throughout the farm, whereas ChlA is assumed to change in relation to the actual ingestion by blue mussels within the farm.In the calibrated DEB-Model, the maximum ChlA-ingestion rate of an individual blue mussel per day (max(in ChlA )) is related to its biomass m bio [g DW ] and the C:ChlA ratio in the phytoplankton (Table_S 2, suppl.Information): As the C:ChlA ratio is naturally highly variable, we applied mean values per month of the year (22 ≤ C:ChlA ≤61) from Jacobsen and Markager (2016) representing the same study area.The theoretical relationship of filtration rate F $ m bio 2 3 was also confirmed in experimental studies that demonstrated exponents very close to 2/3 (Riisgård et al., 2014).
The actual ChlA-ingestion of a mussel is further related to fT and fC: in ChlA mg day Over time, ChlA within the farm is reduced, which again affects fC (Eq.( 5)).The respective reduction rates of ChlA depend on m bio and density of mussels within the farm volume.As discussed above, the density of individual mussels on settling material is assumed to be: . Consequently, the ChlA consumption of each meter spat-collector material Col-in ChlA in the model becomes independent of m bio : The collector density ρ col = l col /V farm , thus, determines volume specific ChlA ingestion rates V-in ChlA : For V farm , we assume that the water column 1 m above and below the collector loops mixes homogeneously within the farm: V farm = A farm * (r dep + 2 m).The spat-collector material and therefore mussels were assumed to be constantly distributed throughout V farm .
We postulate that there is a critical residence time (t res crit ), below which a mussel farm becomes subject to significant within-farm food limitation.In the mussel growth model (Module 2), fC relates linearly to the monthly growth performance factor fTSC of blue mussels (Eq.( 6)).As mentioned above, we assume fT and fS to remain constant.For a given farm setup (r dep , i loop ), Temp, C:ChlA, and starting ChlA 0 , we modelled ChlA(t res n ) and fC(t res n ) over time (see Eq. ( 5)): Resulting v h crit values are then projected on the spatial models of monthly Temp, ChlA, and C:ChlA ratios for a specific farm setup to estimate the spatially explicit required hydrodynamics, which are necessary to maintain predicted mussel growth (Module 2) for this farm setup (Module 3).We compared the resulting values of v h crit with modelled monthly-mean absolute horizontal flow velocities of the BSPAF and FlexSem (Section 2.7).

Module 1: Spatial and temporal habitat factor model
The spatial and temporal modelling of pelagic habitat factors generates geospatial raster files, representing estimates of local long-term mean values of Temp, Sal, and ChlA values, as well as their standard deviations for each month of the year.In Figure_S 6, Temp distribution in July is presented as an example.All underlying maps and spatial datasets can be obtained from the corresponding author upon request.Note that each one of the outlined separated water bodies along the fine-structured coast of the western Baltic Sea (Fig. 2) is represented by constant values in space.Smaller-scale variation is, however, very likely and should be considered locally at higher resolution when looking into detailed implementation.

Module 2: Blue mussel growth model
We found a very significant non-linear relationship between m bio , derived from the DEB-Model, and the monthly-integrated fTSC values (Figure_S 11).A power-function was assumed and we derived the optimal exponent by maximizing the Pearson correlation coefficient over an exponent range from 1 to 4 in intervals of 0.01.Finally, we fit a linear model to predict biomass dry-weight of individual mussels based on total doses of monthly fTSC (Fig. 3).The linear model was forced through zero to avoid predictions of negative mussel biomasses.The The validation delivers reasonable results for all data points acquired from within the first year of the mussel growth season (Fig. 4).Only one of these points falls slightly outside the determined 95% prediction interval of the empirical mussel growth model function (Fig. 3, Eq. ( 18)).Contrarily, most points for months in the 2nd year of a growth season show considerable overestimation of m bio in the model.Note that there is a huge uncertainty related to all four monitored m bio values from January (Figure_S 12).Consequently, the empirical mussel growth  model appears only valid for the first 6 months period of a growth season, even though it is able to reproduce simulations with the DEB-model reasonably up until March/April (Fig. 3).Therefore, all final model outputs shown in this paper are only based on the November harvest scenario.We will subsequently discuss possible explanations for the model overestimation in spring.

Seasonality of blue mussel growth performance
The modelled integrated growth performance factor fTSC is subject to strong seasonality (Figure_S 7) with high values from June until September and a steep drop afterwards until February.Similar plots for the individual growth performance factors fT, fS, and fC are shown in Figures_S 8-S10.

Effects of natural variability in pelagic habitat factors on mussel growth
Each of the three pelagic habitat factors Temp, Sal, and ChlA exhibits natural variability, here represented by the spatial and temporal models of respective standard deviations of the estimated long-term means (e.g.Figure_S 6).The effects of random changes in Temp, Sal and ChlA on the corresponding modelled mussel growth are of non-linear nature (Figure_S 4).Furthermore, cumulated effects of random changes in all three factors require consideration.To evaluate this for the whole model domain, we simulated 500 random scenarios based on the assumed normal distributions of Temp, Sal, and ln(ChlA) in both space and time.We randomly assigned individual multipliers m spg of a normal distribution (R function: rnorm(n = 1,mean = 0,sd = 1)) for each scenario s ∈ {0 : 500}, parameter p ∈ {Temp, Sal, ln(ChlA)}, and growthmonth g ∈ {1 : 9} to derive different realistic growth scenarios for each location.Each m spg was multiplied with its corresponding standard deviation σp g and added to its respective long-term average ∅p g : The resulting values of p sg were then applied in Eq. (2) (fT), Eq. (3) (fS), Eq. (4) (fC), and Eq.(5) (fTSC) to derive modelled m bio (Eq.( 18)) for all 500 scenarios and both harvest times in Nov|MarApr.The process is shown in Fig. 5 for one selected site in the Limfjord.Therein, the red dots illustrate modelled m bio when it is based on the monthly longterm average values ∅p g , only.
As shown in the example for one point only (Fig. 5), it appears that m bio , modelled with the long-term average Temp, Sal, and ChlA conditions, is almost always larger than the calculated median value based on 500 simulated scenarios.We consider that those median values, based on a large number of realistic simulations, more accurately represent the expected long-term average biomass of blue mussels, than  when only taking into account long-term average pelagic habitat factor conditions.Furthermore, we can conveniently use the resulting distributions to derive expectation-ranges of local biomass for individual mussels.All further calculations and visualizations, therefore, are based on these simulated scenarios and we use median values (50-percentile) as well as the 5-95-percentile range of the modelled values (Boxplot-whiskers in Fig. 5).

Spatial model of blue mussel growth
The spatial distribution of modelled biomass of individual blue mussels is shown in Figure_S 13 for harvest in November.The highest 10% of modelled m bio values (0.38-0.71 g DW ) are located in areas within exclusive economic zones of all three countries involved.Key areas are the Limfjord system, the north-western Swedish coast and fjords, the Isefjord and Roskilde Fjord system, as well as a large area extending from Kiel Fjord northwards through the little belt and up to the northern coast of Funen and Horsens Fjord.The expected natural variability of mussel biomass in different growth seasons in these areas is between 0.17 and 0.32 g DW for 99% of the values.

Module 3: Mussel farm model
By empirical analyses of the available field data on mussel length, biomass, mussel density, and biomass density from longline setups (Taylor et al., 2019;Nielsen et al., 2016), we could successfully verify the relationships assumed above (m bio ~lmus 3 ; Eqs. ( 8); ( 9)).Therefore, we fit a respective model function to describe biomass density on the collector-substrate based on the dry-weight of individual mussels (Fig. 6), which resulted in const ρ = 1269.
To estimate the potential range of modelled biomass density, we used the 5th-percentile value of m bio together with the lower boundary of the 95% prediction interval of Eq. ( 9) as the lower end, and the 95thpercentile and upper boundary as the upper end, respectively (Fig. 6).This, we will further refer to as the 'modelled range' of the mussel farm model.The long-term expected middle value is represented by the median m bio in Eq. ( 9).

Spatial distribution of harvest potential for flexible farm setups
As described above, the mitigation farm setup is flexible with respect to r dep and i loop .In Fig. 7, we show the spatial model results for the expected long-term mean harvests for (a) 2-2_0.7 and (b) 2-8_0.7 farm setups.The upper 10% of long-term mean harvests are estimated between 1070-1280 t WW (2-2_0.7)and 3690-4490 t WW (2-8_0.7)per farm.Harvest potentials appear ~3.5 times higher in the second case, which almost represents the respective ratio of maximum l col of 3.55 between both farm setups (Eq.( 1)).Consequently, the bathymetry of shallow near-coastal areas is clearly mirrored in the model results for 2-8_0.7 farms in the spatial distribution.Therefore, the outlines of the areas including the upper 10% of modelled harvests are shifted from shallow to deeper water areas.Maps showing the corresponding modelled ranges are presented in the suppl.Information (Figure_S 14).The key areas are similar to those described for m bio of individual mussels, but areas with b4 m water depth are excluded in the mussel farm model.

Mitigation effect/nutrient extraction potential
In Table 2, a summary of modelled ranges within the upper 10% areas of harvest potentials is shown for different farm setups.According to the model results, the expected long-term mean nutrient extraction

Table 2
Summary of the modelled range for nutrient extraction potentials at different farm setups in the upper 10% areas (Fig. 7).Min, mean and max values refer to the corresponding boundaries of the modelled range (Figure_S 14).potentials per farm at harvest in November are 16.5 t N and 0.94 t P for a 2-2_0.7 setup, as well as 56.7 t N and 3.24 t P for the 2-8_0.7 setup.This corresponds to 0.88 t N /ha and 0.05 t P /ha for a 2-2_0.7 farm and 3.0 t N / ha and 0.17 t P /ha for a 2-8_0.7 farm, respectively.Due to the linear conversion of mussel biomass to nutrient contents in the model, the spatial distribution pattern of potential mitigation effects reflects the one of modelled biomass harvest (Fig. 7).To validate the model results and the above-described derivation of a 'modelled range' based on natural variability and model uncertainty, we compared this modelled range for N-reduction potentials with respective estimates from two field studies (Nielsen et al., 2016;Taylor et al., 2019).These were conducted at three different sites in the Limfjord system (Fig. 8) in the years 2010 and 2017-2018, respectively.Both studies included data based on similar setups to the 2-2_0.7 and 2-2_1.0farms of this model and were harvested in between October and December.We applied the spatial model by using the exact l col values and harvest months as applied in the field experiments.Therefore, the modelled estimates per hectare appear partly higher than those presented before in relation to Table 2.All field estimates of nutrient extraction potential fall well into the modelled range of the modular spatial model.Two field estimates based on 2018 data are close to the upper boundary of the modelled range, while the estimate based on 2010 data falls slightly below the expected long-term middle value.Note, in 2018, the surface water temperatures in Danish fjords and bays were the second highest ever observed, while 2010 was coldest within the last 20 years (Hansen and Høgslund, 2019).

Module 4: Avoidance of food limitation model
After the assumed mussel-and biomass-density functions were validated (const ρ = 1269), the assumption that the volume specific ChlA ingestion V − in ChlA is related to the collector density ρ col (Eq.( 14)) appears valid, too.Therefore, we evaluated the simplified model for critical residence time and flow velocity in order to estimate the required hydrodynamics to avoid food limitation for the modelled mussel growth (Module 2) for specific mussel farm setups (Module 3).
In ) approaches b75% of its starting value fC 0 after t res crit of 210 min.This corresponds to v h crit values of 0.02 and 0.06 m/s for the 250 m and the 750 m axis of the farm, respectively.In Figure_S 15 it can be seen that v h crit is generally decreasing with higher i loop , ChlA 0 , and C:ChlA, while it increases with higher r dep .For Temp, v h crit shows a maximum at around 20 °C, corresponding to Eq. ( 3) and Figure_S 4. The calculated values of t res crit and v h crit can be projected onto the monthly Temp, ChlA, and C:ChlA conditions at given farm setups (i loop , r dep ), resulting in respective spatial distribution models for each month.From July-November (the selected valid model period), v h crit remains relatively constant (Figure_S16).Therefore, we used respective mean values for further evaluation.In Figure_S17, v h crit is shown for the 250 m axis of 2-2_0.7 and 2-8_0.7 farms.In the areas of upper 10% of modelled harvest, the required hydrodynamics expressed as v h crit (250 m) range between 0.02 and 0.10 m/s for the 2-2_0.7 farm and between 0.05 and 0.14 m/s for the 2-8_0.7 farm, respectively.For the long (750 m) axis v h crit is three times higher.For comparing v h crit with modelled values (BSPAF and FlexSem) on v h (EU Copernicus, 2019; Larsen, 2020), we calculated v h ratio = v h /v h crit , which in theory should be N1 to avoid within-farm food limitation.In Table 3, we present selected values of v h ratio from different farm setups and areas, where modelled values on v h were available.The entire upper 10% harvest areas are not continuously covered by modelled v h , in particular not the Limfjord system.There, we used point estimates from FlexSem (Larsen, 2020)  ) for N-reduction potential per hectare with field data from three different locations in the Limfjord system (Nielsen et al., 2016;Taylor et al., 2019).The modelled range for each location, setup, and harvest month is displayed in grey.The modelled and field data both base on setups similar to the 2-2_0.7 and 2-2_1.0farms.Table 3 Average values (Jul-Nov) of v h ratio for the short 250 m axis of different farm setups and locations.For the upper 10% areas (see Fig. 7), we extracted a range covering 90% of the values.The other locations represent mussel growth experiment sites (Fig. 2), where v h data was available.(0.013 m/s).However, based on the underlying modelled values, we derive that v h ratio N 0.5 (minimum for Løgstør farm-scale validation site is 0.44) appears sufficient to avoid significant overestimation of the model due to within-farm food limitation.This condition holds true for 83% of the areas with upper 10% harvest estimates for 2-2_0.7 farms and for 75% for 2-8_0.7 farms, respectively (Fig. 10).The areas with v h ratio b 0.5 are mainly shallow near-coastal waters, fjords and bays.

Validity of the modular spatial modelling approach
The presented modular spatial model for nutrient reduction potentials of blue mussel mitigation farms in the Western Baltic Sea proved valid at all scales and modules.In Module 1, spatial distributions of long-term average values of the pelagic habitat factors Temp, Sal, and ChlA are successfully estimated over the months of a year by geostatistical modelling of operationally acquired marine monitoring data of Denmark (DCE, 2018), Germany (LUNG, 2018;LLUR, 2018), and Sweden (SMHI, 2018).Cross-validation shows that spatial distributions of monthly long-term average Temp, Sal, and ChlA values are estimated with high predictive powers.Module 2 adequately models the growth of individual blue mussels over time by a simplified empirical model function up until harvest in December.This module bases both on data from mussel growth experiments in the field (Taylor et al., 2019) and on runs of a calibrated DEB-model over several years at 50 representative sites throughout the study area.Furthermore, Module 2 provides a reasonable range of year-to-year variability of mussel growth based on the individual variability of pelagic habitat factors.In Module 3, an empirical model function is used to up-scale from individual blue mussels to mitigation farm-scale for various farm setups.The 95% prediction interval of this model function is used on top of the variability estimation of Module 2 to estimate realistic ranges of the model output.These modelled ranges of both biomass harvest and nutrient reduction potentials are in good agreement with results of comparable farm-scale field studies (Fig. 8).Module 4 provides information on the required hydrodynamics, necessary to avoid food limitation of mussel growth (Module 2) for given farm setups (Module 3), based on a simplified modelling approach of ChlA consumption throughout a mitigation farm.Significant food limitation would lead to overestimations in both the mussel growth and mussel farm models.Thus, Module 4 is important to judge the corresponding validity of modelled nutrient mitigation effects.
A study by Bergström et al. (2015) used a statistical multi-parameter modelling approach to attribute parts of the northwestern Swedish coast to three qualitative growth classes for blue mussels.The spatial distribution of these growth classes is generally in very good accordance with the long-term expected middle values of m bio of the mussel growth model (Module 2) at harvest in November (Figure_S 18).Note that considerable parts of the respective compared study area are treated as separated water bodies with constant modelled m bio in this study, while Bergström et al. (2015) could present their results with higher spatial resolution.Kotta et al. (2020) have recently published a similar study on the mitigation potential of mussel farms covering the whole Baltic Sea.In their approach, they use yearly or seasonally averaged modelled environmental input data for estimating mussel growth along with a connectivity model to mussel reefs.Their model results, with connectivity included, deliver zero harvest in many of the upper 10% harvest areas of this study, even though they present a dataset with abundant natural mussel populations in this area, e.g.around Funen and Zealand.This reflects the unsatisfied quality of publicly available regional sediment maps that were used to model populations' connectivities.Without connectivity (data layer personally obtained from J. Kotta) we find high correlation coefficients (Pearson's r = 0.90, Spearman's r = 0.85) and a general agreement between the two models, however, the areas with maximum harvest potential modelled by Kotta et al. (2020), cover almost the whole north-western model domain without any spatial variation.This is because the focus areas of Kotta et al. (2020) were the Baltic Proper and the northeastern basins of the Baltic Sea.They accounted for competition between mussels for space on the longlines, by defining the maximum biomass of mussels on a farm in Kiel Fjord as the maximum biomass yield modelled.Consequently, the model of Kotta et al. (2020) divides large parts of the western Baltic Sea into either highly suitable or non-suitable areas, but does not provide a basis for specific site-selection.Our model, in contrast, offers spatially differentiated estimates along with uncertainty ranges, and is additionally flexible with respect to farm setup and harvest time.Therefore, our model could directly be used for the specific implementation of marine policies (e.g.WFD, MSFD, MSPD) by considering and prioritizing different aspects, respectively targets.An example is demonstrated in the following section.

Model applicability
The spatial model for blue mussel growth and nutrient extraction potential at harvest in November shows sites in the upper 10% of mussel biomass, as well as harvest and nutrient reduction potentials within the exclusive economic zones (EEZ) of all three countries involved .Most of these areas also feature sufficient hydrodynamics to achieve the estimated growth rates.The highest values, however, are modelled for sites in Danish water bodies.With the perspective of this study, the most promising areas for successful and effective mussel mitigation farming in water bodies of Denmark, Germany and Sweden can be selected.
The Danish River Basin Management Plans for the years 2015-2021 (EPA, 2016) define specific N-reduction targets to reach good ecological status in the covered coastal marine water bodies.We made a fictive assumption that these N-reductions could be accomplished by solely applying mussel mitigation farms.Further, we assumed an implementation of only one mitigation farm/km 2 , that all locations with N4 m depth and v h ratio N 0.5 are suitable sites, and that a water body should at least cover 5 km 2 .We compared two farm systems: 2-2_0.7 and 2-8_0.7.Based on the sites with highest modelled values for N-reduction, we calculated for each polygon, how many mussel farms there needed/could be placed to meet the specified N-reduction targets.
In Table 4, we show the top-ten results of this analysis with respect to and sorted after efficiency of the 2-2_0.7 farm setup.It represents those water bodies, where 2-2_0.7 farms would be able to extract at least 100% of the required N-reduction, and where hydrodynamics were sufficient to achieve the modelled N-reduction rates (not considered for the Limfjord system due to lack of continuous model data).Respective locations are shown in Figure_S 20.For Åbenrå Fjord, hydrodynamic conditions were not sufficient (v h ratio b 0.5) for 2-8_0.7 farms.This example illustrates how the results of the modular spatial model for mussel mitigation farms can be applied in implementing specific marine policies.In these water bodies with high mitigation potential (0.9-1.0 t N /(ha*year)), only small areas (b 3.6%) would require to be covered by 2-2_0.7 farms.In this way, decision makers can prioritize areas (e.g. by mitigation requirement) and criteria (e.g.area usage) to identify a suitable set of different setups of mitigation farms, other mitigation measures, and enough space left for additional utilization of marine areas.A recent study (Taylor et al., 2019) has shown ~2 times higher nutrient reduction potentials in the Limfjord with a new technology consisting of floating tubes and nets as collector material.Therefore, mitigation farm efficiency can obviously further be optimized, but the model functions for Modules 2-4 would require re-calibration for the specific conditions within such a farm to provide comprehensive and valid spatial estimates.Other potential in situ mitigation measures, such as sugar kelp farming (0.01-0.03 t N /(ha*year)), translocation of mussels (0.01-0.5 t N /(ha*year)), and eelgrass restoration (0.3 t N / (ha*year)) are estimated less efficient for Danish waters than mussel farming (Bruhn et al., 2020).However, feasibility and efficiency of each mitigation measure are strongly dependent on local habitat conditions.Due to a lack of comparable spatial models, a direct comparison is not yet possible.

Implications for coastal and marine policy implementation and other potential applications
A representative quantitative spatial model of the mitigation potential is a prerequisite for effective implementation of coastal and marine policies with respect to mussel mitigation farming.The presented modular spatial model provides specific quantitative output for a certain location, farm setup, and harvest time.This output is coupled with a modelled range based on natural variability and model uncertainty.Furthermore, locations are identified where significant overestimation of the model is likely due to limited food supply to the farm.These features form a novel quality of data usually used in spatial planning tools for marine aquaculture (Gimpel et al., 2018;Kotta et al., 2020).These are often based on pre-assessed suitability by ecological threshold values for specific aquaculture species (Gimpel et al., 2018), or on single model runs without variability considered (Kotta et al., 2020).The DEB-model based approach of Sainz et al., 2019 considers temporal variability, but is exclusively designed for one specific farm setup and limited by the number of sites and spatial resolution due to the full DEBmodel implementation.Our model approach can be rapidly applied to different farm setups and harvest times at high spatial resolution.Based on several scenarios, specific mitigation concepts can be tested according to intended mitigation targets at high spatial resolution, which is important in the very heterogenically structured western Baltic Sea.
Within the specific use of mussel mitigation farms in integrated multi trophic aquaculture, the habitat requirements of other species need to be considered in addition, as well as potential impacts of e.g.fish farms on ChlA levels (e.g.Petersen et al., 2019a;von Thenen et al., 2020).Therefore, in its current state, this model should only be used with caution in this respect.Potential effects of mussel farming on biodiversity have not been representatively studied, yet and are not included in this model.Mussel farms can affect the food web and act as artificial habitats for other species.The huge additional biomass in the water column can increase biodeposition of faecal matter and shells and attract predators (Petersen et al., 2020).
Other crucial factors for MSP, such as other uses and interests in marine areas, environmental risks to and from mussel mitigation farming, as well as economic considerations are not within the scope of this paper, but will be addressed in a following study.There, we will discuss a multi-criteria GIS tool for optimized site selection of mussel mitigation farms.Resulting layers and methods of this study could also be integrated into already existing MSP products, such as the AquaSpace toolbox (Gimpel et al., 2018).
In general, the modules could in parts be transferred to other aquatic species (marine and freshwater environments), given that ecophysiological or habitat models can adequately describe their respective growth or distribution.The geostatistical approach of Module 1 can theoretically be used to map spatial and temporal distributions of other water quality parameters (e.g.Turbidity, K d , O 2 concentration), given that the respective monitoring density is spatially representative.Consequently, it could contribute to habitat suitability modelling approaches, such as for eelgrass (Staehr et al., 2019).The empirical approach of Module 2 might also work for eco-physiological models of, e.g.other bivalve species and improve the efficiency and performance of large-scale modelling approaches.Harvest potentials of other (mitigation) aquaculture approaches can potentially also be estimated by model functions similar to those described in Module 3. Whenever food-flux can become a limiting factor for farms, the simplified model approach of Module 4 could prove helpful to identify respective critical areas to prevent establishing farms with lower harvest than expected.
In summary, the described approach could prove valuable for implementing environmental policies in aquatic systems involving e.g.nutrient management, in situ mitigation measures, spatial planning, aquaculture in general, as well as habitat suitability mapping of (endangered) aquatic species.

Assumptions and limitations
The model in its current form is based on a number of fundamental assumptions that could partly be addressed for future improvement: Ass. 1. Mussel growth is only dependent on Temp, Sal, and ChlA conditions and all individuals grow identically under the same pelagic habitat conditions.
Mussel growth, survival, and residence on the spat collectors is of course related to many environmental factors, such as e.g.O 2 concentrations, toxic substances (e.g.H 2 S), physical exposure, predators (e.g.starfish and eider ducks), and inorganic resuspended material.Oxygen depletion develops episodically in bottom waters of several areas in the western Baltic Sea and can cause release of toxic H 2 S (e.g.Hansen and Høgslund, 2019).When mussel farms are in contact or proximity to the sediment, anoxia, toxicity and benthic predators can significantly affect harvests (Taylor et al., 2019).A combination of low salinity and high temperatures likely cause loss of mussels from collectors due to weak byssus threads (Buer et al., 2020).The model restricts mussel farms to a minimum distance of 1 m to the sediment to account for sediment related effects.However, overlay of the model with spatial data on the abundance of respective risks could help to select suitable areas more adequately.High loads of inorganic resuspended material can have adverse effects on the ChlA assimilation efficiency.This could be implemented in the DEB-model (Saraiva et al., 2011) and the mussel growth model (Module 2).Coupled with a spatial model on turbidity, the model's applicability could though be extended to high turbid environments.
The currently used mean m bio values of individual mussels might not precisely represent the whole biomass content of a farm, as mussel growth is not constant among individuals.Representative field data on statistical size and biomass distributions in relation to growth could help to further improve the model accordingly.
Ass. 2. Pelagic habitat factors are homogeneously distributed throughout mussel farms.
Studies have shown stratified conditions in suspended mussel culture units, altered within-farm flow conditions (Stevens and Petersen, 2011), and heterogenic patterns of mussel condition across farms (Taylor et al., 2019).Different growth rates can appear on both horizontal and vertical dimensions.Even though in our dataset growth gradients with depth were not observed, under stratified conditions it is important to choose an optimal depth-range for the collectors (Mizuta and Wikfors, 2019).However, increased turbulence within the mussel canopy can also equalize vertical food distribution and compensate for food limitation (Petersen et al., 2019b).A modelling study has shown that differences in growth are considerably more impacted within the horizontal plane (von Thenen et al., 2020).This model's Module 4 intends to avoid the latter.3D numerical ecological models could help to resolve corresponding gradients, but we consider it difficult to integrate this directly into a large-scale spatial model, as these effects are unique for individual farm setups and corresponding fine scale environmental conditions.
Ass. 3. Monthly resolved values of pelagic habitat factors are adequate to describe mussel growth.Indeed, smaller time-intervals may improve the performance of the mussel growth model (Module 2).It is, however, unrealistic to acquire large-scale in situ monitoring data with higher continuous temporal resolution.The alternative use of modelled input data with higher temporal resolution, on the other hand, would introduce additional uncertainties related to the model and increase computing efforts.We also consider other environmental factors, as well as the uncertainty related to upscaling in the mussel farm model (Module 3) to have higher impacts on the estimated harvest and mitigation potentials.Furthermore, the scope of this model is not to estimate individual scenarios precisely, but to estimate realistic long-term expected ranges.The natural variability simulations within Module 2 are intended to account for this.
Ass. 4. Mussel growth on spat-collector substrate starts with ideal settlement all over the spat collector material in July.
Currently, we do not consider limitation of insufficient spat supply and differences in settlement time in the model.This might be reasonable for major parts of the study area (Riisgård et al., 2015), but exceptions are likely.A connectivity model to naturally occurring mussel populations could be integrated, to model settlement intensity and start of the growth season more adequately in space and time.However, with insufficient data quality this step does not necessarily improve model outputs (see discussion of the recent study of Kotta et al., 2020 above).We are currently not aware of large-scale available datasets on mussel populations, respective larvae distribution models, or sitespecific recruitment models with sufficient coverage and quality.
Farms with r dep N 3 m have not been tested in Danish waters, yet.However, similar farming techniques for larger depth-ranges of collector loops are employed in many locations e.g. in the US, Sweden, New Zealand, and Turkey (Mizuta and Wikfors, 2019;Karayücel et al., 2002;Ackerfors and Haamer, 1987).For ideal settlement, it is concluded (Mizuta and Wikfors, 2019;Karayücel et al., 2002) that the collector material should be placed in the depth range with optimal habitat conditions for the mussels.We consider this technologically feasible in the field, by e.g.modifying pre-and post-recruitment substrate depths, and adoption of alternative buoyancy methods to accommodate higher biomass loads.
Ass. 5. Mussels contain constant fractions of nitrogen and phosphorus.
The assumed constant 14% N and 0.8% P contents of biomass dryweight are in reality subject to variability (e.g.Kotta et al., 2020;Smaal and Vonck, 1997;Taylor et al., 2019).Differences in N and P content likely appear in relation to environmental factors (e.g.Sal, ChlA) and seasonality (e.g.spawning) and not on a random statistical basis.To our knowledge, respective functional relationships that could be applied in our model are not yet comprehensibly described in the literature.
Ass. 6. Mussels are homogeneously distributed in size and density throughout the farm volume.
This assumption is used in Module 4, but in reality, aggregated mussel biomass and mussel condition can cluster in a farm on at least two scales.(1) Mussel aggregates may be discontinuous in density along the collector length; (2) within the canopy, heterogeneous clusters can form because of food limitation or husbandry practice (e.g.Taylor et al., 2019).The effect of mussel clusters on small-scale ChlA depletion was described by Nielsen et al. (2016) and Petersen et al. (2019b).However, this variability of aggregate distribution on collector material is contained within the model calibration, as all mussel growth data from the field were acquired on longline systems similar to the one modelled here.Small-scale ChlA depletion gradients integrating with meso-scale depletion reduces the food-uptake efficiency across a farm, but is complex in relation to hydrodynamics and food-filtration feedback mechanisms, and is likely part of the explanation why model results were still valid at v h ratio b 1. Representative field investigations of in situ flow conditions and related ChlA depletion at a large farm-scale would help to specify the threshold of v h ratio more precisely for adequate site selection.

Limitation of the model to account for biomass loss over winter
The empirical mussel growth model overestimated biomass for those mussel samples derived from the period of the growth season from January until May (Fig. 4).In the field data, average m bio of individual mussels drops considerably over winter.Biomass loss, however, cannot be reflected by the empirical model function derived here, as it integrates monthly values between 0 and 1. DEB-model results, however, do not show such considerable biomass drops over the winter months (data not shown), even though it takes biomass loss by respiration and starvation into account.
Other external factors are, therefore, likely to have affected average m bio in the field experiments.There is a very strong peak of m bio standard deviations in January (Figure_S 12).Riisgård et al. (2015) report a typical second peak of mussel larvae in the autumn within the model domain of this study.A second recruitment of mussel larvae in the autumn could cause lower average mussel biomass in the following months if there is significant self-thinning (Lachance-Bernard et al., 2010).Physical exposure, predation and exposure to anoxic/toxic conditions after sediment contact over winter (Taylor et al., 2019) could have affected m bio distributions on test lines.Nevertheless, after thorough analysis of field data from mitigation farm tests in the Limfjord system, Taylor et al., 2019 recommend an early harvest in November/ December for Danish waters and this is the period, where the growth model presented here proved valid.

Conclusion and perspectives
In this study, we developed a modular spatial model for nutrient mitigation potential of blue mussel farms in the western Baltic Sea.Unlike other available approaches, this model includes uncertainty estimation, takes into account temporal variability of habitat conditions, and is flexible with respect to farm setup and harvest time.Furthermore, required hydrodynamics to avoid food limitation are part of the model results.For Danish, German, and Swedish water bodies, we identified key areas for efficient mussel mitigation farms.In efficient Danish areas, mitigation farms (18.8 ha, 90 km collector substrate in loops with 2 m depth-range) required b3.6% of the space to extract the target nitrogen loads for good ecological status.The developed approach is a highly suitable tool for the detailed planning of mussel mitigation farms, e.g.within coastal and marine policy implementation and MSP.Furthermore, the methodology can be transferred to other aquatic species (marine and freshwater environments) to model e.g.aquaculture productivity and habitat suitability of (endangered) species.
In a following paper, we will further implement the specific results of this study into a multi-criteria tool for optimized site selection of mussel mitigation farms in the western Baltic Sea.This will include spatial data on other uses and interests in marine areas, environmental risks to and from mussel mitigation farming, as well as economic considerations.

Fig. 1 .
Fig.1.Flow diagram of the spatial modular modelling approach to estimate harvest and nutrient reduction potentials and locally required hydrodynamics for Mussel Mitigation Farms on trans-national scale in the western Baltic Sea.Black: inputs to the spatial model; green: modules of the spatial model; blue: outputs of the spatial model.

Fig. 2 .
Fig. 2. Overview map of the study area showing the national exclusive economic zones (EEZ), model domain, outlined separated water bodies, monitoring sites used in separated water bodies (SW) and open sea (OS), sites where the DEB-model was run, and sites of in situ longline field experiments.
we derived average constant factors to transfer mussel biomass dry-weight m bio [g DW ] to biomass wet-weight m bio [g WW ], as well as nitrogen m N [g N ] and phosphorus m P [g P ] content of the whole mussel including tissue, shell and byssus threads: m Bio [g DW ] = m Bio [g WW ]/9.68 = m N [g N ]/0.14 = m P [g P ]/ 0.0080.These factors are used to estimate local potentials of nutrient reduction by mussels.
Therein, we modelled all possible combinations of r dep = c(2 : 8); i loop = c(0.7,1.0, 1.5, 2.0); Temp = seq(0 : 25, by = 0.5); C : ChlA Jul Mar =c(61, 53,43,30,24, 23, 22, 23,32); ChlA 0 = c(0.21* 1.1 0:65 ) (R functions).Time intervals for the steps 1:n were continuously increased over the iteration from 60 s to n * 60 s.We defined t res crit as that specific residence time, after which the following criterion was not met: mean(fC 0t res)/fC 0 N 0.75.This means that as long as t res b t res crit , the expected overall growth performance of the farm (fTSC = fT * fS * mean(fC 0 tres )) for a given month would be at least 75% of the one modelled for constant ChlA conditions (fTSC = fT * fS * fC 0 ) (Module 3).From t res crit , we can in turn derive the required hydrodynamics, as the critical minimum flow velocity [m/s], for the flow distances d along both the short (250 m) and the long (750 m) axis of the farm, necessary to avoid significant within-farm food depletion: 95% prediction interval was almost constantly ±0.16 g DW across all biomasses predicted by the DEB-model.

Fig. 3 .
Fig. 3. Linear fit of the monthly-integrated growth performance factors fTSC vs. biomasses of DEB-model runs for two selected harvest times in November and March/April.

Fig. 4 .
Fig. 4. Validation of the Blue Mussel Growth Model with data from seven longline test sites in Danish water bodies for growth seasons from 2017 to 2019 (Taylor et al., 2020).Stars with error bars indicate data and standard deviations from 2nd year of a growth season and labels show corresponding months.

Fig. 5 .
Fig. 5. Exemplary illustration of modelled blue mussel growth for one selected site in the Limfjord, taking the effect of natural variability of pelagic habitat factors into account.(left) Monthly long-term average pelagic habitat factors with their respective monthly standard deviations.(right) Boxplots of the modelled distributions of mussel biomass m bio are based on 500 random scenarios of Temp, Sal, and ChlA.The red dots illustrate modelled m bio based on the monthly long-term average values, only.

Fig. 6 .
Fig. 6.Derivation of fit model functions describing length vs. dry-weight relationship (a), as well as mussel density (b) and biomass density (c) vs. dry-weight on collector substrate in a blue mussel mitigation farm.The latter is used for upscaling from individual mussel dry-weight to biomass content on the farm-scale.Data points with exceptionally low mussel densities (stars) were excluded to only fit the mussel farm model to ideal settling conditions.Dashed grey lines represent the 95% prediction interval of the model functions (solid lines).Derivation of the modelled range of biomass density is illustrated in green.

Fig. 7 .
Fig. 7. Example of spatial distributions of expected long-term mean harvest and nutrient reduction potentials of standard model mitigation farms at harvest in November.(a) 2-2_0.7 farm; (b) 2-8_0.7 farm.The areas of the model domain including the upper 10% of values are highlighted by green polygons.Borders of exclusive economic zones (EEZ) are drawn in black.
Fig. 9, the derivation of t res crit is visualized for Temp of 20 °C, ChlA 0 of 5 μg/L, i loop of 1 m, r dep of 2 m and the July C:ChlA value of 61.The average of the modelled growth performance factor mean(fC 0 tres Fig. 8.Comparison of results of the mussel farm model (Module 3) for N-reduction potential per hectare with field data from three different locations in the Limfjord system(Nielsen et al., 2016;Taylor et al., 2019).The modelled range for each location, setup, and harvest month is displayed in grey.The modelled and field data both base on setups similar to the 2-2_0.7 and 2-2_1.0farms.

Fig. 9 .
Fig. 9. Visualization of the simplified model to derive critical residence time t res crit , above which significant within-farm food depletion becomes likely.(green) ChlA is reduced over time due to ingestion by mussels; (blue) fC over time, as related to ChlA (Eq.(5)); (orange) mean(fC 0 tres )/fC 0 over time.The derivation of t res crit is illustrated where the orange line falls below the threshold of 0.75.

Fig. 10 .
Fig. 10.Spatial distribution maps of v h ratio for the short (250 m) axis of (a) 2-2_0.7 and (b) 2-8_0.7 farm setups.v h ratio is the fraction of modelled horizontal flow velocities v h (BSPAF, EU Copernicus, 2019) vs. modelled critical flow velocities v h crit .The areas of the model domain including the upper 10% of values are highlighted by red polygons, while therein areas with v h ratio N 0.5 are highlighted by green polygons.Validation with farm-scale field data shows that v h ratio N 0.5 is sufficient to avoid significant overestimation of the model due to withinfarm food limitation.Note: Several fjords (including the Limfjord system) not covered by continuous v h from BSPAF.
Module 3 contains a statistical model function to up-scale biomass of individual blue mussel to the farm-scale and delivers respective biomass harvest and nutrient reduction potentials.This function is based on data of blue mussel growth experiments in the field and describes different types of mitigation farm setups in relation to local bathymetry.
• Module 2 contains an empirical blue mussel growth model based on monthly resolved inputs of Temp, Sal, and ChlA values (Module 1).This model was statistically fit to results of eco-physiological DEBmodel runs distributed throughout the model domain and validated with data from blue mussel growth experiments in the field.This module is used to estimate biomass of individual blue mussels at harvest in November throughout the whole model domain.•

Table 1
Summary of the input and output parameters used and produced by each of the four modules.

Table 4
Example application of the spatial nutrient extraction model for mussel mitigation farms on Danish water bodies.2-2_0.7 and 2-8_0.7 farm setups are applied.The selection is limited to the top-ten efficiency results where 2-2_0.7 farms could achieve N100% of required N-reduction and where hydrodynamic conditions were sufficient (not considered for Limfjord and Skive Fjord).The numbers represent the modelled long-term mean values of N-reduction.Corresponding locations are shown in Figure_S 20.
a Water body names and IDs refer to the Danish River Basin Management Plans for the years 2015-2021 (EPA, 2016).b Hydrodynamics not considered.