Drivers of the spatial phytoplankton gradient in estuarine–coastal systems: generic implications of a case study in a Dutch tidal bay

As the primary energy and carbon source in aquatic food webs, phytoplankton generally display spatial heterogeneity due to complicated biotic and abiotic controls; however our understanding of the causes of this spatial heterogeneity is challenging, as it involves multiple regulatory mechanisms. We applied a combination of field observation, numerical modeling, and remote sensing to display and interpret the spatial gradient of phytoplankton biomass in a Dutch tidal bay (the Eastern Scheldt) on the east coast of the North Sea. The 19 years (1995–2013) of monitoring data reveal a seaward increasing trend in chlorophyll-a (chl a) concentrations during the spring bloom. Using a calibrated and validated three-dimensional hydrodynamic–biogeochemical model, two idealized model scenarios were run: switching off the suspension feeders and halving the open-boundary nutrient and phytoplankton loading. Results reveal that bivalve grazing exerts a dominant control on phytoplankton in the bay and that the tidal import mainly influences algal biomass near the mouth. Satellite data captured a post-bloom snapshot that indicated the temporally variable phytoplankton distribution. Based on a literature review, we found five common spatial phytoplankton patterns in global estuarine–coastal ecosystems for comparison with the Eastern Scheldt case: seaward increasing, seaward decreasing, concave with a chlorophyll maximum, weak spatial gradients, and irregular patterns. We highlight the temporal variability of these spatial patterns and the importance of anthropogenic and environmental influences.

1 Introduction

As the most important energy source in aquatic systems, phytoplankton account for 1 % of the global biomass but create around 50 % of the global primary production (Boyce et al., 2010). Located at the land–ocean interface, estuarine–coastal systems, including estuaries, bays, lagoons, fjords, river deltas, and plumes, are relatively productive and abundant in phytoplankton (Carstensen et al., 2015). As the basis of the pelagic food web, phytoplankton have an immense impact on the biogeochemical cycles, water quality, and ecosystem services (Cloern et al., 2014). A sound understanding of the spatial variability of phytoplankton is critical for effective assessment, exploitation, and protection of estuarine–coastal ecosystems but remains a challenge due to the complicated natural and anthropogenic controls (Grangeré et al., 2010; Srichandan et al., 2015).

The standing stock of phytoplankton is a function of sources and sinks that are subject to both biotic and abiotic influences (Lancelot and Muylaert, 2011; Jiang et al., 2015). Phytoplankton growth is regulated by bottom-up factors such as nutrients, light, and temperature (Underwood and Kromkamp, 1999; Cloern et al., 2014), while natural mortality and grazing pressure from zooplankton, suspension feeders, and other herbivores contribute to the loss of phytoplankton biomass (Kimmerer and Thompson, 2014). Physical transport can act as either a direct source or sink, driving algal cells into or out of a certain region (Martin et al., 2007; Qin and Shen, 2017). The hydrodynamic conditions also affect the phytoplankton biomass indirectly; for example, phytoplankton growth is dependent on the transport of dissolved nutrients (Ahel et al., 1996). Tides and waves affect concentrations of light-shading suspended particulate matter (SPM) and, thus, photosynthesis (Soetaert et al., 1994). The efficiency of benthic filtration feeding on the surface phytoplankton is associated with the stratification of the water column (Hily, 1991; Lucas et al., 2016).

For these reasons, the phytoplankton distribution in estuarine–coastal systems relies on the spatial patterns of physical, chemical, and biological environmental factors of each system (Grangeré et al., 2010). For example, phytoplankton variability in one semi-enclosed water body can be dominated by terrestrial input (river-dominated system), oceanic input (tide-dominated system), top-down effects (grazing-dominated system), other inputs, or a combination of the above factors. More often, it is a delicate balance of multiple factors that determine phytoplankton gradients. Under high river discharge, phytoplankton growth can be promoted by increasing nutrient input, whereas advective loss and high riverine SPM loading may inhibit algal enrichment (Lancelot and Muylaert, 2011; Shen et al., 2019). In tide-dominated systems, tides can resuspend SPM, negatively impacting phytoplankton, and concurrently transport regenerated nutrients into the water column, or they can drive upwelling-induced algal blooms from the coastal ocean into estuaries (Sin et al., 1999; Roegner et al., 2002). Nitrate can support more phytoplankton biomass in microtidal estuaries than in macrotidal estuaries (Monbet, 1992). The relative importance of zooplankton and bivalve grazing on phytoplankton varies spatially (Kromkamp et al., 1995; Herman et al., 1999; Kimmerer and Thompson, 2014). These complexities make it challenging to discern the driving mechanisms of the spatial phytoplankton gradient, and comparative studies of different systems are lacking (Kromkamp and van Engeland, 2010; Cloern et al., 2017).

In situ observation, remote sensing, and numerical modeling are common techniques to reveal spatial patterns and detect their biophysical controls (Banas et al., 2007; Grangeré et al., 2010; Srichandan et al., 2015). Shipboard measurements of chlorophyll a (chl a) provide a precise and dynamic assessment of the phytoplankton variability; however, the temporal (usually monthly) and spatial (usually tens of kilometers) resolutions are limited compared with satellite images and numerical models (Soetaert et al., 2006; Valdes-Weaver et al., 2006; van der Molen and Perissinotto, 2011; Cloern and Jassby, 2012; Kaufman et al., 2018). Remote sensing of chl a reveals the surface distribution with a sufficiently high spatial resolution and coverage, although only under favorable weather conditions (Srichandan et al., 2015). In comparison, ecological models are based on simplified assumptions and numerical formulations and cannot simulate every detail of natural processes. However, a properly calibrated and validated model is capable of representing the system of interest at a fine resolution (Friedrichs et al., 2018) and allows one to test hypotheses related to the mechanisms driving the phytoplankton distribution (Jiang and Xia, 2017, 2018; Irby et al., 2018). As a reliable biophysical model must be based on observational and satellite data (Soetaert et al, 1994; Feng et al., 2015; Jiang et al., 2015), a combination of these approaches is optimal to improve our knowledge of the spatial heterogeneity in estuarine–coastal ecosystems.

