Pronounced Seasonal and Spatial Variability in Determinants of Phytoplankton Biomass Dynamics along a Near–Offshore Gradient in the Southern North Sea

Abstract

Marine phytoplankton biomass dynamics are affected by eutrophication, ocean warming, and ocean acidification. These changing abiotic conditions may impact phytoplankton biomass and its spatiotemporal dynamics. In this study, we used a nutrient–phytoplankton–zooplankton (NPZ) model to quantify the relative importance of the bottom-up and top-down determinants of phytoplankton biomass dynamics in the Belgian part of the North Sea (BPNS). Using four years (2014–2017) of monthly observations of nutrients, solar irradiance, sea surface temperature, chlorophyll-a, and zooplankton biomass at ten locations, we disentangled the monthly, seasonal, and yearly variation in phytoplankton biomass dynamics. To quantify how the relative importance of determinants changed along a near–offshore gradient, the analysis was performed for three spatial regions, i.e., the nearshore region (<10 km to the coastline), the midshore region (10–30 km), and the offshore region (>30 km). We found that, from year 2014 to 2017, the phytoplankton biomass dynamics ranged from 1.4 to 23.1 mg Chla m⁻³. Phytoplankton biomass dynamics follow a general seasonal cycle, as is the case in other temperate regional seas, with a distinct spring bloom (5.3–23.1 mg Chla m⁻³) and a modest autumn bloom (2.9–5.4 mg Chla m⁻³). This classic bimodal bloom pattern was not observed between 2003 and 2010 in the BPNS. The seasonal pattern was most expressed in the nearshore region. The relative contribution of factors determining phytoplankton biomass dynamics varied spatially and temporally. Throughout a calendar year, solar irradiance and zooplankton grazing were the most influential determinants in all regions, i.e., they jointly explained 38–65% of the variation in the offshore region, 45–71% in the midshore region, and 56–77% in the nearshore region. In the near- and midshore regions, nutrients were the greatest limit on phytoplankton production in the month following the spring bloom (44–55%). Nutrients were a determinant throughout the year in the offshore region (27–62%). During winter, sea surface temperature was a determinant in all regions (15–17%). By the high-resolution spatiotemporal analysis of the relative contributions of different determinants, this study contributes to a better mechanistic understanding of the spatiotemporal dynamics of phytoplankton biomass in the southern North Sea. This detailed understanding is anticipated to contribute to the definition of targeted management strategies for the BPNS and to support sustainable development in Belgium’s blue economy.

1. Introduction

Marine phytoplankton, which forms the base of marine food webs, is responsible for about half of the world’s total primary production [1]. Net global marine primary production is estimated at 50.7 Gt carbon per year [2]. Principal factors that determine marine primary production are solar irradiance, nutrient availability, and sea surface temperature (SST), and these factors mainly limit phytoplankton’s growth rate and carrying capacity. According to Liebig’s law of the minimum, phytoplankton production will be as high as allowed by the least available resource [3]. However, in many cases, co-limitation by resources is a better description of the factors that influence phytoplankton biomass dynamics [4,5]. Living in the Anthropocene, phytoplankton biomass dynamics are affected by human activities that directly or indirectly alter the abiotic marine environment, such as the burning of fossil fuels, eutrophication, and chemical pollution. To date, we have limited insight into how the combination of changing conditions may affect phytoplankton biomass dynamics at high-resolution spatiotemporal scales.

In temperate marine regions, phytoplankton biomass dynamics follow an annual cycle consisting of spring and autumn phytoplankton blooms, followed by periods of zooplankton grazing. Phytoplankton blooms are triggered by high nutrient availability and sufficient solar irradiance [6]. After a few weeks of rapid growth, phytoplankton biomass becomes restricted by nutrient limitation and zooplankton grazing. As in other parts of the North Sea, the most important factors that determine the phytoplankton biomass in the Belgian part of the North Sea (BPNS; 3454 km² [7]) are nutrient concentrations, SST, and solar irradiance [8,9,10,11]. Everaert et al. [12] is one of the first studies that quantified the relative importance of these conditions in the BPNS. Based on a relatively short time series at one location in the BPNS, it was found that SST and solar irradiance accounted for 20% (summer) to 50% (winter) of the observed seasonal variation [12] and can thus be considered potentially key determinants. Nutrients appeared to be less of a determinant than SST and solar irradiance in the BPNS [12], and this was also found to be the case for the entire North Sea by Llope et al. [13] and McQuatters-Gollop et al. [14]. Nutrients become the dominant determinant of phytoplankton biomass (30%) in the month after the phytoplankton bloom. Besides these bottom-up determinants, there is also a strong top-down control of phytoplankton biomass dynamics by zooplankton grazing, i.e., up to 50% of the phytoplankton growth limitation [12]. However, the BPNS is a heterogeneous and dynamic coastal area, so it is doubtful whether the quantifications in Everaert et al. [12] are generalisable for the entire BPNS. The BPNS is relatively shallow with water depths gradually increasing to 45 m from the southeast towards the northwest [15]. Sea surface temperatures vary seasonally between 5 °C and 20 °C. The salinity is strongly influenced by the river plumes of the Scheldt, Rhine, Seine, and Meuse [16], and it varies between 29 to 35 PSU. Seawater from the English Channel, which contains runoff from the Seine, flows northwards through the BPNS driven by the anti-clockwise current in the North Sea [17]. In the case of nutrients, the Seine plays a major role, except in the vicinity of the Scheldt estuary and in the northern part of the BPNS, i.e., it is mainly influenced by the influx of water from the Atlantic Ocean [18]. Due to its shallowness and riverine influx, there is a permanent vertical mixing of the water column in the BPNS [16]. Overall, in heterogeneous and dynamic coastal areas such as the BPNS, the relative contribution—i.e., how much a determinant contributes to the change in phytoplankton biomass in relation to the cumulative contribution of all determinants—of bottom-up and top-down determinants may be likewise dynamic, both spatially and temporally. The BPNS being a prime example of such a system in combination with the availability of long-term observations with high spatial resolution offers unique possibilities to study the scales at which the relative contributions of bottom-up and top-down determinants may shift. Having a better understanding of the spatial variation in the relative contribution of the determinants of phytoplankton biomass dynamics can lead to a more adjusted and specified management of the BPNS, as well as to the further development of the blue economy in Belgium.

In the present research, the aim is to analyse which factors drive marine phytoplankton biomass dynamics in the BPNS and how their relationship to primary production varies on a spatiotemporal scale. In particular, we analysed how the relative contributions of SST, nutrient regimes, solar irradiance, and zooplankton grazing to the marine phytoplankton biomass change, both spatially and temporally. We used a nutrient–phytoplankton–zooplankton (NPZ) model (obtained from Soetaert and Herman [19] and adjusted by Everaert et al. [12]) to simulate changes in plankton density in the nearshore, midshore, and offshore regions; this simulation was based on the monthly data collected from 2014 to 2017 at ten sampling locations in the BPNS. The novelty of this study lies in its high-resolution spatiotemporal analysis, which examines the varying impacts of potential determinants on phytoplankton biomass. The high spatiotemporal analysis contributes to a better understanding of phytoplankton biomass dynamics in the BPNS. It is anticipated that this comprehensive exploration has the potential to guide more focused management strategies for the BPNS and could bolster sustainable development within Belgium’s blue economy.

2. Materials and Methods

2.1. Input Data

Three regions of interest were studied (Figure 1), i.e., the nearshore region (10 km), the midshore region (10–30 km), and the offshore region (>30 km). These regions were defined based on an integration of information about their distance to the coast, their sediment composition, bathymetry [20,21], and prior knowledge about different abiotic conditions [22]. LifeWatch measuring stations are located in each of these regions (indicated as triangles in Figure 1). LifeWatch is a European research infrastructure within the European Strategy Forum on Research (ESFRI) that focuses on biodiversity research and activities, such as measuring the biotic and abiotic environment. LifeWatch stations in the near- and midshore are visited monthly [23]. We used data from two stations, i.e., 130 for nearshore and 330 for midshore, as well as pooled the data of seven sampling stations for the offshore region (Figure 1). The data of all offshore stations, i.e., seven stations, were pooled in this study due to the lower temporal sampling. In the offshore region, the LifeWatch stations are visited seasonally, whereas stations in the near- and midshore are visited monthly.

Four open-access datasets obtained from LifeWatch sampling campaigns, which are coordinated by the Flanders Marine Institute, at nine locations in the BPNS from 2014 to 2017 were used [23,24,25,26]. A first open-access dataset related to the nutrient concentrations, i.e., NH₄, NO₃, NO₂, PO₄, and SiO₃, was used and the measurements were performed by a Skalar AutoAnalyser system [23,25]. A second open-access dataset consisted of seawater temperature measurements, and this was performed by a CTD [25]. Environmental measurements were taken at a depth of three metres. A third open-access dataset comprised zooplankton abundances, which were obtained through ZooScan analysis [24,26]. A fourth and final dataset used in this research contained information on the in situ pigment concentrations, i.e., chlorophyll-a, and these measurements were performed by HPLC [23,25]. The pigment concentrations were the result of pigment purification and sample analysis, i.e., obtained by taking the cumulative sum of phytoplankton taxa found in the water column. By doing so, the model takes into account the presence of phytoplankton taxa in the BPNS.

The gathered data were assembled in two new datasets: the first containing the data related to nutrient and pigment concentrations, as well as seawater temperature; the second containing the zooplankton abundances. The zooplankton dataset contains the abundances of the following taxa, as they are the most abundant in the BPNS (Appendix C [27]): Calanoida, Noctiluca, Harpacticoida, and Appendicularia. All data were timestamped and were location-specific.

In addition to the SST dataset gathered from LifeWatch, a second SST dataset with a higher temporal coverage was required to infer the daily time series, i.e., the input data for our nutrient–phytoplankton–zooplankton model, which was achieved by means of generalised additive models. The SST data from the Westhinder station, i.e., the data from the Westhinder measuring pile, was complemented with data from the Westhinder–buoy (2% of the dataset). This station is part of the Flemish Banks Monitoring Network [28], and their SST data was used to infer the daily time series for nutrients and the SST. Due to large data gaps in the SST dataset of the Westhinder station in 2018, we have chosen to use the dataset from 2014–2017 in order to avoid increasing random noise when calibrating the model by including an extra dataset of the SST.

2.2. Time Trends for Input Data

Generalised additive models (GAM) were used to infer the daily trends for nutrient concentrations, i.e., N, P, and Si, in the three regions of interest (Appendix C). These daily trends were based on the monthly and seasonal data collected via LifeWatch. The daily time series were used as input for the nutrient–phytoplankton–zooplankton model. Dissolved inorganic nitrogen (DIN) was calculated as the sum of ammonium (NH₄⁺), nitrate (NO₃⁻), and nitrite (NO₂⁻). Measurements of NH₄, NO₃, NO₂, PO₄, and SiO₃ registered in the LifeWatch database were already converted into nitrogen, phosphorus, and silica equivalent weights, respectively. The covariates used in the GAM models were month and year, as in Everaert et al. [12], to infer the daily trends. The model fit of the GAMs was assessed following the method of Zuur et al. [29], i.e., using the Akaike Information Criterion (AIC) and inspecting the homogeneity and normality of the model residuals. The AIC was calculated via the package ‘stats’ [30], and the minimum AIC was used to select the best-fit GAM models (

Table A7. Performance of the different generalised additive models for each variable of interest, i.e., phosphate (PO₄), ammonium (NH₄), nitrite (NO₂), nitrate (NO₃), silicate (SiO₃), and the sea surface temperature (SST), in the nearshore, midshore, and offshore regions using data from 2011 to 2017. The number of basis functions (k) to use for each smooth term (s(day), s(year)) are included. The Akaike information criterion (AIC) and the penalised R² for adding variables to the model (adjusted R²), i.e., a p-value whether the function s() (smooth term) is non-zero, performs well (p < 0.05). The model formulation (the GAMs) is also provided. Models in bold were selected for modelling the corresponding variable.

Nearshore Region
Variable	k-	k-	AIC	Adjusted R2	k Performance	GAM
	s (day)	s (year)
PO4	3	3	58.1	0.19	p-value < 0.05 for day, and 0.38 for year	PO4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	50.26	0.30	p-value < 0.05 for day, and 0.23 for year
	5	5	46.96	0.35	p-value < 0.05 for day, and 0.07 for year
	6	6	47.6	0.36	p-value < 0.05 for day, and 0.06 for year
NH4	3	3	203.54	0.04	p-values > 0.05 for both smoothers	NH4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	201.27	0.12	p-values > 0.05 for both smoothers
	5	5	202.06	0.11	p-values > 0.05 for both smoothers
NO2	3	3	44.04	0.08	p-value < 0.05 for day, and 0.86 for year	NO2~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	33.44	0.25	p-value < 0.05 for day, and 0.79 for year
	5	5	34.31	0.25	p-value < 0.05 for day, and 0.81 for year
NO3	3	3	336.38	0.25	p-values < 0.05 for day, and 0.56 for year	NO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	334.55	0.31	p-values < 0.05 for day and 0.26 for year
	5	5	335.36	0.32	p-values < 0.05 for day, and 0.24 for year.
SiO3	3	3	481.08	0.28	p-value < 0.05 for day, and 0.22 for year	SiO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	476.2	0.35	p-value < 0.05 for both smoothers
	5	5	471.5	0.41	p-value < 0.05 for both smoothers
	6	6	472.28	0.41	p-value < 0.05 for both smoothers
Midshore Region
Variable	k-	k-	AIC	Adjusted R2	k Performance	GAM
	s(day)	s(year)
PO4	3	3	10.82	0.21	p-values < 0.05 for day, and 0.67 for year	PO4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	4.53	0.30	p-values < 0.05 for day, and 0.71 for year
	5	5	4.96	0.31	p-values < 0.05 for day, and 0.69 for year
NH4	3	3	136.188	0.09	p-values > 0.05 for both smoothers	NH4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	136.19	0.09	p-values > 0.05 for both smoothers
	5	5	136.188	0.09	p-values > 0.05 for both smoothers
NO2	3	3	−27.8	0.37	p-values < 0.05 for day, and 0.68 for year	NO2~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	−26.36	0.37	p-values < 0.05 for day, and 0.72 for year
	5	5	−25.81	0.37	p-values < 0.05 for day, and 0.72 for year
NO3	3	3	271.38	0.37	p-values < 0.05 for day, and 0.31 for year	NO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	273.25	0.38	p-values < 0.05 for day, and 0.41 for year
	5	5	274.05	0.38	p-values < 0.05 for day, and 0.39 for year
SiO3	3	3	362.68	0.28	p-value < 0.05 for day, and 0.54 for year	SiO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	363.23	0.28	p-value < 0.05 for day, and 0.56 for year
	5	5	363.47	0.28	p-value < 0.05 for day, and 0.55 for year
Offshore Region
Variable	k-	k-	AIC	Adjusted R2	k Performance	GAM
	s(day)	s(year)
PO4	3	3	−281.69	0.80	p-value < 0.05 for both smoothers	PO4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	−288.22	0.82	p-values < 0.05 for both smoothers
	5	5	−306.43	0.84	p-values < 0.05 for both smoothers
	6	6	−327.19	0.87	p-values < 0.05 for both smoothers
NH4	3	3	232.78	0.05	p-values ≥ 0.05 for both smoothers	NH4~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	230.73	0.09	p-values > 0.05 for both smoothers
NO2	3	3	−46.83	0.37	p-values < 0.05 for both smothers	NO2~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	−72.92	0.51	p-values < 0.05 for both smoothers
	5	5	−92.81	0.60	p-values < 0.05 for both smoothers
	6	6	−125.85	0.70	p-values < 0.05 for both smoothers
NO3	3	3	430.44	0.62	p-values < 0.05 for day, and 0.10 for year	NO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	431.57	0.62	p-values < 0.05 for day, and 0.14 for year
SiO3	3	3	401.81	0.26	p-values < 0.05 for day, and 0.07 for year	SiO3~s(day, k = ks(day)) + s(year, k = ks(year))
	4	4	399.83	0.28	p-value < 0.05 for day, and 0.37 for year
	5	5	396.58	0.31	p-value < 0.05 for day, and 0.47 for year
	6	6	396.99	0.31	p-value < 0.05 for day, and 0.47 for year

). The corresponding R² of the best-fit GAM model provides an indication of how well the model fits the observational data (

Table A6. Performance of the generalised additive models in creating the daily time series of the nutrient data, i.e., dissolved inorganic nitrogen (DIN), phosphate (PO₄), silicate (SiO₃), and the sea surface temperature (SST) data in the three regions of interest. The root mean square error (RMSE) and the coefficient of determination (R²) were calculated for each nutrient and the SST against field observations per region.

Region of Interest	Nutrient and SST	RMSE	R2
Nearshore region	DIN (mmol N m−3)	9.61	0.30
	PO4 (mmol P m−3)	0.30	0.32
	SiO3 (mmol Si m−3)	6.43	0.39
	SST (°C)	1.82	0.94
Midshore region	DIN (mmol N m−3)	5.69	0.37
	PO4 (mmol P m−3)	0.23	0.30
	SiO3 (mmol Si m−3)	3.65	0.29
	SST (°C)	0.88	0.96
Offshore region	DIN (mmol N m−3)	2.33	0.55
	PO4 (mmol P m−3)	0.05	0.87
	SiO3 (mmol Si m−3)	1.30	0.30

).