In this study, we combined satellite observations, long-term monitoring, and numerical modeling to investigate the potential drivers of the spatial phytoplankton gradients in a well-mixed tidal bay, the Eastern Scheldt (SW Netherlands). In this case study, we identified the main environmental drivers of spatial phytoplankton distribution in the bay and used some sensitivity model tests to quantify the impact of these drivers. Via this mechanistic investigation into the spatial phytoplankton gradient, our case study was then used as a prototype for comparing spatial phytoplankton gradients among global estuaries and coastal bays. Based on a literature review, five main types of spatial heterogeneity in phytoplankton biomass were identified along with examples and dominant controls.

2 The study site

The Eastern Scheldt is located in the Southwest Delta of the Netherlands (Fig. 1). Due to the flood-protection construction projects, known as Delta Works, since the 1980s, the delta region has changed from an interconnected water network to individual water basins isolated by dams and sluices (Ysebaert et al., 2016). The confluence of the Rhine and Meuse rivers flows into the North Sea through a narrow channel (Fig. 1) with a combined discharge of over 2000 m³ s⁻¹ (Ysebaert et al., 2016). The Western Scheldt (Fig. 1) is the only remaining estuary in the delta region covering fresh to saline waters (Ysebaert et al., 2016). With the substantially reduced freshwater input, the Eastern Scheldt has been filled with saline water (salinity 30–33), lost the characteristics of an estuary, and developed into a tidal bay (Nienhuis and Smaal, 1994; Wetsteyn and Kromkamp, 1994). The northernmost end of the northern branch has a salinity level that fluctuates between 28.5 and 30.5, caused by a small freshwater inflow through the Krammer sluice. However, the occasional freshwater flux controlled by the sluice is mostly below 10 m³ s⁻¹ (Ysebaert et al., 2016). This is negligible compared with the tidal exchange, which is ∼ 2 × 10⁴ m³ s⁻¹, estimated from a typical tidal prism of 9 × 10⁸ m³ over a 12 h tidal cycle. The tidal prism is about one-third of the basin volume, ∼ 2.76 × 10⁹ m³ (Nienhuis and Smaal, 1994). As part of the Delta Works, a partially open storm surge barrier was implemented at the mouth of the Eastern Scheldt, which is occasionally closed during severe storms. Since the implementation of the barrier, the tidal basin still experiences a semidiurnal tidal regime, but the average tidal range has been reduced by ∼ 13 % to 2.5–3.4 m from west to east, the tidal flat area has been reduced, and the current velocity has decreased (Nienhuis and Smaal, 1994; Vroon, 1994). In the post-barrier decades, the entire basin has been dominated by the tidal exchange with the North Sea, causing a net import of phytoplankton biomass and seston; the water residence time (RT) of the bay ranges from 0 to 150 d from the western to eastern ends (Jiang et al., 2019a).

Figure 1 Geographic location of the Eastern Scheldt (denoted by “Oosterschelde” in the left panel) and the GETM-FABM model grid, domain, and bathymetry. Green, pink, and red dots in the right panel indicate the distribution of three dominant bivalve species in the Eastern Scheldt: oysters, mussels, and cockles (data source: Wageningen Marine Research).

The phytoplankton composition in the Eastern Scheldt has also changed since the Delta Works: the previously dominant diatoms have decreased, while the small flagellates and weakly silicified diatoms have become more abundant, especially in summer (Bakker et al., 1994). The annual cycle of phytoplankton biomass is characterized by a spring bloom and a much weaker late summer peak (Wetsteyn and Kromkamp, 1994). The Eastern Scheldt has been extensively used for aquaculture of Pacific oysters (Crassostrea gigas) and blue mussels (Mytilus edulis) over the past decades, and their annual production is approximately 3 and 20–40 kt fresh weight, respectively (Smaal et al., 2009; Wijsman et al., 2019). Oysters, mussels, and wild cockles (Cerastoderma edule) are the main benthic suspension feeders in the basin (Fig. 1). Strong pelagic–benthic coupling has been reported for the Eastern Scheldt ecosystem: benthic filtration very likely accounts for the declining annual primary production (Smaal et al., 2013). In addition, abundant benthic suspension feeders make the Eastern Scheldt an important feeding ground and international conservation zone for wading birds (Tangelder et al., 2012).

3 Methods

3.1 Field observations