The daily SST for the near- and midshore stations were calculated using GAM models that were based on the temperature observations from LifeWatch. For the offshore region, the SST data of the Westhinder station (Figure 1) were used. In the case that multiple temperature loggings were available per day, the median daily SST was used as the input for the model.

2.3. Ecological Model

A nutrient–phytoplankton–zooplankton (NPZ) model as presented in Everaert et al. [12] (Figure 2), which was adjusted from Soetaert and Herman [19], was used to simulate spatiotemporal changes in plankton density in the BPNS from 2014 to 2017. The model time step is expressed in mmol N m⁻³ d⁻¹, and the state variables are expressed in mmol N m⁻³. The necessary input data for the NPZ model are (i) nutrient concentrations, (ii) SST, and (iii) solar irradiance (Photosynthetically Active Radiation (PAR)), which is corrected for the diffusion attenuation coefficient (the K_d is calculated at a three-metre depth; Equation (10)). The PAR and nutrients, i.e., DIN, PO₄, and SiO₃, were implemented in the NPZ model as saturating Michaelis–Menten equations, as shown in Arndt et al. [8]. As such, the Michaelis–Menten equations describe the determining factor of a variable of interest, i.e., PAR, DIN, PO₄, or SiO₃, for each time step. The outcome of these equations varies from one to zero, with a value of one indicating no limitation for plankton growth, and a value of zero indicating a complete limitation of the plankton growth [19]. Both the vertical and horizontal mixing of nutrients were incorporated in the model as they were inherently included in the in situ observations used to calibrate and fit our model. The amount of PAR available was corrected for the diffuse attenuation coefficient by means of the Lambert–Beer law [31,32]. The same fixed sinusoidal surface irradiance for all three regions was used as the angle of the sunlight and surface, and light hours were assumed to be the same across all regions in our small study area. The influence of the SST on plankton growth followed the Thomann and Mueller [33] equation (Equation (6)) [12]. The NPZ model and corresponding calculations were performed in R [30] (version 3.4.4; R packages ‘doParallel’ [34], ‘dplyr’ [35], ‘foreach’ [36], ‘ggplot2′ [37], ‘ggpubr’ [38], ‘lubridate’ [39], ‘parallel’ [30], ‘plyr’ [40], ‘RColorBrewer’ [41], ‘reshape2′ [42], ‘stats’ [30], ‘viridis’ [43], and ‘xts’ [44]).

The equation-based description of the model is presented following the state variables of the model (set according to Soetaert and Herman [19]). Nutrients were depicted as regression-based models, i.e., for nitrogen GAM_DIN, phosphorus GAM_PO4, and silica GAM_SiO3 (Equation (1)).

$$\mathrm{N}\mathrm{u}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}={\mathrm{G}\mathrm{A}\mathrm{M}}_{\mathrm{D}\mathrm{I}\mathrm{N}}+{\mathrm{G}\mathrm{A}\mathrm{M}}_{{\mathrm{P}\mathrm{O}}_4}+{\mathrm{G}\mathrm{A}\mathrm{M}}_{{\mathrm{S}\mathrm{i}\mathrm{O}}_3}$$

(1)

The differential equations determined the rates of change in the biomass of phytoplankton (dPHYTO; Equation (2)) and zooplankton (dZOO; Equation (3)).

$$\frac{\mathrm{d}\mathrm{P}\mathrm{H}\mathrm{Y}\mathrm{T}\mathrm{O}}{\mathrm{d}\mathrm{t}}={\mathrm{N}}_{\mathrm{u}\mathrm{p}\mathrm{t}\mathrm{a}\mathrm{k}\mathrm{e}}-\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}$$

(2)

$$\frac{\mathrm{d}\mathrm{Z}\mathrm{O}\mathrm{O}}{\mathrm{d}\mathrm{t}}=\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}-\mathrm{F}\mathrm{a}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{s}-\mathrm{E}\mathrm{x}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}-\mathrm{M}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}$$

(3)

where t indicates time in days; N_uptake, the phytoplankton’s nutrient uptake; Grazing, the rate of grazing of zooplankton on phytoplankton; Faeces, the faeces production of zooplankton; Excretion, the zooplankton excretion; and Mortality, the mortality of zooplankton. N_uptake is the phytoplankton’s nutrient uptake, which was expressed as the maximum uptake of nutrients (maxUptake) and is limited by the following: solar irradiation (Par_lim), SST (Temp_lim), the concentrations of nitrogen (DIN_lim), phosphorus (P_lim), and silica (Si_lim) expressed as limitation factors, and the phytoplankton biomass (PHYTO; Equation (4)). The equations used to define the limitation factors follow the saturating Michaelis–Menten equations [19]:

$${\mathrm{N}}_{\mathrm{u}\mathrm{p}\mathrm{t}\mathrm{a}\mathrm{k}\mathrm{e}}=\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{U}\mathrm{p}\mathrm{t}\mathrm{a}\mathrm{k}\mathrm{e}·{\mathrm{P}\mathrm{a}\mathrm{r}}_{\mathrm{l}\mathrm{i}\mathrm{m}}·{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{i}\mathrm{m}}·{\mathrm{P}}_{\mathrm{l}\mathrm{i}\mathrm{m}}·{\mathrm{D}\mathrm{I}\mathrm{N}}_{\mathrm{l}\mathrm{i}\mathrm{m}}·{\mathrm{S}\mathrm{i}}_{\mathrm{l}\mathrm{i}\mathrm{m}}·\mathrm{P}\mathrm{H}\mathrm{Y}\mathrm{T}\mathrm{O}$$

(4)

$${\mathrm{P}\mathrm{a}\mathrm{r}}_{\mathrm{l}\mathrm{i}\mathrm{m}}=\frac{\mathrm{P}\mathrm{A}\mathrm{R}}{\mathrm{P}\mathrm{A}\mathrm{R}+{\mathrm{k}\mathrm{s}}_{\mathrm{P}\mathrm{A}\mathrm{R}}}$$

(5)

$${\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{i}\mathrm{m}}={\mathsf{\theta }}^{\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}-{\mathrm{T}}_{\mathrm{o}\mathrm{b}\mathrm{s}}}$$

(6)

$${\mathrm{D}\mathrm{I}\mathrm{N}}_{\mathrm{l}\mathrm{i}\mathrm{m}}=\frac{\mathrm{D}\mathrm{I}\mathrm{N}}{\mathrm{D}\mathrm{I}\mathrm{N}+{\mathrm{k}\mathrm{s}}_{\mathrm{D}\mathrm{I}\mathrm{N}}}$$

(7)

$${\mathrm{P}}_{\mathrm{l}\mathrm{i}\mathrm{m}}=\frac{\mathrm{P}}{\mathrm{P}+{\mathrm{k}\mathrm{s}}_{\mathrm{P}}}$$

(8)

$${\mathrm{S}\mathrm{i}}_{\mathrm{l}\mathrm{i}\mathrm{m}}=\frac{\mathrm{S}\mathrm{i}}{\mathrm{S}\mathrm{i}+{\mathrm{k}\mathrm{s}}_{\mathrm{S}\mathrm{i}}}$$

(9)

where ks_PAR, T_obs, ks_DIN, ks_P, and ks_Si are parameters of the model. The equation for PAR (Equation (10)) is defined based on the Lambert-Beer law [31,32]. where I₀ is the surface irradiance (µEinst m⁻² s⁻¹), z is the depth (m), and k_d is the diffuse attenuation coefficient (m⁻¹).

$$\mathrm{P}\mathrm{A}\mathrm{R}={\mathrm{I}}_0·{\mathrm{e}}^{-\mathrm{k}\mathrm{d}·\mathrm{z}}$$

(10)

$$\theta =1.185-0.00729·Temp$$

(11)

In the model, the surface irradiance is modelled as photosynthetically active radiation, as defined in Soetaert and Herman [19]:

$${\mathrm{I}}_0=0.5·\left(540+440·\sin \left(\frac{2·\mathsf{\pi }·\mathrm{t}}{365-1.4}\right)\right)$$

(12)

where t indicates the day of the year. The diffuse attenuation coefficient (k_d) describes the rate at which light diminishes with depth due to absorption and scattering in the water column [31,32,45]. K_d is used as a proxy of the influence of Suspended Particulate Matter (SPM) on the PAR [45]. The K_d values were calculated based on the PAR and depth data recorded in the Marine Information and Data Acquisition System (MIDAS). The depth (z) used in Equation (10) was at three metres as the LifeWatch data are measured at a three-metre depth [23].

Grazing represents the grazing of zooplankton (ZOO) on phytoplankton (PHYTO), i.e., how much biomass is eaten, and it is expressed as follows:

$$\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}={\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}}·\frac{\mathrm{P}\mathrm{H}\mathrm{Y}\mathrm{T}\mathrm{O}}{\mathrm{P}\mathrm{H}\mathrm{Y}\mathrm{T}\mathrm{O}+{\mathrm{k}\mathrm{s}}_{\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}}}·\mathrm{Z}\mathrm{O}\mathrm{O}$$

(13)

where max_grazing is the maximum phytoplankton biomass that can be eaten by zooplankton, and ks_grazing is the half-saturation constant for grazing.

Zooplankton faeces production (Faeces) is expressed as a constant fraction (pFaeces) of the total grazing.

$$\mathrm{F}\mathrm{a}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{s}=\mathrm{p}\mathrm{F}\mathrm{a}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{s}·\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{z}\mathrm{i}\mathrm{n}\mathrm{g}$$

(14)

Zooplankton excretion (Excretion) is estimated as the excretion rate (excretionRate) multiplied by the zooplankton biomass.

$$\mathrm{E}\mathrm{x}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}=\mathrm{e}\mathrm{x}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e}·\mathrm{Z}\mathrm{O}\mathrm{O}$$

(15)

Zooplankton mortality is determined by multiplying the mortality rate (mortalityRate) with the zooplankton biomass squared.

$$\mathrm{M}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}=\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\mathrm{R}\mathrm{a}\mathrm{t}\mathrm{e}·\mathrm{Z}\mathrm{O}{\mathrm{O}}^2$$

(16)

The phytoplankton biomass is converted to a chlorophyll-a concentration by multiplying the phytoplankton biomass with the chlorophyll-a-to-nitrogen ratio.

$$\mathrm{C}\mathrm{h}\mathrm{l}\mathrm{o}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{h}\mathrm{y}\mathrm{l}\mathrm{l}=\mathrm{c}\mathrm{h}\mathrm{l}\mathrm{N}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}·\mathrm{P}\mathrm{H}\mathrm{Y}\mathrm{T}\mathrm{O}$$

(17)

The maxUptake, excretionRate, maxGrazing, ks_Grazing, pFaeces, mortalityRate, and chlNratio are parameters of the model of which the seasonal ranges are provided for each region (

Table 1. The seasonal minimum and maximum values of the thirteen parameters used as the input for the nutrient–phytoplankton–zooplankton model for each region.

Parameter	Unit	Period	Nearshore Region	Midshore Region	Offshore Region
maxUptake	day−1	Spring	0.38–0.66	0.38–0.61	0.50–1.12
		Autumn	0.40–0.78	0.38–0.80	0.38–0.90
excretionRate	day−1	Spring	0.16–0.18	0.12–0.16	0.11–0.15
		Autumn	0.11–0.17	0.11–0.15	0.11–0.14
maxGrazing	day−1	Spring	0.87–0.96	0.85–0.92	0.88–0.97
		Autumn	0.88–0.96	0.85–0.93	0.89–0.97
ksGrazing	mmol N m−3	Spring	2.15–3.27	1.54–2.15	1.48–2.22
		Autumn	1.31–2.27	1.19–1.59	1.25–1.94
pFaeces	day−1	Spring	0.29–0.41	0.27–0.40	0.27–0.40
		Autumn	0.25–0.40	0.24–0.32	0.25–0.37
mortalityRate	(mmol N m−3)−1 day−1	Spring	0.28–0.39	0.28–0.41	0.29–0.41
		Autumn	0.32–0.44	0.35–0.42	0.33–0.45
ChlNratio	mg Chla (mmol N)−1	Spring	7.00–7.86	6.78–7.60	6.65–7.47
		Autumn	6.62–7.55	5.33–6.84	4.33–6.61
ksPAR	µEinst m−2 s−1	Spring	126–227	133–224	103–210
		Autumn	121–205	126–200	115–210
Tobs	°C	Spring	9.86–13.84	10.11–13.29	9.54–13.62
		Autumn	10.41–12.83	9.66–13.86	10.08–13.62
ksDIN	mmol N m−3	Spring	1.33–4.21	1.62–3.94	1.92–3.70
		Autumn	1.17–3.64	2.22–4.11	2.07–4.29
ksP	mmol P m−3	Spring	0.30–0.43	0.28–0.44	0.30–0.44
		Autumn	0.33–0.40	0.29–0.46	0.30–0.44
ksSi	mmol Si m−3	Spring	0.41–0.66	0.43–0.67	0.34–0.67
		Autumn	0.35–0.63	0.40– 0.67	0.39–0.65
Kd *	m−1	Spring	0.73–0.90	0.44–0.60	0.28–0.38
		Autumn	0.77–0.92	0.45–0.60	0.28–0.38

* Calculated at a three-metre depth.

).

2.4. Selection of Model Parameters and Validation

The NPZ model requires parameterisation of thirteen parameters (Table 1 and

Table A4. The body mass per individual (mg C ind⁻¹) of the most common species of taxa found in the Belgian part of the North Sea.

Taxon	Body Mass (mg C ind−1)	Species	Reference
Calanoida	0.0006	Acartia clausi, Temora longicornis, Paracalanus parvus, Centropages hamatus, Pseudocalanus elongatus, Centropages typicus and Calanus helgolandicus	[116]
Noctiluca	0.0003	Noctiluca scintillans	[117]
Harpacticoida	0.001	Euterpina acutifrons	[118]
Appendicularia	0.002 to 0.006	Oikopleura dioica	[119]

). For the model parameterisation, we followed a two-step approach to find the most optimal model calibration for each region of the BPNS to mimic the biogeochemical processes in each respective region (Figure 3). In a first step, 5000 unique sets of parameters, which were selected using a Monte Carlo approach, were run for each region of interest. The initial minimum and maximum parameter values used to define these unique sets of parameters were based on values reported in the literature (

Table A1. Parameterisations of the nutrient–phytoplankton–zooplankton (NPZ) model, i.e., the minimum and maximum values, units, and references. These intervals are used as a reference to create the different sets to calibrate the NPZ model. * Minimum and maximum values for kd were calculated for each region of interest based on CTD data recorded in MIDAS.

Parameter	Minimum Possible Value	Maximum Possible Value	References
maxUptake	0.25 (day−1)	1.5 (day−1)	[19,108,109,110,111]
excretionRate	0.1 (day−1)	0.2 (day−1)	[19,109,111]
maxGrazing	0.8 (day−1)	1 (day−1)	[19,108,109,112,113]
ksGrazing	1 (mmol N m−3)	4 (mmol N m−3)	[19,113,114]
pFaeces	0.2 (day−1)	0.5 (day−1)	[19]
mortalityRate	0.25 ((mmol N m−3)−1 day−1)	0.5 ((mmol N m−3)−1 day−1)	[19,108,109,110,111,113,115]
ChlNratio	1 (mg Chla/mmol N)	8 (mg Chla/mmol N)	[50]
ksPAR	30 (µEinst m−2 s−1)	250 (µEinst m−2 s−1)	[19]
Tobs	7 °C	15 °C	[8,113,114]
ksDIN	0.25 (mmol N m−3)	5 (mmol N m−3)	[8,19,108,109,110,111,112,113,114,115]
ksP	0.2 (mmol P m−3)	0.5 (mmol P m−3)	[8,113]
ksSi	0.2 (mmol Si m−3)	0.8 (mmol Si m−3)	[49]
Kd *	0.6 (m−1)	1 (m−1)	Nearshore station
	0.27	0.67	Midshore station
	0.21	0.44	Offshore stations

). In a second step, a new set of 5000 unique parameterisations, which were selected using a Monte Carlo approach, were run for each region of interest and for two seasons, i.e., spring and autumn. In this second step, the minimum and maximum parameter values were based on the 10% best parameterisation for either spring or autumn conditions from the first step. To rank these 10% best models, we calculated the Root Mean Square Error (RMSE) by comparing the NPZ model predictions of phyto- and zooplankton density with the observed pigment chlorophyll-a [25], i.e., a proxy for phytoplankton biomass and zooplankton densities [26], respectively. We calculated the total RMSE for each unique parameterisation as the cumulated error for phytoplankton and zooplankton. For each unique model parameterisation that we tested, we cumulated the error over the different time steps. We considered the best parameterisations as those with the lowest 10% RMSE. In this second step, we retained a set of parameters (Table 1) to describe the spring conditions, i.e., the situation after winter solstice and before summer solstice, and a set of parameters to describe the autumn conditions, i.e., the situation after summer solstice and before winter solstice. The two parameter sets were used to describe the changes in the environmental conditions in each season.