From 1995 to 2013, the Royal Netherlands Institute for Sea Research (NIOZ) conducted shipboard monitoring of the Eastern Scheldt on a biweekly to monthly basis. The monitoring campaign routinely collected water samples at eight stations in the basin (OS1–OS8, Fig. 2) for nutrient measurements and filtered them for measurements of chl a and SPM. The dataset has been applied in several previous studies (Cloern and Jassby, 2010; Smaal et al., 2013), and the sampling method was the same as that described by Kromkamp and van Engeland (2010). The Dutch government agency Rijkswaterstaat (RWS) has monitored nutrients and chl a in the Eastern Scheldt at different locations (e.g., RWS1–RWS4, Fig. 2), and these monthly data are freely accessible on the RWS data portal (https://waterinfo.rws.nl, last access: 13 August 2020). Compared with the NIOZ data, the RWS data include two offshore stations: RWS1 and RWS2 (Fig. 2). Since the study region is mostly well mixed (Wetsteyn and Kromkamp, 1994), both datasets used surface samples to represent the water column at each station.

Figure 2 Average spring (March to May) chlorophyll-a (chl a) concentration from 1995 to 2013 at the NIOZ (OS1–OS8) and Rijkswaterstaat (RWS1–RWS4) monitoring stations. The error bars for each of the stations indicate standard deviations. The map shows the bathymetry in the GETM-FABM model domain denoted in Fig. 1 and marks all NIOZ and RWS stations within the domain. In this study, RWS1–RWS4 are short names for Walcheren 2 km, Walcheren 20 km, Wissenkerke, and Lodijkse Gat in the RWS database, respectively.

Primary production was estimated by ¹⁴CO₂ uptake (mg C h⁻¹) during 2 h incubation experiments for all 16 NIOZ sampling dates in 2010 at stations OS2 and OS8 (Fig. 2). Incubation experiments were conducted following a previously described method (Kromkamp and Peene, 1995; Kromkamp et al., 1995). A PI curve linking irradiance (I, $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) to the chl-a-normalized carbon fixation rate (P, $\mathrm{mg}\mathrm{C}\mathrm{mg}\mathrm{chl}{\mathit{\text{a}}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) was mathematically represented by a maximum carbon fixation rate (P_m, $\mathrm{mg}\mathrm{C}\mathrm{mg}\mathrm{chl}{\mathit{\text{a}}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$), an initial slope of the curve (α, $\mathrm{mg}\mathrm{C}\mathrm{mg}\mathrm{chl}{\mathit{\text{a}}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}(\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$), and the irradiance at which P_m occurred (I_opt, $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) according to the model of Eilers and Peeters (1988). Light intensity at multiple water depths was measured in the field with LI-COR LI-192SB cosine-corrected light sensors connected to a LI-COR LI-185B quantum radiometer photometer to estimate the light extinction coefficient K_d (m⁻¹) and generate the light attenuation curve $I_{\mathrm{z}}=I_{\mathrm{0}}\exp (-z{K}_{\mathrm{d}})$, where I₀ and z are surface light intensity and water depth, respectively. With the hourly photosynthetically active radiation (PAR) measured at the NIOZ as I₀, the hourly PAR (I_z) throughout the water column was computed. For a full description, see Kromkamp and Peene (1995). Then, with P_m, α, I_opt, and I available, the hourly photosynthetic rate at each water depth (P_z) was calculated and integrated over depth to obtain the primary production of the entire water column and during the whole day. We used the measured primary production data without estimating the respiratory losses as respiration will not affect the nitrogen content of the algae. Over a short incubation time, the ¹⁴C method is often thought to reflect gross primary productivity (GPP). However, results by Halsey et al. (2010, 2013) showed that even a 30 min ¹⁴C incubation experiment can reflect GPP at low growth rates and net primary productivity (NPP) at high growth rates. Hence, as growth rates are generally high during the main growing seasons (Underwood and Kromkamp, 1999), we assumed that our ¹⁴C method reflected NPP measurements. The phytoplankton turnover time (PT) was calculated by dividing the observed phytoplankton biomass by the NPP.

3.2 Numerical modeling

A three-dimensional hydrodynamic–biogeochemical model – the GETM-FABM (General Estuarine Transport Model coupled with the Framework for Aquatic Biogeochemical Models) – was applied in a 2-year (2009–2010) simulation to identify drivers of spatial phytoplankton dynamics in the Eastern Scheldt. The GETM and FABM are open-source models and are available from https://getm.eu/ and https://github.com/fabm-model/fabm (last access: 13 August 2020), respectively. The model was implemented on a 300 m × 300 m rectangular grid with 10 equally divided vertical layers, covering the Eastern Scheldt and part of the North Sea (Fig. 1). The hydrodynamic model using GETM version 2.5.0 was driven by realistic meteorological forcing (winds, irradiance, air pressure, etc.), and tides and the output water level, temperature, salinity, and current velocity were calibrated and validated with observational data (Jiang et al., 2019a). Jiang et al. (2019a) provided a detailed description of the GETM setup and model validation. The validation of FABM is presented in Sect. 4.2.

The biogeochemical model was coupled online with GETM on the FABM platform (Bruggeman and Bolding, 2014). The physical and biogeochemical simulations were conducted simultaneously with a time step of 8 s. At each time step, GETM provided FABM with the environmental variables, such as temperature, water elevation, and irradiance. The transport and mixing of nutrients, detritus, and plankton biomass was represented by the same equation as that for salinity except that phytoplankton and detritus sank at a speed of 0.2 and 1.0 m d⁻¹, respectively (Eppley et al., 1967; Soetaert et al., 2001).

Table 1 Formulations used in the biogeochemical model in this study. Parameters and variables in each equation are described in Table 2.

Table 2 The main variables (bold) and parameters (underlined, followed by values) in equations in Table 1. The parameter values are based on ranges found in prior literature (Soetaert et al., 2001; Jiang and Xia, 2017; Wijsman and Smaal, 2017) and were tuned for our application.

Our biogeochemical model was nitrogen-based and consisted of a pelagic and benthic module (Fig. 3). The pelagic module was a typical NPZD framework comprising the state variables nutrient (DIN, dissolved inorganic nitrogen), phytoplankton, zooplankton, and detritus (in mmol nitrogen m⁻³), while the benthic variables (which were in mmol nitrogen m⁻²) included benthic detritus, microphytobenthos, and the three dominant bivalve species in the Eastern Scheldt: mussels, oysters, and cockles. All mass transfer processes are shown using arrows in Fig. 3, and the main formulations, variables, and parameters for calculating these processes are summarized in Tables 1 and 2. The climatological data in December and January were averaged using the 19-year observations and were used as the initial model condition. The shellfish distribution (see Fig. 1) and annual biomass in 2009 and 2010 were estimated by Wageningen Marine Research (Smaal et al., 2013; Jiang et al., 2019a). The model output was compared with available observational DIN, chl a, and NPP, as described in Sect. 3.1. Given that NPP was measured using carbon-based methods, the nitrogen-based simulation results were converted to carbon using the Redfield ratio (C : N = 6.625 mol C mol N⁻¹). Phytoplankton biomass was measured in chlorophyll units. We prescribed a chl : N ratio (chl : N = 2 mg chl mmol N⁻¹, Soetaert et al., 2001) to compare our model output to the chlorophyll data.

Figure 3 Conceptual diagram of the nitrogen-based seven-variable biogeochemical model structure in FABM. Boxes and arrows denote state variables and fluxes of nitrogen, respectively.

In order to investigate the roles of coastal influx and benthic grazing in shaping the spatial phytoplankton patterns in the basin, we conducted two idealized numerical scenarios in addition to the realistic (baseline) run. One scenario was halving the DIN concentration and phytoplankton biomass at the open boundary (i.e., the North Sea including the Western Scheldt and Rhine river plumes, see Fig. 1). The other scenario switched off the bivalve state variables. Based on our assessment, the effects of freshwater input on DIN and chl a were minimal, local, and far less significant than the above two factors. Thus, the sensitivity runs of the freshwater input will not be elaborated in the following.

3.3 Satellite remote sensing

A clear-sky Sentinel-2 Multispectral Instrument (MSI; 10 m spatial resolution) satellite image of 11 May 2018 (10:55 UTC) for tile 31UET was downloaded as Level-1C data from the Copernicus Open Access Hub (https://scihub.copernicus.eu, last access: 13 August 2020). The ACOLITE processor (version Python 20190326.0) developed by the Royal Belgian Institute of Natural Sciences (RBINS; Vanhellemont and Ruddick, 2016) was applied using default settings to correct for atmospheric (aerosol) effects based on a dark spectrum fitting (Vanhellemont and Ruddick, 2018; Vanhellemont, 2019), to flag clouds and land, and to retrieve chl-a concentrations, using the red edge algorithm defined by Gons et al. (2002) with a mass specific chl-a absorption set to 0.015. Data were extracted in the Sentinel Application Platform (SNAP version 7.0.0) and converted to a GeoTIFF raster for further processing in ArcGIS. Georeferencing of the raster was enhanced using an affine transformation to a detailed topographic map of Rijkswaterstaat. The satellite image was acquired during high water: the water level at the Rijkswaterstaat tide gauge station of Stavenisse (https://waterinfo.rws.nl/#!/nav/index, last access: 13 August 2020) was +1.12 m NAP incoming tide during overpass. A Sentinel-2 MSI image of 21 April 2019 (10:56 UTC) was acquired during low water conditions (such as −1.58 m NAP incoming tide) and was processed in the same way. “Land” flags obtained from this low water image were used to further flag shallow waters (i.e., the inundated tidal flats) in the high water image in order to avoid potential bottom reflectance.

4 Results

4.1 Field observations

The 19-year chl-a time series illustrates the seasonal pattern of phytoplankton biomass in the Eastern Scheldt (Fig. 4). The spring bloom takes place in March or April during conditions of increasingly favorable temperature, light, and nutrients and lasts less than a month. The peak biomass varies dramatically interannually (Fig. 5), with smaller peaks during different months, especially in 2010 (Fig. 4). Likely due to nutrient limitation and grazing pressure, the summer biomass stays mostly below 10 mg m⁻³ (Figs. 4, 5). In the well-mixed Eastern Scheldt with limited freshwater input, temperature and light constrain algal growth in winter, when nutrients accumulate and phytoplankton biomass falls below 3 mg m⁻³ (Figs. 4, 5). These seasonal controls of phytoplankton variability are substantiated with the numerical model in Sect. 4.2.

Figure 4 Time series of chlorophyll-a (chl a) concentrations from 1995 to 2013 at (a) NIOZ stations OS1, OS3, and OS8 and (b) Rijkswaterstaat stations RWS1–RWS4. The intervals between grid lines indicate 2 months. See Fig. 2 for station locations.

Figure 5 Monthly average chlorophyll-a (chl a) concentrations from 1995 to 2013 at (a) NIOZ stations OS1, OS3, and OS8 and (b) Rijkswaterstaat stations RWS1–RWS4. The error bars of each station indicate standard deviations. See Fig. 2 for station locations.

To better display the spatial chl-a gradients, monthly averages and standard deviations of the 19-year observations are displayed (Fig. 5). A decreasing chl-a gradient from the mouth (OS1) to head (OS8) of the basin is observed mainly during the spring bloom, whereas the spatial phytoplankton gradient is not as pronounced in summer and winter (Fig. 5a). The station RWS1 that is close to the mouth of the Western Scheldt Estuary usually has a higher chl-a concentration than further offshore (RWS2) and in the Eastern Scheldt (RWS3 and RWS4; Fig. 5b). Despite interannual variability in the timing of the bloom and the different sampling time every year, the period from March to May mostly covers the initiation, development, and wane of the spring bloom. The 19-year average phytoplankton biomass during this season demonstrates a clear gradient in the bay and adjacent coastal sea (Fig. 2). The chl-a concentration decreases from the Western Scheldt plume region (RWS1) offshore (RWS2) and further into the eastern and northern ends of the Eastern Scheldt (Fig. 2).

4.2 Numerical modeling

The model results compared to the observed concentrations of DIN and chl a from a 2-year simulation are shown in Fig. 6. Most DIN consumption occurs during the spring bloom, and the regenerated DIN accumulates over winter until the next bloom sets off. The simulated chl a during the bloom demonstrates the same gradient between the western and eastern bay as observed (OS1 > OS3 > OS8, Fig. 6d–f). The model skill is quantified by the r² value (r²=0.89 for DIN and 0.66 for chl a) and the root mean square errors (RMSE = 6.0 mmol m⁻³ for DIN and 3.9 mg m⁻³ for chl a) between the simulation and observation. Despite capturing the major seasonal and spatial patterns, the model seems to miss some details, such as overestimating the recycled DIN at OS8 and showing a slower collapse of spring blooms than observed. Meanwhile, the daily time series of the model output exhibits spring–neap biweekly fluctuations (Fig. 6) that cannot be substantiated by the observations owing to a low sampling frequency.

Figure 6 Comparison between simulated and observed dissolved inorganic nitrogen (DIN) and chlorophyll a (chl a) in the years 2009 and 2010 at stations OS8, OS3, and OS1. See Fig. 2 for station locations.

The modeled NPP, the product of the phytoplankton biomass and growth rate, is in general agreement with the measurements (Fig. 7, black line). According to Eq. (2) in Table 1, the growth rate is a function of temperature, nutrient, and light factors in the model. Here, we decompose the seasonal cycle of these three factors and use their product to assess the growth rate (Fig. 8). Before the bloom, both modeled phytoplankton biomass and growth rates are low, resulting in a low NPP (Figs. 7, 8). The fast-growing period, around Day 100 as a consequence of increased temperature and light (Figs. 7, 8), triggers the increase in simulated biomass that leads to the bloom. The spring bloom is terminated due to enhanced nutrient limitation around Day 125 in the model (Figs. 7, 8). In the low-biomass post-bloom summer (Fig. 6), both modeled and measured NPP is only slightly lower than that in the bloom (Fig. 7) and the environmental conditions are still favorable for growth (Fig. 8). The summer growth rate is mainly fueled by regenerated nutrients, while the low biomass primarily results from grazing. The model underestimated DIN and, thus, NPP after the spring bloom (days 490–540 in Fig. 6 and days 125–175 in Fig. 7) and overestimated the recycled DIN at OS8 in fall 2010 (days 600–650 in Fig. 6a), which explains the overestimation of NPP during this period (days 235–285 in Fig. 7a). The observed and simulated NPP at OS8 (902.6 ± 928.4 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 1033.9 ± 1084.3 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively) is generally higher than that at OS2 (722.5 ± 794.6 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and 606.0 ± 499.5 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, respectively). According to a one-tailed t test, the difference between the observed NPP values at these two stations is not significant (t=0.59, p>0.05, n=16), whereas the simulated NPP value at OS8 is significantly higher than that at OS2 (t=6.85, p<0.05, n=365) due to the overestimation. This is in contrast to the chl a, which is higher at OS2 (Fig. 2). Based on the observed chl a and measured NPP, PT (March to October) is 0.92–5.2 d and 0.13–4.3 d during the warm months at OS2 and OS8, respectively.

Figure 7 Comparison between modeled and measured depth-integrated net primary production (NPP, $\mathrm{mg}\mathrm{carbon}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in 2010 at stations OS8 and OS2. The three model scenarios include the baseline scenario, halving the open boundary nutrient and phytoplankton loading, and switching off bivalve filtration feeders.

Figure 8 Modeled surface temperature, nutrient, and light factors affecting the phytoplankton growth rate at (a) OS8 and (b) OS2. The temperature factor is calculated according to Eq. (1) in Table 1. The nutrient factor = $\text{DIN}/(\text{DIN}+\text{ksDIN})$, and the light factor = $\text{PAR}/(\text{PAR}+\text{ksPAR})$, as shown in Eq. (2) in Table 1. According to Eq. (2), the product of these three factors is positively related to the phytoplankton growth rate.

The calibrated and validated model was used to map the 15 d average chl-a concentration during the peak bloom in 2009 (Fig. 9). The North Sea exhibits significantly higher algal biomass than the Eastern Scheldt (Fig. 9). Inside the bay, phytoplankton biomass is clearly low over the shellfish-colonized area (cf. Figs. 1 and 9). The north–south and east–west chl-a gradients observed in field monitoring data are reproduced in the model results (Figs. 2, 9).

Figure 9 The 15 d (5–19 March) average of modeled chlorophyll a (chl a) during the peak spring bloom in 2009. Gray squares indicate the locations of wild and cultured shellfish, as in Fig. 1.

When switching off bivalve activities, the modeled phytoplankton biomass significantly increases, especially at the eastern station OS8 (Fig. 10). At this station, the chl a during the bloom is nearly tripled, whereas it doubled at OS3 and increased by 20 % at OS1, respectively (Fig. 10). The east–west spatial chl-a gradient is weakened in spring and even reversed in summer, i.e., concentrations decrease seaward (Fig. 10). Remarkably, without bivalves, the summer NPP at OS8 is not greatly affected (Fig. 7a) despite increased algal biomass, which implies a reduction in the growth rate (Eq. 2, Table 1). Given the unchanged light and temperature in the bivalve-free scenario, the reduced growth rate results from diminished nutrient regeneration. The summer NPP at OS2 is increased when bivalves are turned off (Fig. 7b), which is a consequence of increased phytoplankton biomass.

Figure 10 Modeled and observed chlorophyll a (chl a) in 2009 at stations OS8, OS3, and OS1. The three model scenarios include the baseline scenario, halving the open boundary nutrient and phytoplankton loading, and switching off bivalve filtration feeders. See Fig. 2 for station locations.

Halving the DIN and phytoplankton loading from the North Sea hardly has any influence on the NPP in the Eastern Scheldt model (Fig. 7). This indicates that allochthonous coastal nutrients are not a major source of inner-bay primary production, which relies mainly on recycling. With halved coastal import, the modeled peak phytoplankton biomass is nearly halved at OS1, but the reduction is lower at OS3 (∼ 35 %) and OS8 (∼ 20 %); see Fig. 10. Therefore, tidal import has its modeled impact almost exclusively near the bay mouth during the bloom. Thus, tidal import contrasts with the benthic bivalves that appear to exert grazing pressure all over the bay. Bivalves also stimulate primary production by replenishing inorganic nutrients into the water column, which is crucial in nutrient-depleted seasons.

4.3 Satellite remote sensing

Remote sensing images with sufficient spatial resolution, in this case the Sentinel-2 MSI data, are utilized to complement the spatial patterns shown in observational and modeling data. In an attempt to find images during the spring bloom and high tide (to avoid interference from bottom reflectance), we only found one post-bloom snapshot under clear-sky conditions (Fig. 11a). This provides additional insight into the observed and modeled spatial chl-a pattern. On 11 May 2018, the chl-a concentration was highest in the central basin and reduced eastward and northward into the highly bivalve-populated areas (Figs. 1, 11a), which is consistent with the chl-a gradient described in Sect. 4.1 and 4.2. However, the post-bloom chl-a concentration was low in the North Sea so that low import was shown in the southwestern bay near the mouth (Fig. 11a). Such a spatial chl-a pattern with higher concentrations in the central basin (Fig. 11a) is often present in the model results. For example, in the post-bloom period in 2008 and 2009, the chl-a concentration at OS3 is higher than OS1 and OS8 at times (Fig. 5a). Likewise, in a post-bloom model snapshot during high tide on 1 May 2010 (Fig. 11b), the phytoplankton distribution exhibits a similar spatial gradient to that in Fig. 11a (r²=0.118, n=2817, p<0.0001).

Figure 11 Chlorophyll a (chl a) in the Eastern Scheldt retrieved from (a) a high tide Sentinel-2 MSI image of 11 May 2018 at 10:55 UTC masking tidal flats from a low tide Sentinel-2 MSI image of 21 April 2019, and (b) the model on 1 May 2010 at 17:00. Both snapshots are during high tide. The coordinate system in (a) is Amersfoort/RD New.

5 Discussion

The approaches applied in this case study, including field observation, numerical modeling, and satellite remote sensing, each have their drawbacks. The monitoring data are not frequent enough to capture the peak bloom that lasts only a couple of weeks and misses details in the spatial distribution between stations. The temporal resolution of the satellite data is even lower, but the spatial detail is very high. The model, while resolving spatial and temporal scales at a high resolution, is based on simplified assumptions. The NPZD model considers nitrogen only and assumes no phosphorus or silicon limitation in phytoplankton growth. In late spring, phosphorus or silicon may become limited in the Eastern Scheldt (Wetsteyn and Kromkamp, 1994; Smaal et al., 2013). This likely explains the faster DIN consumption in the simulated data compared with the observations (Fig. 6a–c), resulting in the accelerated nitrogen limitation and underestimation of the post-bloom NPP (Figs. 7, 8). Additionally, our model does not account for the shellfish harvest, mostly in late summer, which can contribute to the overestimation of the regenerated DIN and, hence, the NPP, especially in the eastern region (e.g., Figs. 6a, 7a). When converting the nitrogen-based model results for comparison with chl a and the carbon-based NPP, the chl : N and Redfield ratios were applied, without taking the variable stoichiometry in natural phytoplankton groups into account. A model accounting for the phytoplankton physiological plasticity (Faugeras et al., 2004), e.g., low (high) chl-a content under high (low) light intensity and high C : N ratios under nitrogen limitation, should be considered in further studies. Despite these simplifications and limitations, the approaches complement each other in the spatiotemporal resolution and coverage and offer insight into the phytoplankton distribution in the Eastern Scheldt, as well as the underlying mechanisms.

Grazing by filtration feeders is found to be the dominant factor shaping the spatial and seasonal phytoplankton patterns in the Eastern Scheldt. In the eastern and northern bay, RT is relatively long (> 100 d, calculated using the e-folding method and the remnant function by Jiang et al., 2019b), water depth and cross-bay area are much smaller (Fig. 12c), and the tidal amplitude and mixing are stronger due to geometric convergence (Jiang et al., 2020), which contributes to stronger pelagic–benthic coupling and creates favorable feeding conditions for suspension feeders (Hilly, 1991). Thus, over and near the shellfish habitat, the phytoplankton biomass is usually low, even during the bloom (Figs. 1, 9). Smaal et al. (2013) attributed the decline in the annual primary production and chl-a concentration in the Eastern Scheldt to overgrazing, as found in the Bay of Brest (Hilly, 1991) and many Danish estuaries (Conley et al., 2000). Our findings support the predominant top-down control on phytoplankton distribution and standing stocks (Fig. 10), as well as on primary production, particularly in the post-bloom seasons (Fig. 7b). It has been shown that a recruitment failure of mussels and cockles promotes primary production and algal accumulation in the Dutch Wadden Sea (Beukema and Cadée, 1996), consistent with our numerical experiment removing bivalves. Although bivalves accelerate nutrient remineralization, this positive feedback on phytoplankton growth does not compensate for the grazing loss (Figs. 7, 8, and 10). Optimization of the bivalve stock size and culture locations based on these scientific insights could enhance phytoplankton proliferation and increase the shellfish carrying capacity.

Figure 12 The 15 d (5–19 March) average of modeled chlorophyll a (chl a) during the peak spring bloom in 2009 (b) along a transect over the southern channel of the Eastern Scheldt and (c) cross-sectionally integrated from the bay mouth landward. The transect location is shown in (a). The blue line in (c) denotes the cross-sectional area used for integration. The northern branch (a) is excluded from the calculation in (c) due to a different orientation of channels. The distance on the x axis of panels (b) and (c) is from west to east. The three model scenarios include the baseline scenario, halving the open boundary nutrient and phytoplankton loading, and switching off bivalve filtration feeders.

The strong top-down control by shellfish can result from cultured or natural populations. In the Eastern Scheldt, the high shellfish biomass consists of both wild cockles and oysters as well as commercial Pacific oysters and blue mussels (see the data in Jiang et al., 2019a). Worldwide, shellfish culture can be an important source of benthic grazing. Examples include the farmed oyster in Willapa Bay (Banas et al., 2007) and cultivated mussels in the Baie des Veys Estuary (Grangeré et al., 2010). Invasive shellfish species can also exert significant grazing pressure. After their introduction, the invasive clam Potamocorbula amurensis in San Francisco Bay (Lucas et al., 2016) and the exotic dreissenid mussels (Dreissena spp.) in the Hudson River (Strayer et al., 2008) and Laurentian Great Lakes (Higgins and Vander Zanden, 2010) quickly dominated the benthos, strongly increased the filtration capacity, and, thus, extensively changing the lower food web. Also in the Eastern Scheldt, it has been noted that great care should be taken with regards to invasive shellfish, such as Ensis americanus (Smaal et al., 2013).

Ecologists use three timescales to assess the carrying capacity of shellfish-dominated ecosystems: clearance time (CT), the time it takes for suspension feeders to filter the entire basin; RT; and PT (Dame and Prins, 1998). A CT : RT ratio of less than one suggests that the rate at which the system is replenished is outpaced by the filtration rate, and that the pelagic ecosystem is controlled by benthic grazing. A CT : PT ratio less than one reveals that the food (phytoplankton) reproduction is slower than filtration so that the system may collapse. In the Eastern Scheldt, the grazing pressure is immense (CT : RT < 1, CT = 19.6 d, average RT > 30 d, Jiang et al., 2019a), although the system is still sustainable (CT : PT > 1). Thus, the strong benthic filtration capacity consumes considerable pelagic production and puts high pressure on the pelagic food web (Smaal et al., 2013). Compared with other estuarine–coastal systems used for shellfish farming (e.g., the western Wadden Sea, Beatrix Bay, Narragansett Bay, etc.), the indices in the Eastern Scheldt suggest relatively high exploitation, including a high ratio of the overall shellfish biomass to basin volume, i.e., the relative shellfish density, and low CT : RT and CT : PT ratios (Jiang et al., 2019a; Smaal and van Duren, 2019).

Compared with benthic feeding, tidal import mainly influences the phytoplankton biomass in the southern channel near the mouth during the spring bloom. The Southern Bight of the North Sea, more specifically the water in river plumes (e.g., the Western Scheldt and Rhine plumes), is relatively productive compared with other shelf seas (van der Woerd et al., 2011). The spring bloom in the Eastern Scheldt is usually not as strong as that in the adjacent North Sea (Figs. 2, 5, and 9); therefore, the tidal import of phytoplankton from the North Sea sets the upper limit of the phytoplankton in the Eastern Scheldt (Fig. 10). However, in other seasons, phytoplankton biomass is very similar inside and outside of the bay, and coastal import is not as high as in spring (Fig. 10). In spring, the measured PT is shorter than 5.2 d (Sect. 4.2), which is comparable to the RT near the bay mouth (Jiang et al., 2019b) so that the phytoplankton biomass import is of similar magnitude to local production in the southwestern Eastern Scheldt. The role of tidal import decreases further into the bay, as supported by the numerical experiment shown in Fig. 12, and is minimal in the eastern bay (e.g, at OS8, Fig. 10a). Thus, the effect of tidal import on the phytoplankton biomass in the Eastern Scheldt is subject to seasonal and spatial variability, and it depends on two conditions: (1) significantly different phytoplankton biomass inside and outside of the bay and (2) a sufficiently short RT compared with the phytoplankton turnover time.

While accounting for only 45.0 % of the basin volume excluding the northern branch, the western part of the Eastern Scheldt, defined as 0–10 km east to the bay mouth, contains 60.0 % of the phytoplankton biomass during the peak spring bloom (data in Fig. 12c). Therefore, the impact of halving tidal import in the model, which is most pronounced in the western Eastern Scheldt (Fig. 12c), has a strong influence on the overall phytoplankton standing stocks, reducing the total phytoplankton biomass by 38.5 % during the peak bloom, from 58.7 to 36.1 t.

Note that the observed seaward increasing phytoplankton biomass is not as pronounced in each of the 19 years, such as in 1997 and 2012 in Fig. 4a. The interannual variability of the spring bloom timing and magnitude is among the potential causes of this difference. In winter and early spring, the shallower and landlocked Eastern Scheldt may be warmed up a few days faster than the adjacent North Sea, which may result in earlier spring bloom at the landward end in some years, such as 1996, 1997, and 2005 (Fig. 4b). Thus, the resultant spatial chl-a distribution may differ if the sampling activity is conducted during or between these two bloom windows. Given the relative low sampling frequency (every month or 2 weeks), different observational activities may change the spatial chl-a distribution (e.g., the year 2010 in Fig. 5a and b). Additionally, if the bloom is of similar magnitude inside and outside of the bay, the spatial gradient may not be as strong, such as in 2012 in Fig. 4a. Therefore, in consideration of the interannual variability of the spring blooms and the possible mismatch with sampling campaigns, the observational spatial phytoplankton gradient may be inconspicuous in certain years but evident over the long term (Figs. 2, 5).

6 Synthesis

The Eastern Scheldt represents a land–ocean transitional system that is shallow, dynamic, and driven by pelagic–benthic coupling and exchange with the sea. The grazing pressure increases into the bay because of reduced water depth and increasing bivalve biomass, RT, tidal mixing, and, thus, pelagic–benthic coupling, while the North Sea with its higher phytoplankton biomass is a phytoplankton source. As a result, the phytoplankton concentration during the spring bloom consistently declines from the seaward to the landward end. When halving the nutrient and algal loading from the North Sea, the phytoplankton gradient in spring is not as pronounced, although it is still decreasing toward the landward end (Fig. 12). Without the grazing sink, however, the phytoplankton distribution tends to be spatially uniform (Fig. 12). Given the temporal variability of dominant environmental factors, the phytoplankton gradient also changes over time. In the post-bloom period for instance, chl a may exhibit a central maximum, or it may exhibit a constant concentration in winter. This shows that the spatial gradient of phytoplankton biomass in estuarine–coastal systems depends on the relative importance of the main drivers of phytoplankton accumulation. Therefore, we compare the Eastern Scheldt case with other estuarine–coastal systems, including different spatial phytoplankton patterns and their reported dominant environmental drivers (Fig. 13).

Figure 13 Common spatial patterns of phytoplankton biomass in estuarine–coastal systems. For comparison with the Eastern Scheldt, example ecosystems for each type are given along with references, the dominant flushing mechanisms, and main drivers of phytoplankton accumulation.

Increasing phytoplankton biomass from the landward to the seaward ends (Fig. 13) can often be ascribed to an increasing source from the seaside, an increasing sink landward, or more favorable growth conditions towards the sea. The Eastern Scheldt is a typical example during the spring bloom (Fig. 8), shaped by both marine input and increasing benthic filtration landward. The seaward increasing gradient is also common in estuaries and bays open to coastal upwelling zones (e.g., the Rías Baixas of Galicia and Tomales Bay), where algal blooms generated during upwelling events are transported into bays via multiple physical mechanisms including tidal stirring and gravitational and wind-driven circulation (Figueiras, et al., 2002; Hickey and Banas, 2003; Martin et al., 2007). In contrast, phytoplankton in the Chilika Lagoon is mostly light-limited due to the massive riverine sediment loading. Here, it is a seaward increase in water transparency that leads to increasing chl-a concentrations (Srichandan et al., 2015). Hence, similar phytoplankton gradients may be driven by distinct mechanisms in various systems.

Contrary to the Eastern Scheldt case, phytoplankton biomass decreases in the seaward direction in some systems (Fig. 13). A typical example is the Scheldt River and Western Scheldt Estuary, a eutrophic and turbid estuary with salinity ranging from 0 to 30 (Soetaert et al., 2006). Numerical models and field observation reveal that the chl-a concentration is highest in the tidal fresh portion, reduces sharply between a salinity of 5 and 10, and is maintained at a lower level towards the polyhaline region (Soetaert et al., 1994; Kromkamp et al., 1995). The high phytoplankton biomass in the upper reach is a result of tributary import, high nutrient levels, and a lack of zooplankton grazers, whereas the increasing salinity stress on the freshwater species and grazing pressure in the mesohaline zone suppress the phytoplankton proliferation (Soetaert et al., 1994; Kromkamp and Peene, 1995; Muylaert et al., 2005). A similar seaward decreasing phytoplankton gradient is also found in many river and estuarine plumes (Fig. 13), where the nutrient gradient controls the phytoplankton distribution (Gomez et al., 2018; Jiang and Xia, 2018).

A chl-a maximum zone (CMZ) occurs in many estuaries with substantial freshwater input (Fig. 13). Taking Chesapeake Bay as an example, the upper bay is characterized by high terrestrial sediment concentrations and strong light limitation for phytoplankton growth (Son et al., 2014). A turbidity maximum zone is located near the front of salt intrusion (North et al., 2004). Beyond this location, the CMZ appears in the middle reach (Jiang and Xia, 2017), while nitrogen limitation is constantly detected in the lower bay (Miller and Harding, 2007). The CMZ is a combined consequence of the optimal light conditions and abundant terrestrial nutrients, and the CMZ location and coverage shift with river discharge and weather (Fisher et al., 1988; Miller and Harding, 2007). In some other estuaries with a CMZ (e.g., the Neuse–Pamlico Estuary and York River), owing to a narrow river channel and high discharge, the flushing rate in the upper estuary can be faster than the phytoplankton turnover rate, which, rather than light, limits phytoplankton accumulation (Sin et al., 1999; Valdes-Weaver et al., 2006). In these systems, the CMZ is always in wider reaches with sufficiently long RTs (Valdes-Weaver et al., 2006).

If the transport loss is higher than the growth rate throughout the basin, the phytoplankton biomass is low everywhere and is negatively correlated with the flow velocity (Fig. 13). The Hudson River estuary is one such estuary, with high nutrient loading but low and hardly spatially variable chl a (Howarth et al., 2000). After the colonization of the invasive zebra mussel (Dreissena polymorpha) in the 1990s, grazing and transport losses have been two dominant sink terms maintaining a low basin-wide phytoplankton standing stock in the estuary (Strayer et al., 2008). Similarly, due to the nonindigenous P. amurensis, the San Francisco Bay witnessed a five-fold drop in chl a and the suppression of zooplankton, as well as higher trophic levels (Cloern and Jassby, 2012; Lucas et al., 2016). The chl-a distribution has shifted from maximizing at the middle bay to being spatially uniform (Cloern et al., 2017), which is generally associated with strong sink factor(s) distributed throughout the system.

The dominant sink (or source) factor is not always distributed uniformly nor does it follow consistent gradients in estuarine–coastal systems, generating an irregular phytoplankton distribution (Fig. 13). For instance, in the Baie des Veys Estuary, benthic grazing by cultivated oysters results in an area of low chl-a concentrations over the oyster bed, and this patch of low chl a is imposed onto a seaward decreasing chl-a gradient, forming an irregular spatial pattern (Grangeré et al., 2010). In the Krka Estuary, an untreated sewage discharge acts as a DIN point source, increasing the phytoplankton production downstream. Without the point source, phytoplankton seem to decrease seaward (Ahel et al., 1996). In the St. Lucia Estuary, controls of primary production include nutrient stoichiometry, temperature, irradiance, and hydrological changes which all vary in different subregions and render complex spatial heterogeneity in phytoplankton distribution (van der Molen and Perissinotto, 2011).

7 Summary

In the Eastern Scheldt, a tidal bay along the North Sea, we detect a seaward increasing phytoplankton gradient in spring in the 2-decade monitoring data. This spatial chl-a pattern was also reproduced with a nitrogen-based NPZD model calibrated and verified using observational data. In an effort to understand the main drivers of such a phytoplankton gradient, two experimental model runs were performed: switching off bivalve filtration and halving the nutrient and phytoplankton concentrations in the North Sea boundary, respectively. Results indicate that the landward increasing benthic grazing pressure is the primary cause of the spatial phytoplankton gradient, while import from the North Sea tends to strengthen the gradient. The satellite image implies that tidal import is mainly influential in the southwestern bay. With the variation of these two drivers, the spatial phytoplankton distribution varies seasonally.

In a synthesis of the literature, we compared the Eastern Scheldt with other estuarine–coastal systems focusing on how the spatial phytoplankton gradient is shaped by the distribution of the main environmental drivers. Common spatial phytoplankton patterns include seaward increasing, seaward decreasing, concave with a chlorophyll maximum, weak spatial gradients, and irregular patterns. It should be noted that the spatial phytoplankton pattern is subject to temporal changes and cannot be discussed without specifying the temporal window. For example, the spatial chl-a gradient in this study is different before, during, and after the spring bloom and is subject to substantial interannual variability; in river-dominated systems (Fig. 13), the phytoplankton distribution is usually regulated by the episodic, seasonal, and interannual variations of river discharge (Kromkamp and van Engeland, 2010). In addition to natural changes, phytoplankton abundance and spatial heterogeneity can reflect how the lower trophic levels are affected by anthropogenic influences and stressors, such as aquaculture (this study), invasive species (Cloern et al., 2017), and coastal engineering works (Wetsteyn and Kromkamp, 1994).