To assess the model fit, we compared the NPZ phytoplankton biomass predictions (expressed in mmol N m⁻³) with the observed phytoplankton biomass data (expressed in mg Chla m⁻³). To enable the comparison, we converted the unit of the NPZ phytoplankton biomass predictions and the unit of the LifeWatch phytoplankton biomass data to mg Chla m⁻³. To do so, we used the Chl:N ratio parameter of the corresponding simulation (Table 1). To compare the NPZ zooplankton production predictions (expressed in mmol N m⁻³) with the LifeWatch zooplankton observations (expressed in ind m⁻³), we converted the latter to the same unit. To do so, we calculated the taxon-specific body mass per individual (mg C ind⁻¹), converted this mass to a molar mass (mmol C m⁻³), and used a taxon-specific C:N ratio to finally convert the molar mass to mmol N m⁻³. Details about this conversion are available in Appendix C (Table A4 and

Table A5. The C:N ratio for each of the most common taxa present in the Belgian part of the North Sea.

Taxon	C:N Ratio	Species	Reference
Calanoida	5.5–7	Acartia spp., Temora sp., Centropages, Oithona sp., Pseudo/Paracalanus spp.	[120]
Noctiluca	2.3–4.4	Noctiluca scintillans	[121]
Harpacticoida	4.26–4.74 7.7–8.1	Euterpina acutifrons	[122] [123]
Appendicularia	4.08	Oikopleura dioica	[119]

). In addition, the model predictions were visually inspected by plotting the observed values against the predicted values (Figure A5 and Figure A6). The visual inspection was used to evaluate the seasonality aspect of the model fit and to what extent that the models that obtained a good RMSE also corresponded to the ecological interactions and dynamics in one growth season. This visual inspection was performed to unravel the over- or underestimation of the model.

2.5. Relative Contributions

The relative importance of the SST, PAR, DIN, PO₄, SiO₃, and zooplankton grazing for phytoplankton biomass dynamics were calculated as was performed in Everaert et al. [12]. To do so, we made use of the forcing functions that were integrated in the model (cfr. 2.3). For each determinant, the absolute limitation was calculated as one minus the limitation factor. Then, the relative contribution of each determinant was calculated as the absolute limitation divided by the sum of all absolute limitations. Afterwards, the monthly relative contribution of each determinant was calculated based on the average relative contribution of the 5% best simulations, i.e., the simulations with the lowest RMSE, for each region of interest. The 5% best simulations were selected for calculating the monthly relative contribution as they offered a good compromise between computational cost and performance. The normality and homogeneity of the relative contribution data were tested, using the packages ‘stats’ [30] and ‘lawstat’ [46], by means of the Shapiro–Wilk test (p < 0.05) and the Levene’s test (p < 0.05), respectively. Potential differences in the determinants between regions were examined using the Kruskal–Wallis test (‘stats’ package; [30]) and the Dunn test (‘dunn.test’ package; [47]) in R [30].

3. Results

3.1. Model Fit

The GAMs were used to create a daily time series of the nutrients and the SST, and were based on the monthly and seasonal observed data (Figure A4). The SST had the highest R² values, and the parameter with the lowest R² varied between regions (Table A6).

To assess the model fit, the model predictions of the phyto- and zooplankton were compared to the field observations, i.e., the RMSE was calculated for each unique parameterisation. The nearshore region had a total RMSE of 1.24–1.38. The mid- and offshore regions had a total RMSE of 0.39–0.46 and 0.32–0.42, respectively. We found that the RMSE for the phytoplankton in the nearshore region was 1.09–1.31, midshore was 0.30–0.40, and offshore was 0.26–0.38. The model tended to overestimate the phytoplankton density at low observed densities (Figure A5). At high observed phytoplankton densities, the model performs better (Figure A5). For the zooplankton, the RMSE was 0.07–0.23 in the nearshore, 0.06–0.13 in the midshore, and 0.04–0.09 in the offshore region. A similar pattern was observed for the zooplankton densities (Figure A6). In the offshore region, it was observed that the underestimation at high zooplankton densities was more pronounced (Figure A6).

3.2. Phytoplankton Time Trends

We found clear spatial differences in the phytoplankton biomass dynamics in the BPNS. The closer to the coastline, the higher the amplitude of the spring phytoplankton bloom and the more distinct the seasonal pattern (Figure 4). In the nearshore region, the phytoplankton biomass ranged from 3.5 to 23.1 mg Chla m⁻³ and followed a clear seasonal trend, with the highest chlorophyll-a concentration being observed in spring and the lowest chlorophyll-a concentration found in the winter (Figure 4). In the midshore region, the maximum phytoplankton biomass was estimated to be 8.8 mg Chla m⁻³, and the minimum phytoplankton biomass was 1.8 mg Chla m⁻³ (Figure 4). In the offshore region, the spring blooming periods were still noticeable but less prominent (1.4–5.3 mg Chla m⁻³) in terms of absolute phytoplankton biomass when compared to the nearshore region. In each of the selected regions, the autumn blooming periods were noticeable, but they became relatively more pronounced with an increasing distance from the coastline (Figure 4). Overall, the spring blooms were observed for each region and were followed by a smaller peak at the end of summer (Figure 4).

We observed a response in the zooplankton population to the spring phytoplankton bloom (Figure 5 and Figure A10). As theoretically expected from a classic predator–prey pattern, there is a time lag between the peak in phytoplankton density and zooplankton density (Figure 5 and Figure A10). Note that higher, i.e., almost four times higher, zooplankton densities in the spring blooms were observed in the nearshore region as compared to the offshore region. In the nearshore region, the zooplankton density ranged from 0.002 to 0.41 mmol N m⁻³ (Figure A7); in the mid- and offshore regions, these ranges were 0.003–0.12 mmol N m⁻³ and 0.002–0.10 mmol N m⁻³, respectively (Figure A7).

3.3. Relative Contributions

We found that solar irradiance and zooplankton grazing are the most important determinants of phytoplankton biomass throughout the year (Figure 6). Together they contribute 38% to 77% to the phytoplankton biomass dynamics in the BPNS. The contribution of zooplankton grazing is the highest after the spring blooms (31%) and during the autumn months (35%). PAR plays a more important role in the phytoplankton biomass during, autumn (43%). We found that the nutrients and SST play a relatively less important role in phytoplankton biomass dynamics. The total contribution of nutrients is, at maximum, 44%. During spring bloom, PO₄ plays an important role, while SiO₃ and DIN take over during early summer (Figure 6 and Appendix G). PO₄ is the most limiting nutrient (12–29%), followed by DIN (3.6–26%) and SiO₃ (1.2–15%), respectively. The SST only plays a limiting role during winter (max. 17%, Figure 6).

The relative contributions gradually changed along the nearshore–offshore transect (Figure 6 and Figure 7). For example, the PAR’s relative contribution is the highest in the nearshore region (30–43%) and decreases towards the open sea (14–34%). Zooplankton grazing is less important in the offshore (22–31%) than in the nearshore region (21–35%) or the midshore region (23–34%). There is also a nearshore–offshore gradient in terms of the relative contribution of nutrients (Figure 7). For each of the individual nutrients, i.e., DIN, PO₄, and SiO₃, we found them to be more limiting in the offshore region (8.5–26%; 15–29%; 3.9–15%, respectively) than in the nearshore region (3.7–11%; 12–27%; 1.2–13%, respectively). The SST was equally limiting in all three regions (Figure 7).

4. Discussion

Using an NPZ model, we reproduced phytoplankton and zooplankton dynamics in the BPNS (Figure 4 and Figure 5), and we quantified the relative contribution of the key determinants of the phytoplankton biomass, i.e., the nutrients, solar irradiance, SST, and zooplankton grazing (Figure 6 and Figure 7). This was conducted for three regions, i.e., the near-, mid-, and offshore, to examine the spatiotemporal variation. Only a few studies, e.g., Everaert et al. [12], Llope et al. [13], and McQuatters-Gollop et al. [14], have quantified the temporal relative contribution of phytoplankton biomass’ key determinants. We found clear regional differences and seasonal patterns in the relative contribution of the key determinants of phytoplankton biomass dynamics.

4.1. Generalised Additive Modelling Performance

The GAMs were created based on the available monthly data and their model fit was assessed following the method of Zuur et al. [29]. Using monthly data meant that the performance of the models was restricted by the availability of data. Although the R² was not high, the smoothers created based on the GAMs captured the expected seasonality changes that were observed in the nutrient data. It was possible to capture the seasonality changes in the nutrient input data but within a reduced variance of the possible range for the input data. This could have affected our phytoplankton biomass dynamics if the nutrients were underestimated in the GAMs. Phytoplankton blooms may have been simulated smaller when the nutrient concentrations were estimated lower than they were in reality. This may have led to a flattening of the phytoplankton biomass dynamics, i.e., the extreme peak values may have been missed.

4.2. Comparison Phyto- and Zooplankton Modelling Results

The in situ chlorophyll-a concentrations of this study (nearshore: 1.1–57.5 mg Chla m⁻³, midshore: 0.6–13.1 mg Chla m⁻³, and offshore: 0.2–10.3 mg Chla m⁻³) were close to other field studies in the BPNS (Figure 8), such as Desmit et al. [9] (nearshore: 1–20 mg Chla m⁻³, and offshore: 0.5–15 mg Chla m⁻³) and Muylaert et al. [48] (nearshore: 0.2–60 mg Chla m⁻³, midshore: 0.2–20 mg Chla m⁻³, and offshore: 0.2–15 mg Chla m⁻³). The NPZ model predictions (nearshore: 3.5–23.1 mg Chla m⁻³, midshore: 1.8–8.8 mg Chla m⁻³, and offshore: 1.4–5.3 mg Chla m⁻³) based on the chlorophyll-a concentrations observed in this study are close to the modelling results of Arndt et al. [8] (0.5–40 mg Chla m⁻³) and Lancelot et al. [49] (0.3–23 mg Chla m⁻³; Figure 8). However, our model predictions did not show extreme, i.e., high and/or low, phytoplankton densities, thus resulting in the phytoplankton biomass predictions being close to the modelled-based quantifications found in the literature (Figure 8). In neighbouring regions, similar chlorophyll-a values were observed by, e.g., Alvarez-Fernandez et al. [50], Lancelot et al. [49], EEA [51], Colella et al. [52], and Lundsør et al. [53] (Figure 8). Although the chlorophyll-a concentrations are in line with previous studies, we noticed that there is some variability between previous studies and our observations. A potential reason for this is that the data used in this study are more recent. Indeed, with the inclusion of recent observations, Desmit et al. [9] found a decrease in the annual mean chlorophyll concentration for offshore regions over a time span of 40 years. This is supported by Xu et al. [54], thereby demonstrating a decreasing trend in the chlorophyll-a in the offshore region of the central North Sea.

In the BPNS, our model indicated a clear seasonal pattern with low phytoplankton biomass in winter (min. 1.4 mg Chla m⁻³) and increasing phytoplankton biomass during spring. This spring bloom (max. 23.1 mg Chla m⁻³), typically consisting of diatoms and Phaeocystis spp. [48], occurs in March and April and is followed by a smaller bloom in autumn (max. 7.7 mg Chla m⁻³). We found a decrease in phytoplankton biomass overall, and a decrease of the amplitude of the spring bloom (23.1 mg Chla m⁻³ nearshore to 5.3 mg Chla m⁻³ offshore) with increasing distances from the coast (Figure 4). These regional differences were also observed by Desmit et al. [9] (20 mg Chla m⁻³ nearshore to 12 mg Chla m⁻³ offshore) and Muylaert et al. [48] (60 mg Chla m⁻³ nearshore to 15 mg Chla m⁻³ offshore) in the BPNS. Muylaert et al. [48] and Desmit et al. [9] have also observed that the seasonal pattern, i.e., a spring bloom followed by a smaller autumn bloom, was more distinct when closer to the coast (Figure 4). Jiang et al. [55] found a high interannual variability in the peak biomass, which is also observed in our field observations but is less expressed in our modelling results. Nevertheless, in nearshore regions we see that the autumn bloom is more modest, i.e., three to four times smaller in amplitude, than the spring bloom.

The classic bimodal bloom pattern that we (Figure 5) and others (e.g., Muylaert et al. [48] and Lancelot et al. [49]) have observed in the BPNS was not found by Nohe et al. [56] between 2003 and 2010. They found that the spring bloom was more intense and extended, and that there was no autumn bloom. Nohe et al. [56] suggested that increased SST and water transparency, as well as changes in nutrient concentrations and ratios, are potential reasons for the lack of an autumn bloom. This lack of an autumn bloom is in strong contrast with our findings with more recent data, i.e., 2014–2017. Similar to our study, Speeckaert et al. [57] found two blooms, i.e., a spring bloom followed by a smaller autumn bloom in 2016, whereas Nohe et al. [56] observed a mean diatom cell density of 3.9 × 10⁵ cell L⁻¹ in the autumns in the period of 2003–2010. In addition, Speeckaert et al. [57] found a peak cell density of 2.0 × 10⁶ cell L⁻¹ with Guinardia spp., which is the dominant diatom in the autumn bloom in 2016. This autumn bloom is also observed in the data from the LifeWatch Flowcam [58,59] for 2017 and later (Figure A12, Figure A13 and Figure A14). This suggests that the BPNS may have shifted back to a two-bloom pattern.

The seasonal dynamics of zooplankton correspond largely with the data found by van Ginderdeuren et al. [27], Mortelmans et al. [24], and Deschutter et al. [60] for copepods, which largely dominate the zooplankton community in the BPNS [27,61]. The delayed increase in zooplankton density corresponding with the spring bloom, illustrating the zooplankton grazing on phytoplankton in our model (Figure 5 and Figure A10 [49]), depicts the population dynamics of a classic predator–prey relationship [62]. Zooplankton grazing is an important determinant of phytoplankton biomass dynamics, but the predator–prey relationship is complex and changes under different abiotic conditions [63].

The LifeWatch observations and our model results suggest higher zooplankton density in nearshore regions (Figure 4), and this agrees well with the findings of Mortelmans et al. [64] and of Deschutter et al. [60] for the BPNS. This was also observed in other areas in the world [65,66]. Van Ginderdeuren et al. [27] observed that the highest copepod densities occurred in their ‘midshore region’ of the BPNS, but that region overlaps with our nearshore region (Figure 1).

4.3. Relative Contribution of the Key Determinants

This study provides a better understanding on the spatiotemporal variation in the relative contribution of the key determinants of marine phytoplankton biomass dynamics in the BPNS through using an NPZ model that is driven by in situ observation. The relative contributions of key determinants to changes in the phytoplankton biomass dynamics are highly seasonal and spatially related in the BPNS.

The relative contribution of the determinants changes seasonally (Figure 6). During autumn and winter, when there is the least amount of light in the BPNS (±8 h per day), the PAR is the major determinant (max. 43%). Obviously, sufficient solar irradiance is needed for photosynthesis [67]. The SST becomes more of a limiting factor for phytoplankton growth during mid-winter and early spring (max. 17%). This was also found by Everaert et al. [12]. They found that the combined relative contribution of the SST and PAR varied from 20% during summer to 50% during winter, which is similar to our results (summer: 14% and winter: 51%). In shallow coastal areas like the BPNS, the SST is a key factor affecting phytoplankton bloom dynamics [68] as the low temperature decreases the growth rate of marine phytoplankton [69]. In early spring, when solar irradiance and temperature increases, phosphorus becomes depleted due to the spring bloom [8,49,70] and becomes the major determining factor (29%). Nitrogen becomes increasingly limiting in early summer (26%), together with silica in the nearshore (13%) and midshore (10%) regions. The increased phytoplankton density during spring bloom results in a higher zooplankton density after the spring bloom in early May [64]. Together with zooplankton density, the grazing pressure on phytoplankton increases. During summer and autumn, zooplankton grazing is an important determinant with a high relative contribution on the phytoplankton biomass (31% and 35%, respectively); this is because zooplankton biomass remains high [64], whereas phytoplankton biomass decreases, thus changing the phytoplankton–zooplankton ratio. This corresponds with the findings of Gowen et al. [71] for the Irish Sea, who found that the percentage of phytoplankton biomass grazing is highest after spring bloom in May. From the winter onwards, zooplankton grazing’s relative contribution decreases as their abundance has decreased together with their food supply, i.e., phytoplankton, which is mainly limited again by the low solar irradiance completing the seasonal cycle of the determinants (Figure 6).

Besides their temporal variability (cfr. previous paragraph), the key determinants affecting phytoplankton biomass also vary spatially. The influence of nutrients on the phytoplankton biomass dynamics increases with increasing distance from the coastline (Figure 7). The main reason for this is related to the fact that nearshore regions tend to be nutrient-rich due to riverine discharges. As the river runoff creates a nearshore–offshore gradient in nutrients with lower nutrient availabilities in offshore regions (Figure A4, [8,70], all nutrients have a higher relative contribution in limiting phytoplankton biomass (Figure 6 and Figure 7), i.e., nitrogen (8.5–26%), phosphorus (15–29%), and silica (3.9–15%). The difference in relative contributions between nearshore and offshore regions is most obvious for nitrogen, i.e., DIN is being more limiting in the offshore region (max. 26%) than in the nearshore (max. 11%) region. Phosphorus, on the other hand, limits phytoplankton biomass in the offshore region for a longer period. This could be related to the water depth as phosphorus is largely regenerated from the sediment; additionally, for the offshore region, there is a much higher water mass-to-sediment ratio [70]. Silica becomes a limiting factor earlier in the offshore region than in the near- and midshore regions, as the low SiO₃ reserve in the offshore region is depleted quickly during the spring bloom [48]. The gross part of riverine nutrient resources (the Scheldt, Rhine, Meuse, and Seine [18]) is depleted before it reaches the offshore region [8]. Our results largely agree with the findings of Burson et al. [72], i.e., nearshore regions in the North Sea are P limited and, the mid- and offshore regions of the BPNS are N and P co-limited. The impact of zooplankton grazing on phytoplankton biomass is smaller in offshore regions, which could be due to the lower zooplankton density (Figure A7), and thus the lower grazing pressure. In the nearshore region, PAR is a major determinant year-round (30–43%); whereas, in the mid- and offshore regions (14–34%), this is more in balance with the other determinants (Figure 7). The high turbidity in the nearshore region [18,73] restricts the phytoplankton’s access to sunlight and causes PAR to be a major determining factor throughout the year. The more offshore the region, the less influential PAR is in terms of phytoplankton biomass dynamics.

Overall, the relative contributions of the key determinants of phytoplankton biomass dynamics showed a clear spatiotemporal variation. The spatiotemporal variation is partially driven by meteorological conditions, which affects currents, temperature, the levels of SPM, light availability, river discharge, and nutrient dynamics in the BPNS [74]. Nutrients from various rivers such as the Seine and Somme are carried by the influx of Atlantic water to the BPNS [18]. Local rivers, primarily the Scheldt, Rhine, and Meuse—which flow through industrialised and densely populated areas, thus providing them with high nutrient loads—further contribute to the nutrient status of the BPNS [75]. Consequently, there is a noticeable southeast–northwest gradient in nutrient concentrations, which was also observed in the relative contributions of the nutrients (Figure 6 and Figure 7). Similar to the nutrient gradient, a salinity gradient from east to west is created due to the Scheldt [16]. The BPNS also exhibits a gradient in PAR, primarily due to spatial variations in the SPM and in the dissolved organic matter. In the vicinity of the coastline, shallow waters tend to be more turbid and have higher SPM concentrations [75,76]. These factors mainly drive the spatial gradient in phytoplankton biomass dynamics (Figure 4) in the BPNS.

Figure 8. An overview of the phytoplankton biomass concentrations in the Belgian part of the North Sea (BPNS, yellow) [8,9,48,49,57,76], Greater North Sea (green) [9,49,50,53], Mediterranean Sea (blue) [52], and the North East Atlantic and Baltic Sea (purple) [51]. The phytoplankton biomass concentrations are subdivided by region, the near-, mid-, and offshore regions in the BPNS, as well as the coastal and open sea regions in the Greater North Sea, Mediterranean Sea, and North East Atlantic and Baltic Sea. Citations on the y-axis refer to the corresponding studies. For the present study, both model-based simulated phytoplankton biomass (This study (simulated)), as well as in situ observations (This study (observed)) have been integrated.

Figure 8. An overview of the phytoplankton biomass concentrations in the Belgian part of the North Sea (BPNS, yellow) [8,9,48,49,57,76], Greater North Sea (green) [9,49,50,53], Mediterranean Sea (blue) [52], and the North East Atlantic and Baltic Sea (purple) [51]. The phytoplankton biomass concentrations are subdivided by region, the near-, mid-, and offshore regions in the BPNS, as well as the coastal and open sea regions in the Greater North Sea, Mediterranean Sea, and North East Atlantic and Baltic Sea. Citations on the y-axis refer to the corresponding studies. For the present study, both model-based simulated phytoplankton biomass (This study (simulated)), as well as in situ observations (This study (observed)) have been integrated.

Whereas the temporal variation in phytoplankton biomass, i.e., the blooms (Figure 4), are mainly determined by temperature and PAR (Figure 6) [63]. Increased temperature resulting from more solar irradiation promotes phytoplankton growth. Additionally, the increased settling of SPM in spring increases the light availability. The temporal changes in PAR are influenced by the resuspension and movement of particulate matter, which are further modified by winds and tides. For example, the stronger winds during winter lead to a higher attenuation of PAR. This settling process affects the timing of the onset of spring phytoplankton blooms [75,76].

The temporal variation in phytoplankton biomass is also, to a lesser extent, determined by the mixing of the water column, nutrient conditions, grazing, and phytoplankton community composition [63]. For example, in most areas of the region, the BPNS exhibits a constant mixing of the water column throughout the year due to a combination of robust tidal currents and shallow depth [15,77].

Additionally, there are also factors that could affect the phytoplankton biomass dynamics, such as anthropogenic activities [78,79,80,81,82] and biological agents [83,84,85]. Biological agents, such bacteria and viruses can affect the bloom development and termination [83,85]. For example, the rapid proliferation of viruses can have a strong influence on the populations of phytoplankton hosts, thus playing a crucial role in marine biogeochemistry and ecology [83].

Quantifying the relative contribution and identifying the spatiotemporal variation offer a better understanding of how the key determinants’ limitations to phytoplankton biomass will change under changing conditions, e.g., those related to climate change. However, every method has its advantages and limitations, and below we consider those of our modelling approach before addressing the potential implications of this study.

4.4. Modelling with Field Data: Advantages and Limitations

Often multiple driver-related research has been performed in laboratories. A big advantage of this is that model species can be kept in optimal conditions, e.g., temperature, light, nutrients, etc., which aids in isolating the effects of the stressor in question. However, these optimal conditions are rarely experienced by organisms in their natural environment [86], thus hampering the conversion of the laboratory-based conclusions towards field conditions. Therefore, we have selected a modelling approach using field data, which has the advantage that the natural background variation of bottom-up drivers of marine ecosystems are implicitly included in the results [87]. Therefore, by using field data to quantify the relative contribution, we can assess the impact of the determinants under their natural and continuously changing conditions.

Even though the phytoplankton biomass dynamics modelled in this study agree well with the field observations from LifeWatch (Figure 4), we acknowledge that our modelling approach contains limitations: (i) The periods with low phytoplankton growth in the offshore region correspond less with the field observations. This is likely attributable to the seasonal sampling strategy in that region, i.e., no monthly measurements [23]. As such, we had less data available for the offshore region to calibrate our model. (ii) In this study, only one station, i.e., station 130, for the nearshore region, and one station, i.e., 330, for the midshore region were selected to represent their region. Data exploration on the publicly accessible LifeWatch server shows that the patterns seen in the data from station 130, chosen to represent the nearshore area, were also evident in the data from other nearshore stations, such as station 120 in the west and station 700 in the east. This conclusion also applies to the midshore region, currently represented by station 330, as similar dynamics were found in station ZG02 and station W07bis. (iii) For the NPZ model, different functional groups of phytoplankton and zooplankton were grouped, and—by doing so—we may have missed species-specific limitations. It is known that reactions to changes in determinants, such as with climatic changes, are quite species-specific [88,89]. (iv) Additionally, only the most abundant and dominant zooplankton species, i.e., copepods and Appendicularia, were aggregated. However, they account for 76% of the total zooplankton density [27], and, as such, represent the main biomass dynamics of zooplankton. (v) Chl:N ratios are known to have a seasonal trend [90]. This seasonality was accounted for by selecting two Chl:N ratios, i.e., one ratio to reflect spring conditions and one ratio for autumn conditions. These ranges were selected based on minimising the RMSE for each season. As such, after calibration, model simulations were performed on fixed Chl:N ratios for each season. Using the fixed Chl:N ratios could have caused a flattening of the curve, in both high (spring bloom peak) and low (in winter) extreme values in phytoplankton densities. Even though Chl:N ratios were selected for each region individually, the pattern of over- (during winter) and underestimating (during spring) was found at each region (Figure A5). The over- (during winter) and underestimating (during spring) was most pronounced in the nearshore region (Figure A5). Additionally, the overestimation of phytoplankton density at low observed densities, i.e., in winter (Figure A5), may have resulted from the selected Chl:N ratio for the autumn conditions (Table 1) (which was in the upper part of the range found in the literature (Table A1)). (vi) The NPZ model does not include a horizontal and vertical mixing of the nutrients in the water column. Nevertheless, these were accounted for: firstly, the BPNS is permanently mixed [16]; and secondly, the model was calibrated and fitted using in situ observations, which contain vertical and horizontal mixing characteristics [23]. However, when the model would be applied to another area where stratification does occur, vertical and horizontal mixing should be accounted for. (vii) Using a fixed sinusoidal surface irradiance in our small study area, i.e., 67 km to 65 km, did not impact our calculations as the difference on such a small scale is negligible. When applying the model to large areas, the spatial properties of this factor should be considered.

Although other coastal and regional sea models exist for the North Sea, e.g., ERSEM, MIRO, and ECOHAM [49,91,92], we used the NPZ model to make the best use of the LifeWatch datasets. The complexity of the aforementioned models was not suited for our goal, that is, to develop a model to assess which factors drive marine phytoplankton biomass dynamics in the BPNS and how their relationship to primary production varies on a spatiotemporal scale.

Using this NPZ model allowed us to integrate and simulate with forcing functions. In turn, integrating the forcing functions gave us the possibility to determine the relative contributions, and it permits us to add forcing functions in the model. In addition, this NPZ model can be applied on a small scale, which is needed as the difference in the parameterisation of the regions (Table 1) illustrate. Models covering larger areas produce more general patterns. These general patterns have less detailed information on local differences that is required to enable policymakers to make local decisions. We provided here the first step towards such a model that adheres to the need to understand the consequences of blue economy investments in the BPNS. An advantage of the NPZ model, being less complex, is that it requires less computing power and time.

In addition, some components, e.g., the vertical mixing and depth levels, which are included in the more complex models such as ERSEM [91], are less applicable to the BPNS. The BPNS is permanently mixed [16], and photosynthesis is restricted to the top layer as available light usually does not penetrate the water column more than ten meters [93]. More complex models can provide a more detailed and comprehensive set of outputs than NPZD models as they include a wider range of physical, chemical, and biological variables. However, this could also make it more complex for the end users, such as industry workers and policymakers, to interpret and analyse the model results. Overall, using an NPZ model was selected as the optimum between computing time, interpretability, and complexity in terms of addressing our research topic.

The novelty of this study lies in five aspects: (1) performing a high-resolution spatiotemporal analysis that contributes to a better understanding of phytoplankton biomass dynamics in the BPNS (cfr. 4.2), (2) quantifying the relative contribution of the determinants of phytoplankton biomass dynamics along a near-offshore gradient (cfr. 4.3), and (3) the observation of a bimodal bloom pattern in the BPNS after being absent from 2003 until 2010 (cfr. 4.2). Furthermore, this study has the potential to support dedicated management strategies for a sustainable development of Belgium’s blue economy (see Section 4.5). Finally, as the fifth aspect, the modelling scripts and data are available according to the FAIR principles (see Data Availability Statement Section).

4.5. Future Perspectives and Implications

Our model is well suited to simulating local phyto- and zooplankton dynamics, and it could be used to complete data gaps in phytoplankton and zooplankton observations. Upon the availability of abiotic input variables, these abiotic input variables could be used to drive the model to obtain phytoplankton and zooplankton biomass, thus completing the data gaps. The obtained results should be validated with the available observations from before and after the data gap. In addition, our model and results could be used to evaluate how the key drivers of phytoplankton biomass change under climate change conditions and blue economy activities. As such, this model is suitable for understanding and predicting the consequences of compound events, which are defined as and generated by short time changes (weather related) combined with climate events [94].

During the Anthropocene, human activities have caused ocean warming [96], acidification [95], and eutrophication [96], resulting in an on-average increasing SST, and increased nutrient concentrations compared to the pre-industrial times. To date, it is not clear which long-term changes in abiotic conditions are most important in the context of marine phytoplankton biomass dynamics. It is hypothesised that the relative contribution of the dominant determinants will change due to changing environmental conditions [97]. Modelling and quantifying the relative contribution of the phytoplankton’s determinants provide a more holistic view rather than a one-to-one relation that is obtained from directly measuring phytoplankton biomass. The quantification of relative contribution could give more insight into the underlying mechanisms, e.g., changes in the relative contribution of the key determinants of phytoplankton biomass may indicate a disruption of the phytoplankton community [98,99]. Modelling climate change scenarios should always be performed carefully as (i) the NPZ model does not include different species and could miss changes in the community composition, e.g., due to a poleward shift [98], and (ii) the natural variability of the phytoplankton dynamics may not be fully covered in the four years analysed in this study [100].

In the context of the blue economy, our model contributes to a better understanding and quantification of the dynamics of the local ecosystem and the potential effects of human activities on the marine environment. Here, we present two concrete blue economy cases in which the model can be used to assess the potential consequences of aquaculture development and wind farms on the local ecosystem. In the case of bivalve aquaculture, it is anticipated that nutrient fluxes will be adapted through bottom-up effects [78,79], as well as the limitation of phytoplankton production through filter feeding [80,81]. Once a wind farm is operational, a well-known effect that is generated is the turbidity plume [82]. The decrease in light availability may lead to a decrease in phytoplankton productivity locally. These changes in abiotic conditions can be translated via the NPZ model to phytoplankton biomass dynamics. Overall, our model has the potential to evaluate the effectiveness of various management strategies while being less complex, which could translate to it being user friendly and having more users, e.g., scientists and policymakers.

The model could benefit from (i) using different sources for validation, (ii) expanding the scope of the model, and (iii) by incorporating a greater amount and more accurate data. Firstly, for the validation of the model, high-resolution remote sensing data could be used. This would work well with climate change conditions [101]. Secondly, based on the results of this local model, it would be possible to integrate more parameters, and more complexity could also be built, such as hydrodynamics. This is conducted, for example, in the ERSEM and MIRO models, where the NPZD model is incorporated into a complex model. Lastly, we provided the first step in this direction with our NPZ model, which has the potential to be further developed to a more detailed taxonomical level, e.g., using CHEMTAX as this allocates chlorophyll-a to phytoplankton groups by marker pigments and Chla ratios [102,103,104]. It would also be possible to include, in a further iteration, more detailed data, e.g., data on a more detailed taxonomic level, to have a better understanding of the effects of each of the key determinants on the various subgroups of phytoplankton to predict whether the ecosystem will change [105] (even though phytoplankton biomass dynamics are resilient to a certain extent [106]).

5. Conclusions

In this study, a first step towards a model that adheres to the need to understand the consequences of blue economy investments in the BPNS, is provided. We quantified phytoplankton biomass dynamics and the relative contribution of its determinants using an NPZ model along a near–offshore gradient over a period of four years (2014–2017). By doing so, a better understanding of the spatiotemporal variations in these contributions was provided. After not being observed between 2003 and 2010, we again found a classic bimodal bloom pattern. Further, we found that the relative contributions of phytoplankton determinants alter spatially and temporally. Solar irradiance (up to 43%) and zooplankton grazing (up to 35%) are the most influential determinants for phytoplankton biomass, and this is the case throughout the year. A clear spatial gradient was observed for most of the determinants, e.g., nutrients and zooplankton grazing, which are greater limiting factors offshore, while the opposite is true for solar irradiance. This comprehensive study contributes to a better understanding of phytoplankton biomass dynamics and its determinants along a near-offshore gradient in the BPNS. Additionally, it has the potential to support blue economy strategies and management in the BPNS. In the scope of future climate scenarios and blue economy, it is important to understand the key determinants of phytoplankton biomass dynamics to develop management scenarios that ensure a sustainable marine ecosystem.