Influences of Nutrient Sources on the Alternation of Nutrient Limitations and Phytoplankton Community in Jiaozhou Bay, Southern Yellow Sea of China

Abstract

A marine ecosystem box model was developed to reproduce the seasonal variations nutrient concentrations and phytoplankton biomasses in Jiaozhou Bay (JZB) of China. Then, by removing each of the external sources of nutrients (river input, aquaculture, wastewater discharge, and atmospheric deposition) in the model calculation, we quantitatively estimated its influences on nutrient structure and the phytoplankton community. Removing the river input of nutrients enhanced silicate (SIL) limitation to diatoms (DIA) and decreased the ratio of DIA to flagellates (FLA); removing the aquaculture input of nutrients decreased FLA biomass because it provided less dissolved inorganic nitrogen (DIN) but more dissolved inorganic phosphate (DIP) as compared to the Redfield ratio; removing the wastewater input of nutrients changed the DIN concentration dramatically, but had a relatively weaker impact on the phytoplankton community than removing the aquaculture input; removing atmospheric deposition had a negligible influence on the model results. Based on these results, we suppose that the change in the external nutrients sources in the past several decades can explain the long-term variations in nutrient structure and phytoplankton community. Actually, the simulations for the 1960s, 1980s, and 2000s in JZB demonstrated the shift of limiting nutrients from DIP to SIL. A reasonable scenario for this is the decrease in riverine SIL and increase in DIP from aquaculture that has reduced DIA biomass, promoted the growth of FLA, and led to the miniaturization of the phytoplankton.

1. Introduction

Coastal areas are regions of transition from terrestrial to ocean environments. The primary production is often high, and the transfer of this productivity through the food web could yield prolific fisheries [1,2]. The growth of phytoplankton in coastal areas can be limited by several factors, including light, temperature, and nutrients [3,4]. Among these, nutrients have attracted the most attention because their input can be human-controlled to reverse eutrophication [5,6,7].

Nutrient elements such as nitrogen (N), phosphorus (P), and silicon (Si) have been reported to be limitation factors for phytoplankton growth and their element ratio may affect community composition [8,9,10,11]. Usually, there are several external sources of nutrients in coastal waters, such as river inputs, exchanges with the marginal seas, and atmospheric deposition. The nutrient loads from these sources are different not only in the amount but also the element ratios, both of which have dramatically changed in past decades under the influences of anthropogenic activities and climate changes [12]. The changes in external sources of nutrients naturally regulate phytoplankton growth, biomass, and species composition, causing fluctuations in the higher trophic species [13], which have been already observed in coastal seas [14,15].

Jiaozhou Bay (JZB) is a semi-enclosed and well-utilized bay in the southwestern Yellow Sea (Figure 1). It has an area of 390 km² and an average water depth of approximately 7 m [16]. Many observations in JZB have documented the basic characteristics of the seasonal, interannual, and long-term variations in nutrient concentrations and phytoplankton biomasses [17,18,19,20,21]. It is suggested that the important factor controlling the phytoplankton growth in JZB is nutrients [22]. By comparing the nutrient loads and the standing stocks of nutrients in the bay, previous studies indicated that external nutrient sources had decisive effects on the multi-time scale variations of nutrient structures [19,23,24]. However, most of the observations focused on the bay itself, with limited quantitative information about the effects of different nutrient sources on the phytoplankton community composition.

There are several external nutrient sources for JZB. JZB receives fresh water and nutrients from more than ten rivers, among which Dagu River is the largest, accounting for 84% of the total freshwater discharge into the bay [25]. Three of the sewage treatment plants in Qingdao City discharge into JZB directly, i.e., the Haibohe, Licunhe, and Tuandao plants [23]. The fish bait and the excretion of aquaculture species are sources of dissolved inorganic nitrogen (DIN) and dissolved inorganic phosphate (DIP), which can accelerate the eutrophication [26]. Another source of nutrients in JZB is the atmospheric deposition [27]. In addition, the nutrient concentrations are generally higher in JZB than in the Yellow Sea (YS), and the water exchange with the YS removes nutrients from JZB [28]. Many studies calculated the loads of nutrients from these external sources and tried to estimate their contributions to primary production and phytoplankton compositions [23,24]. However, assumption is always necessary in such discussion, because field observation essentially cannot distinguish the nutrients from different sources.

Figure 1. The location and water depths (in meters) of Jiaozhou Bay. The estuary of Dagu River is marked by a star. The observations, including concentrations of three nutrient elements and chlorophyll-a used to validate model results, are shown. Stations of monthly observations from June 2003 to May 2004 carried by [17] are marked by diamonds. Stations of the monthly observations from January 2003 to December 2003 carried by [18] are marked by triangles. Jiaozhou Bay Marine Ecosystem Research Stations from 2004 to 2008 carried by [19] and from January 2011 to December 2011 carried by [29] are marked by dots.

Figure 1. The location and water depths (in meters) of Jiaozhou Bay. The estuary of Dagu River is marked by a star. The observations, including concentrations of three nutrient elements and chlorophyll-a used to validate model results, are shown. Stations of monthly observations from June 2003 to May 2004 carried by [17] are marked by diamonds. Stations of the monthly observations from January 2003 to December 2003 carried by [18] are marked by triangles. Jiaozhou Bay Marine Ecosystem Research Stations from 2004 to 2008 carried by [19] and from January 2011 to December 2011 carried by [29] are marked by dots.

Over the last half-century, the external sources of nutrients for JZB have changed under the anthropogenic perturbations [20,21,30]. The diminished river runoff due to the construction of water conservancy projects, the increase in industrial and domestic wastewater discharge, the heavy aquaculture activities, together, have had great impacts on the nutrient structures in JZB over the past few decades [31]. The unbalanced nutrient structure induced serious eutrophication and high frequent occurrences of the harmful algal bloom in JZB [24]. In addition, more DIN and DIP loads but less SIL load into JZB were suggested to decrease the biomass of large-size phytoplankton (DIA) and increase the biomass of small-size phytoplankton (FLA) [19,32]. The miniaturization of the phytoplankton community subsequently hinders energy transport to the higher trophic levels. Some toxic FLA have also caused the collapse of fisheries, along with the loss of several fish species [32,33]. To quantitatively understand the relationship between decadal shifts of nutrient loads of each external source and the nutrient structures and phytoplankton community inside the bay, the correlation analysis used by previous studies is not sufficient, because we need more information to understand the mechanism behind such phenomena.

Numerical modeling is an effective way to reveal the mechanisms linking the external loads of nutrients and the change in nutrient structure and phytoplankton community inside the bay. There have been some related numerical studies [28,34,35]. However, they usually considered one group of phytoplankton (DIA) and failed to reproduce the structure of phytoplankton community. Few studies have considered the role of each external source of nutrients. Therefore, it is still unclear how each external source of nutrients influences the variation in the phytoplankton community of JZB. Considering the small size of the bay and the vertically homogeneous distribution of physical and biogeochemical variables, we choose a box model in this study to quantify the influences of each external source of nutrients on the nutrient concentration and phytoplankton community. The understanding on the temporal changes over the entire bay is also helpful to understand the observations over the bay in the past several decades.

In this study, based on the published estimation of the nutrient loads from various external sources, we establish an ecosystem box model including three kinds of nutrients and two groups of phytoplankton with different cell sizes to quantify the influences of each external source of nutrients on nutrient structure and variation in the phytoplankton community, indicated by the ratio of DIA/FLA. The nutrient limitation conditions for three decades (the 1960s, the 1980s, and the 2000s) are compared to explain the effects of anthropogenic activities on the unbalanced nutrient structure and miniaturization of phytoplankton community.

2. Methods

2.1. Model Descriptions

The state variables in this box model included three kinds of nutrients (DIN, DIP, and SIL), two types of phytoplankton (DIA and FLA), and two types of detritus (DSi comprising Si and DNSi not comprising Si). Zooplankton grazing on phytoplankton could be taken into account by increasing the mortality rate of phytoplankton in the numerical model. The biogeochemical processes (Figure 2) among state variables have been studied previously and used in many ecosystem models [36,37]. DIA absorbs DIN, DIP, and SIL for photosynthesis, and releases three types of nutrients through respiration. The mortality of DIA contributes to DSi, which returns to the three types of nutrients through mineralization. The biogeochemical processes of FLA are similar to those of DIA. However, FLA grows depending on DIN and DIP rather than SIL. The mineralization of DNSi releases DIN and DIP. Sources and sinks for different state variables are provided in Appendix A. The parameters used in this model were based on previous studies [28,37] and are given in

Table 1. Parameters used in the model.

Parameter	Value	Unit	Description of Parameter
${\mathrm{Pmax}}_{\mathrm{dia}}$	1.3	${\mathrm{day}}^{-1}$	Maximum growth rate of DIA
${\mathrm{Pmax}}_{\mathrm{fla}}$	0.3	${\mathrm{day}}^{-1}$	Maximum growth rate of FLA
${\mathrm{resp}}_{\mathrm{bas}}$	0.138	-	Percentage of basic respiration
${\mathrm{resp}}_{\mathrm{pho}}$	0.05	-	Percentage of photorespiration
${\mathrm{r}}_{\mathrm{mort}}$	0.12	${\mathrm{day}}^{-1}$	Phytoplankton mortality rate
$\mathrm{al}\_\mathrm{temp}1$	0.055	${°\mathrm{C}}^{-1}$	Temperature-dependent growth rate
$\mathrm{al}\_\mathrm{temp}2$	0.054	${°\mathrm{C}}^{-1}$	Temperature-dependent respiration rate
${\mathrm{Iopt}}_{\mathrm{dia}}$	120	${\mathrm{W}\mathrm{m}}^{-2}$	Optimum light intensity of DIA
${\mathrm{Iopt}}_{\mathrm{fla}}$	70	${\mathrm{W}\mathrm{m}}^{-2}$	Optimum light intensity of FLA
${\mathrm{Kdin}}_{\mathrm{dia}}$	1.4	${\mathsf{\mu }\mathrm{molN}\mathrm{L}}^{-1}$	DIN half saturation constant for DIA
${\mathrm{Kdip}}_{\mathrm{dia}}$	0.12	${\mathsf{\mu }\mathrm{molP}\mathrm{L}}^{-1}$	DIP half saturation constant for DIA
${\mathrm{Ksil}}_{\mathrm{dia}}$	2.4	${\mathsf{\mu }\mathrm{molSi}\mathrm{L}}^{-1}$	SIL half saturation constant for DIA
${\mathrm{Kdin}}_{\mathrm{fla}}$	1.2	${\mathsf{\mu }\mathrm{molN}\mathrm{L}}^{-1}$	DIN half saturation constant for FLA
${\mathrm{Kdip}}_{\mathrm{fla}}$	0.1	${\mathsf{\mu }\mathrm{molP}\mathrm{L}}^{-1}$	DIP half saturation constant for FLA
${\mathrm{Rmine}}_{\mathrm{dia}}$	0.1	${\mathrm{day}}^{-1}$	Mineralization rate of DIA detritus
${\mathrm{Rmine}}_{\mathrm{fla}}$	0.15	${\mathrm{day}}^{-1}$	Mineralization rate of FLA detritus
N/P/Si	16/1/22	-	Atomic ratio of N, P, and Si

. The external sources of nutrients included river input, wastewater discharge, aquaculture, and atmospheric deposition. The exchange amounts of nutrients between the JZB and the YS depend on the differences in concentrations and the average residence time (52 days) (A19). It has been reported that the fluxes from sediments vary by nutrient, site, and redox conditions and directly observed benthic nutrient fluxes are lower than diffusive fluxes calculated from pore waters, indicating that the nutrients released from sediment pore waters are probably utilized in benthic production [27]. Therefore, the nutrients from benthic fluxes were not taken into consideration due to the limited observational data.

2.2. Model Configuration and Data Sources

2.2.1. CONTROL Run (the 2000s)

The simulation in the 2000s is named CONTROL. The external sources of nutrients were based on observations in the 2000s (

Table 2. The external sources of nutrients in the model.

	1960s	1980s	2000s	Unit	Reference
Dagu River runoff (F)	7.5	5.0	3.0	1011 L y−1	[38,41]
Riverine DIN concentration (${\mathrm{C}}_{\mathrm{river}-\mathrm{DIN}}$)	60	260	400	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[23,41]
Riverine DIP concentration (${\mathrm{C}}_{\mathrm{river}-\mathrm{DIP}}$)	3	5.3	7.0	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[23,41]
Riverine SIL concentration (${\mathrm{C}}_{\mathrm{river}-\mathrm{SIL}}$)	600	420	150	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[23,41]
DIN from wastewater ($\mathrm{L}{(\mathrm{DIN})}_{\mathrm{waste}}$)	0	80	200	106${\mathrm{mol}\mathrm{y}}^{-1}$	[23,26]
DIP from wastewater ($\mathrm{L}{(\mathrm{DIP})}_{\mathrm{waste}}$)	0	0.8	2.0	106${\mathrm{mol}\mathrm{y}}^{-1}$	[23,26]
SIL from wastewater ($\mathrm{L}{(\mathrm{SIL})}_{\mathrm{waste}}$)	0	3.2	8.0	106${\mathrm{mol}\mathrm{y}}^{-1}$	[23,26]
DIN from aquaculture ($\mathrm{L}{(\mathrm{DIN})}_{\mathrm{aqu}}$)	0	1.8	54	106${\mathrm{mol}\mathrm{y}}^{-1}$	[24,39]
DIP from aquaculture ($\mathrm{L}{(\mathrm{DIP})}_{\mathrm{aqu}}$)	0	0.12	3.8	106${\mathrm{mol}\mathrm{y}}^{-1}$	[24,40]
DIN from atmospheric deposition ($\mathrm{L}{(\mathrm{DIN})}_{\mathrm{atm}}$)	56	56	56	106${\mathrm{mol}\mathrm{y}}^{-1}$	[27]
DIP from atmospheric deposition ($\mathrm{L}{(\mathrm{DIP})}_{\mathrm{atm}}$)	0.42	0.42	0.42	106${\mathrm{mol}\mathrm{y}}^{-1}$	[27]
SIL from atmospheric deposition ($\mathrm{L}{(\mathrm{SIL})}_{\mathrm{atm}}$)	0.93	0.93	0.93	106${\mathrm{mol}\mathrm{y}}^{-1}$	[27]
DIN concentration of YS (${\mathrm{C}}_{\mathrm{YS}-\mathrm{DIN}}$)	2.2	5.0	6.0	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[18,40]
DIP concentration of YS (${\mathrm{C}}_{\mathrm{YS}-\mathrm{DIP}}$)	0.17	0.3	0.3	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[18,40]
SIL concentration of YS (${\mathrm{C}}_{\mathrm{YS}-\mathrm{SIL}}$)	6.0	5.5	4.0	${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$	[18,40]

). The amounts of nutrients from rivers depended on both the riverine concentrations of nutrients and the discharge of fresh water (A18). In this model, only Dagu River was taken into consideration, and the other small rivers were ignored. The river discharge (in ${\mathrm{L}\mathrm{s}}^{-1}$) of the Dagu River in the 2000s was derived from the model results of the Soil and Water Assessment Tool (SWAT) [38]. The river discharge reached its peak in summer (July to September), up to 9.0 × 10⁴ ${\mathrm{L}\mathrm{s}}^{-1}$, and remained low in other seasons, especially in winter (December to February), with an annual total discharge of 3.0 × 10¹¹ L. The riverine concentrations of DIN, DIP, and SIL were based on observations in 2002 [23]. The annual amount of nutrients from wastewater was obtained from investigations in 2002 [23]. The annual aquaculture input was estimated from aquaculture conditions in 2008 [39]. Atmospheric depositions of nutrients, including wet and dry deposition, were based on rainwater and aerosol samples collected at the coastal site of JZB from April 2004 to March 2005 [27]. Nutrients from these three external sources (wastewater, aquaculture, and atmosphere) were assumed to have little seasonal variations, and the annual amounts were averaged to each time step for model calculation. The concentrations of nutrients in the YS were set to 6.0, 0.3, and 4.0 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ for DIN, DIP, and SIL, respectively, according to observations [18,40].

The physical forcing includes the annual cycles of sea surface light intensity (I₀, in W m⁻²) and seawater temperature (Temp, in °C), expressed by the following cosine functions of time (t in Julian Day) [42]:

$${\mathrm{I}}_0=240-140\times \cos \left[2\mathsf{\pi }\times \frac{\mathrm{t}+10}{365}\right],$$

(1)

$$\mathrm{Temp}=13.5-11.0\times \cos \left[2\mathsf{\pi }\times \frac{\mathrm{t}-40}{365}\right].$$

(2)

The average light intensity (I) in JZB is expressed as I = 0.4 × I₀, considering the light attenuation due to pure water and suspended particulate materials in the water column, which is derived in Appendix A.

The initial conditions of state variables were based on observed data from the winters of the 2000s [17,18]. The initial concentrations of DIN, DIP, and SIL were set to 13.0, 0.8, and 6.0 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$, respectively. The initial biomasses of DIA and FLA were 0.72 and 0.08 ${\mathsf{\mu }\mathrm{g}\mathrm{L}}^{-3}$, respectively, in the form of chlorophyll-a (Chl-a) concentrations. After running the model for three years, the outputs were close to the climatological state of the model. The outputs of the fifth year were used in the following analysis.

The data used to validate the model results were the averaged values of several observations in the 2000s (Figure 1) [17,18,19,29].

2.2.2. NORI, NOWA, NOAQ, and NOAT Runs

In order to examine the different effects of the external nutrient sources, four numerical experiments were carried out and named NORI, NOWA, NOAQ, and NOAT, in which nutrients were removed from river input, wastewater discharge, aquaculture activities, and atmospheric deposition, respectively. Except for removing one external source of nutrients, the other conditions used in these experiments were the same as CONTROL.

2.2.3. The 1960s and the 1980s Runs

As shown in Table 2, the runoff of the Dagu River decreased from 7.5 × 10¹¹ ${\mathrm{L}\mathrm{y}}^{-1}$ in the 1960s to 5.0 × 10¹¹ ${\mathrm{L}\mathrm{y}}^{-1}$ in the 1980s [38,43]. The riverine concentration of SIL in the 1960s and the 1980s was supposed to be proportional to the sediment loads of Dagu River, which were approximately four and three times the amount of that in the 2000s, respectively [41,43]. The SIL concentration in the river decreased from 600 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ in the 1960s to 420 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ in the 1980s (Table 2). The DIN (DIP) concentration in the river increased from 60 (3.0) ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ in the 1960s to 260 (5.3) ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ in the 1980s [23,41]. It was reported by the Qingdao Statistical Yearbook in 2005 that the amount of industrial and domestic wastewater increased dramatically in the past 20 years. The amount of wastewater in the 2000s was nearly twice that in the 1980s [19]. The nutrient loads in the 1980s were estimated to be 40% of the values in the 2000s, proportional to the amount of wastewater (Table 2). In the 1960s, industry activities were rarely thriving; therefore, the nutrients from wastewater discharge can be neglected. The amounts of DIN and DIP from aquaculture activities in the 1980s were estimated to be 1.8 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ and 0.12 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ due to the variation in aquaculture areas [24,44], while nutrients from aquaculture were ignored in the 1960s. Nutrients from atmospheric deposition in the 1960s and 1980s were regarded to be the same as those in the 2000s due to the data limitation.

Although the human-induced changes in the coastline position-configuration and nearshore bathymetry have resulted in substantial changes to the residual current patterns, especially in northeastern JZB, the average residence time did not change apparently from 1966 to 2008 [45]. The same residence time of 52 days was used to calculate the exchange with YS, but the nutrient concentrations of YS changed a lot for different times (Table 2). The annual cycles of sea surface light intensity and seawater temperature were treated the same as the CONTROL Run.

3. Results

3.1. Parameter Sensitivity Analysis

Model sensitivity experiments were conducted to examine how the predicted variables vary with the changes of the model parameters within a given range (

Table 3. Sensitivities of DIA and FLA annual mean biomass to model parameters.

Parameter	Variation of Parameter (%)	Sensitivities of Annual DIA Biomass	Sensitivities of Annual FLA Biomass
${\mathrm{Pmax}}_{\mathrm{dia}}$	+50/−50	0.71/0.93	0.13/0.12
${\mathrm{Pmax}}_{\mathrm{fla}}$	+50/−50	0/0	0.99/0.99
${\mathrm{resp}}_{\mathrm{bas}}$	+50/−50	0.34/0.34	0.22/0.22
${\mathrm{resp}}_{\mathrm{pho}}$	+50/−50	0.06/0.06	0.06/0.06
${\mathrm{r}}_{\mathrm{mort}}$	+50/−50	0.62/1.66	0.65/1.91
${\mathrm{Iopt}}_{\mathrm{dia}}$	+50/−50	0.54/0.63	0.11/0.12
${\mathrm{Iopt}}_{\mathrm{fla}}$	+50/−50	0/0	0.31/0.17
${\mathrm{Kdin}}_{\mathrm{dia}}$	+50/−50	0/0	0/0
${\mathrm{Kdip}}_{\mathrm{dia}}$	+50/−50	0/0	0/0
${\mathrm{Ksil}}_{\mathrm{dia}}$	+50/−50	0.31/0.39	0/0.01
${\mathrm{Kdin}}_{\mathrm{fla}}$	+50/−50	0/0	0/0
${\mathrm{Kdip}}_{\mathrm{fla}}$	+50/−50	0/0	0.17/0.20
${\mathrm{Rmine}}_{\mathrm{dia}}$	+50/−50	0.05/0.14	0/0.02
${\mathrm{Rmine}}_{\mathrm{fla}}$	+50/−50	0/0	0/0

). The predicted state variables in these sensitivity experiments were the annual mean biomasses of DIA and FLA. Taking the parameters listed in Table 1 as the control run, the value of any selected parameter was increased and decreased by 50%. The model was run for three years in each sensitivity experiment, and the results of the third year were used. The sensitivity of a predicted state variable to the selected parameter was quantified as a factor $\mathrm{s}=\frac{\left|\mathsf{\Delta }\mathrm{A}/\mathrm{A}\right|}{\left|\mathsf{\Delta }\mathrm{B}/\mathrm{B}\right|}$, where A is the mean biomass of DIA or FLA and B is the value of one model parameter. ∆A is the variation of A when the parameter varies by ∆B.

The annual mean biomasses of DIA and FLA are sensitive to their own maximum growth rate and mortality rate (s > 0.5). The sensitivity factors about the basic respiration were 0.34 and 0.22 for DIA and FLA, respectively. The annual mean biomass of DIA was more sensitive to the half saturation constant of SIL (s = 0.31−0.39) than those of DIN (s = 0) and DIP (s = 0), while that of FLA was more sensitive to its half-saturation constant of DIP (s = 0.17−0.20) than that of DIN (s = 0). The sensitivity of the half saturation constant of DIN was low, indicating that the growth of DIA and FLA is not limited by DIN. The sensitivity analysis also quantified the competition between the two types of phytoplankton. The annual mean biomass of FLA had a sensitivity to some parameters related to the growth of DIA, such as the maximum growth rate of DIA (s = 0.12−0.13), the optimum light intensity for DIA (s = 0.11−0.12), but not vice versa. It is suggested that DIA is more competitive than FLA. On the other hand, the limiting nutrients for both DIA and FLA did not change when changing the model parameters with a large range (i.e., 50%). Evidently, the model-based conclusions on limiting nutrient and the phytoplankton community are robust.

3.2. Seasonal Variations of Nutrient Concentrations and Phytoplankton Biomasses

The simulated concentrations of the three types of nutrients in JZB exhibited apparent seasonal variations (Figure 3). All three nutrient elements accumulated in winter, when the growth of phytoplankton is limited by low water temperature and weak solar radiation. Then, the concentrations of DIN, DIP, and SIL decreased to about half of their peak values in spring along with the vigorous primary production. In summer, the concentrations increased gradually due to the large river input, and the concentrations of DIN and SIL increased more significantly than DIP. In autumn, the concentrations of DIP and SIL continuously increased in response to the reduction of phytoplankton photosynthesis. However, the concentrations of DIN decreased slightly due to the low river input.

The biomasses of both DIA and FLA experienced an apparent bloom in summer and had their lowest value in winter. The biomasses of DIA and FLA reached the peak values of 1.52 ${\mathsf{\mu }\mathrm{g}\mathrm{L}}^{-1}$ and 0.52 ${\mathsf{\mu }\mathrm{g}\mathrm{L}}^{-1}$ in August (Figure 3d), respectively, consistent with the maximum river runoff in summer. The biomass of DIA was higher than that of FLA throughout the year, especially in summer. The annual mean ratio of DIA to FLA was 2.3, with the maximum value of 2.9 in August and the minimum value of 1.6 in January. Therefore, DIA is the dominant phytoplankton species in JZB, which is consistent with the observations by [46], accounting for approximately 75% of the total phytoplankton in August.

It has been reported that there is a significant correlation between Chl-a and phytoplankton density when the dominant phytoplankton species are stable [47]. Since the biomasses of DIA and FLA represented by Chl-a concentrations could not be observed, the observed monthly density of DIA and FLA were used to validate the modeled seasonal variations of Chl-a concentrations. DIA density was obviously higher than that of FLA all throughout the year, and they both experienced similar seasonal variations with the peak in summer, which could also be captured by the modeled Chl-a concentrations of DIA and FLA [29]. Observations in JZB showed two phytoplankton blooms occurring in February and August [48,49]. Appearance of some cold-water species is suggested to be responsible for the February bloom [18]. Because the two types of phytoplankton in our model did not represent the cold-water species, we could not reproduce the bloom in February. Since the cold-water species of phytoplankton did not appear every February [19], the model was able to simulate the main seasonal variations of nutrient concentrations and phytoplankton biomasses in JZB.

3.3. Seasonal Variations of Nutrient Limitations

The nutrient limitations on the growth of phytoplankton are given by the minimum value of the three factors (${\mathrm{DIN}}_{\lim }=\mathrm{DIN}/(\mathrm{DIN}+{\mathrm{K}}_{\mathrm{DIN}}),{\mathrm{DIP}}_{\lim }=\mathrm{DIP}/(\mathrm{DIP}+{\mathrm{K}}_{\mathrm{DIP}}),{\mathrm{SIL}}_{\lim }=\mathrm{SIL}/(\mathrm{SIL}+{\mathrm{K}}_{\mathrm{SIL}})$), where ${\mathrm{K}}_{\mathrm{DIN}},{\mathrm{K}}_{\mathrm{DIP}},{\mathrm{K}}_{\mathrm{SIL}}$ are half-saturation constants (in Appendix A). The limiting nutrient in JZB was determined through comparing the value of three factors. The simulated annual mean concentrations of the DIN, DIP, and SIL were 19.8, 0.56, and 5.73 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$), respectively. The ratio of DIN/DIP was 38/1, much larger than the Redfield ratio. Both the high concentration of DIN and the ratio of DIN/DIP suggest that DIN was in surplus compared with DIP. As shown in Figure 4, the values of ${\mathrm{DIN}}_{\lim }$ for DIA and FLA were above 0.9 during the whole year, suggesting that DIN never served as the limiting nutrient for the growth of phytoplankton. The limiting nutrient for DIA was SIL all year around (Figure 4). The SIL_lim decreased from January to June, and the limitation effect was strongest in June, when its value was 0.56. Then, the limitation effect gradually weakened due to the riverine SIL input in summer. The growth of FLA was limited by DIP (Figure 4) all year round. The mean DIP_lim for it was about 0.82, which was always higher than the SIL_lim for DIA. Given that DIA is the dominant species in JZB, SIL_lim appears favorable for controlling the ecological equilibrium of JZB.

3.4. Annual Budgets of Nutrients

Figure 5 shows the annual budgets of DIN, DIP, and SIL. The four external sources of nutrients (river, wastewater, aquaculture, and atmosphere) were treated as forcing conditions in the model. Their amounts are listed in Table 2. The budgets of biogeochemical processes represent the renewal and cycle rate of nutrients in JZB. Phytoplankton plays an important role in the renewal of nutrients. The net production of phytoplankton (photosynthesis minus respiration) is the main sink, with 316.89 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ of DIN, 19.81 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ of DIP and 382.37 × 10⁶${\mathrm{mol}\mathrm{y}}^{-1}$ of SIL. The maximum net production of DIA and FLA occurred in summer, consistent with their high biomasses. The mineralization of detritus released 317.19 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ of DIN, 19.83 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ of DIP, and 382.76 × 10⁶ ${\mathrm{mol}\mathrm{y}}^{-1}$ of SIL, which are approximately 70% of the amount of photosynthesis absorption. According to the model results, the annual mean concentrations of the three nutrients in JZB are higher than those in YS. Annually, 432.09 × 10⁶ mol of DIN, 8.36 × 10⁶ mol of DIP, and 54.70 × 10⁶ mol of SIL are transported out from JZB to YS through water exchange. The nutrient budgets are in balance in the whole year, which indicates that the model followed the mass conservation.

3.5. Estimation of the Contributions of Nutrient External Sources on Nutrient Limitations and Phytoplankton Community

The responses of the nutrients and phytoplankton structures in JZB to each external nutrient source were examined by comparing the results of NORI, NOWA, NOAQ, and NOAT with those of CONTROL (Figure 6 and Figure 7). The nutrients flux of atmospheric deposition was the smallest among the four external sources. The DIA and FLA biomasses decreased slightly by 0.2% and 0.5% in NOAT, respectively. Removing atmospheric deposition had a negligible influence on the model results of both phytoplankton biomasses and nutrient concentrations. The results of the other three experiments are discussed in detail bellow.

3.5.1. The Effects of River Input

The Dagu River inputs 45.79 × 10⁶ mol SIL into JZB per year, contributing up to 80% of the annual amount from the four external sources. Meanwhile, it is a large source of DIN (122.10 × 10⁶ mol y⁻¹), just less than that from wastewater discharge. The amount of DIP from river input is 2.14 × 10⁶ mol y⁻¹, which is similar as that from wastewater but much less than that from aquaculture.

Comparing the results of NORI with those of CONTROL, the annual mean concentrations of DIN, DIP, and SIL decreased by 19%, 11%, and 25%, respectively. The decreases are apparent in all these three nutrient elements (Figure 6). Moreover, the nutrient supply from the river exhibited an obvious seasonal variation due to the significantly different runoff during the wet and dry seasons. As a result, nutrient concentrations hardly changed in the dry seasons and changed a lot in the wet seasons. In summer, the peak of DIN disappeared in NORI, and the concentration of DIN increased from July, and then began to decline after reaching the peak in November (Figure 6a). The concentrations of DIP and SIL decreased continuously in summer in NORI rather than increase in CONTROL (Figure 6b,c). After reaching the trough in July, the DIP and SIL concentrations in NORI increased from August, and then began to decrease after reaching the peak in November. The peak time was one month later than other numerical experiments. The maximum decrease of SIL and DIP also occurred in summer, when the values were only half those in CONTROL.

The values of SIL and DIP limiting factors in NORI decreased evidently in summer, enhancing the SIL limitation for DIA and DIP limitation for FLA (Figure 7). The SIL limiting factor in spring and winter changed relatively slightly. SIL still acted as the limiting nutrient for the growth of DIA all year round, but SIL limitation was strengthened in NORI (Figure 7).

The annual biomass of DIA decreased by 17%, while the decreasing percentage of FLA was only 3%. The biomasses of DIA and FLA changed the most in summer, which were reduced by 30% and 6%, respectively. The annual mean ratio of DIA to FLA decreased from 2.3 to 2.0, indicating that the river input of nutrients has stronger influences on DIA than on FLA. The enhanced SIL limitation leads to a phytoplankton structure favored to FLA.

3.5.2. The Effects of Aquaculture Activities

Aquaculture provides 54 × 10⁶ mol y⁻¹ DIN to JZB, which is similar to the amount from atmospheric deposition, but only a quarter of that from wastewater discharge. However, it is the main source of DIP, with an annual amount nearly twice that from wastewater or river. The ratio of DIN to DIP from aquaculture is 14, indicating that aquaculture is the only source with a lower ratio of DIN/DIP than the Redfield ratio.

DIN concentration in NOAQ exhibited similar seasonal variations as that in CONTROL, but the value was about 2.0 ${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$ less than that in CONTROL (Figure 6a). Among the results of CONTROL and four experiments, the DIP concentration in NOAQ almost maintained the lowest throughout the year (Figure 6b), only slightly larger than that in NORI in summer. No SIL comes from aquaculture, so the SIL concentration did not change in NOAQ.

The DIP limiting factor for FLA decreased by 0.05 compared with that in CONTROL (Figure 7). The FLA biomass decreased by 5% annually and changed the most in June, which was decreased by 8%, while the DIA biomass had nearly no change all year round without this source. The ratio of DIA to FLA increased from 2.3 in CONTROL to 2.4 in NOAQ. Therefore, FLA is more sensitive to this source than DIA, and aquaculture activities have effects on the miniaturization of phytoplankton size.

3.5.3. The Effects of Wastewater Discharge

Wastewater is the largest source of DIN, with an annual amount of 200 × 10⁶ mol y⁻¹, which is nearly twice that from rivers. However, it only contributes half of the amount of DIP from aquaculture and one sixth of the amount of SIL from rivers. The ratio of DIN/DIP/SIL is approximately 100/1/4 in wastewater (Figure 5), which suggests that wastewater brings more DIN into JZB compared with DIP and SIL.

The seasonal variations of nutrient concentrations, nutrient limiting factors, and phytoplankton biomasses in NOWA were similar to those in CONTROL (Figure 6). However, the values were different, especially for DIN concentration. DIN concentration in NOWA remained the lowest among the four experiments nearly all year round, and only larger than that in NORI in summer (Figure 6a). DIP concentration in NOWA was also lower than that in CONTROL all through the year, however there was even no considerable change in SIL concentration. Wastewater discharge brings more excessive DIN into JZB, therefore, the limiting factors of SIL for DIA and DIP for FLA in NOWA became smaller than those in CONTROL, indicating a strengthened SIL limitation for DIA and DIP limitation for FLA.

FLA biomass decreased by 3%, while DIA biomass decreased by 2% annually in NOWA. Seasonally, the DIA and FLA biomasses decreased more in winter and spring than in summer and autumn, because the river input relieved the unbalance of nutrient structures caused by the wastewater discharge in summer. In general, the ratio of DIA to FLA increased slightly annually.

3.5.4. The Effects of External Sources on the Cycling Rate of Nutrients in JZB and the Exchange with the YS

Besides the changes in the nutrient concentrations and phytoplankton biomasses in JZB when removing each external source of nutrients, the cycling rates of nutrients among different state variables exhibited even more apparent variations, which could not be revealed directly by in-site observations. Comparing the annual budgets of SIL in NORI and CONTROL (

Table 4. Simulated nutrient budgets (Unit: × 10⁶ mol) during the whole year in CONTROL, NORI, NOWA, NOAQ, and NOAT.

		DIA Photosynthesis	DIA Respiration	FLA Photosynthesis	FLA Respiration	DSi Mineralization	DNSi Mineralization	Exchange with YS
CONTROL	DIN	−395.83	+117.75	−49.91	+11.10	+278.37	+38.82	−432.09
	DIP	−24.74	+7.36	−3.12	+0.69	+17.40	+2.43	−8.36
	SIL	−544.27	+161.90	0	0	+382.76	0	−54.70
NORI	DIN	−266.29	+91.83	−46.50	+10.61	+174.62	+35.90	−309.99
	DIP	−16.64	+5.74	−2.90	+0.66	+10.92	+2.24	−6.22
	SIL	−366.15	+126.25	0	0	+240.10	0	−8.92
NOWA	DIN	−384.07	+115.40	−47.81	+10.79	+268.96	+37.03	−232.09
	DIP	−24.00	+7.21	−2.98	+0.67	+16.81	+2.31	−6.36
	SIL	−528.10	+158.69	0	0	+369.82	0	−46.70
NOAQ	DIN	−395.83	+117.75	−45.40	+10.44	+278.37	+35.00	−378.09
	DIP	−24.74	+7.36	−2.84	+0.65	+17.40	+2.19	−4.56
	SIL	−544.27	+161.90	0	0	+382.76	0	−54.70
NOAT	DIN	−394.49	+117.48	−49.49	+11.04	+277.29	+38.47	−376.09
	DIP	−24.66	+7.34	−3.09	+0.69	+17.33	+2.40	−7.94
	SIL	−542.42	+161.54	0	0	+381.29	0	−53.77

), the annual photosynthesis of DIA was 366.15 × 10⁶ mol SIL when removing river input, which is about 33% less than that in CONTROL. Meanwhile, the respiration of DIA and the mineralization of DSi decreased by about 37%. The decrease in each biochemical process of SIL cycle is even more apparent than that in the standing stocks of SIL in JZB when removing river input. Similar situations occurred in the DIP cycle when removing nutrients from aquaculture. The annual amount of DIP related process (the photosynthesis and respiration of FLA, the mineralization of DNSi) decreased up to 10% when removing aquaculture, while the decrease in annual FLA biomass was 5%. These variations in NOWA are not as much as those in NORI and NOAQ. That is because wastewater discharges mainly DIN to JZB, which is not the limiting nutrient. Therefore, it can be concluded that removing external sources of limiting nutrients obviously slows down the cycling rate of nutrients in JZB.

In addition, there were obvious decreases in the exchanging amounts of nutrients with the YS in these four numerical experiments. When removing river input, the transporting amount of SIL from JZB to the YS decreased from 54.70 × 10⁶ mol y⁻¹ to 8.92 × 10⁶ mol y⁻¹. Wastewater discharge has less influence on the phytoplankton biomasses than river input and aquaculture activities. However, the exchanging amount of DIN in NOWA decreased to only half of that in CONTROL. So, besides the effects on JZB itself, the external sources of nutrients have more remote influences on the adjacent seas through water exchange controlled by complicated hydrodynamic processes.

3.6. Long-Term Variations of the Nutrient Limitation and Phytoplankton Community

As shown in Section 3.5, the external sources of nutrients have various influences on the nutrient structure and phytoplankton community in the JZB. Since the nutrient loads from the external sources have dramatically changed by human activities in the past several decades, it may result in changes in nutrient limitation and then the phytoplankton community in JZB. Considering different amounts of nutrients from the four external sources in different decades (Table 2), the model reproduced the nutrient structure and phytoplankton community in the 1960s and the 1980s.

The annual cycles of nutrient concentrations and phytoplankton biomasses in these two decades showed similar seasonal patterns to the results in the 2000s. The simulated concentrations of DIN and DIP in the 1960s were in good agreement with the observations (

Table 5. The simulated and observed nutrient concentrations (${\mathsf{\mu }\mathrm{mol}\mathrm{L}}^{-1}$) in the 1960s, 1980s, and 2000s [17,18,31].

	1960s		1980s		2000s
	Simulation	Observation	Simulation	Observation	Simulation	Observation
DIN	0.98–5.43	1.00–4.00	12.05–23.65	2.55–19.50	17.24–37.74	3.60–44.21
DIP	0.18–0.46	0.10–0.22	0.08–0.52	0.20–0.72	0.09–1.37	0.09–2.00
SIL	8.38–62.83	-	5.59–27.66	0.20–11.50	2.44–12.51	1.10–20.39

). The concentration of SIL in the 1960s could not be evaluated due to there being too few observations. The simulated concentrations of DIN, DIP, and SIL in the 1980s almost fell into the ranges of observations in this decade (Table 5). The concentrations of DIN and DIP in the 2000s were about six times and triple those in the 1960s, respectively. However, the SIL concentration in the 2000s was only half of that in the 1980s and one quarter of that in the 1960s.

The increase in DIN was mainly caused by the river input and wastewater discharge. The river input was determined by both the runoff and the riverine concentrations. Although the annual river runoff decreased from 7.5 × 10¹¹ L in the 1960s to 3.0 × 10¹¹ L in the 2000s, the annual amount of DIN from the river still increased from 45 × 10⁶ mol in the 1960s to 120 × 10⁶ mol in the 2000s due to the dramatic increase of the riverine DIN concentration. The continuously prosperous industry brought about more and more wastewater, which made wastewater become the biggest source of DIN until the 2000s (Table 2). Consequently, the annual DIN input from the four external sources increased from 101 × 10⁶ mol in the 1960s to 268 × 10⁶ mol in the 1980s, and then reached 430 × 10⁶ mol in the 2000s. The expansion of aquaculture activities was mainly responsible for the rising input of DIP to JZB. The total DIP input increased from 2.67 × 10⁶ mol in the 1960s to 3.99 × 10⁶ mol in the 1980s, then to 8.32 × 10⁶ mol in the 2000s. Overall, the DIN concentration in JZB experienced a faster increase than that of DIP. On the contrary, the concentration of SIL showed an apparently decreasing trend from the 1960s to the 2000s (Table 5), which was definitely caused by the decreases in both the river runoff and the riverine SIL concentration.

The long-term variations in nutrient concentrations changed the nutrient limitation situations and the phytoplankton compositions in JZB. In the 1960s, the limiting nutrient for both DIA and FLA was DIP (Figure 4). The DIP limiting factor for DIA was smaller than that for FLA due to the larger DIP half-saturation constant for DIA, indicating the stronger limitation for DIA. However, the dominant species was still DIA, due to its larger maximum growth rate. The annual mean ratio of DIA to FLA was approximately 2.7. The nutrient limitation in the 1980s changed dramatically compared with the situation in the 1960s. DIP was still the limiting nutrient for FLA. However, DIP and SIL alternatively limited the growth of DIA in the 1980s (Figure 4a). From January to June, the limiting factor of SIL (${\mathrm{SIL}}_{\lim }$) was smaller than that of DIP ($DI{P}_{lim}$), indicating the concentration of SIL controlling the growth of DIA. From July to November (the flood season), the limiting nutrient shifted from SIL to DIP. Although the concentration of SIL decreased from the 1960s to the 1980s, the increasing wastewater discharge and aquaculture activities brought a large amount of DIP and relieved the DIP limitation, making nutrient limitation in the 1980s weaker than the situation in the 1960s. As a result, the annual mean biomass of DIA increased from the 1960s (0.64 ${\mathsf{\mu }\mathrm{g}\mathrm{L}}^{-1}$) to the 1980s (0.80 ${\mathsf{\mu }\mathrm{g}\mathrm{L}}^{-1}$). The annual ratio of DIA to FLA in the 1980s was approximately 2.9, which was larger than that in the 1960s. The ecosystem in JZB changed from DIP limitation to SIL limitation under the influence of reduced SIL input, and relatively increased DIP input from the 1980s to the 2000s. The enhanced SIL limitation from the 1980s to the 2000s (Figure 4a) led to a more apparent increase in FLA biomass and a comparatively insignificant increase of DIA biomass. The annual ratio of DIA to FLA greatly reduced from 2.9 to 2.3, which resulted in the miniaturization of the phytoplankton community.

4. Implications and Limitations

The nutrient limitation may be caused by low concentrations or alternations due to changing ratios. The more obvious changes in phytoplankton are not in the biomass, but in the community composition in JZB. This shift in phytoplankton species was mainly caused by the change in nutrient ratios. From the 1980s to the 2000s, the DIN/DIP/SIL changed from 31/1/26 to 35/1/10, according to the model results. The increases in DIN and DIP concentrations (especially DIN) and the decrease in SIL concentration resulted in a more unbalanced nutrient structure [50]. DIP limitation for DIA growth has almost been eliminated. Meanwhile, strengthened SIL limitation has led to the decrease of DIA, and the smaller species have taken over the niche vacated by the larger species [20].

River input is the main source of SIL, and therefore becomes the only source to support a diatom dominant phytoplankton community in JZB. Damming usually decreases SIL load, which makes SIL a potential limiting nutrient in coastal waters [51,52], such as Guadiana estuaries in Mediterranean [53], Sweden, Finland [54], and the Baltic Sea [14]. The sharp decrease of SIL flux caused higher frequency of occurrence of harmful algal bloom [55,56]. In addition, low SIL concentration limited the silicification ability of DIA, and then reduced their absorption of Cd in the seawater, which finally resulted in heavy metal pollution [57] and toxin enrichment through the aquatic food chain [58]. In the North Sea, more reduction in riverine P loads than N loads caused large offshore gradient of DIN/DIP ratio in spring, from 375:1 nearshore toward 1:1 in the central North Sea and the phytoplankton growth was limited by P and Si nearshore, co-limited by N and P in a transitional region, and limited by N in the outer-shore waters [59].

Nutrients from aquaculture, with a lower DIN/DIP ratio compared with those from other nutrient sources, led to the more obvious growth of FLA than that of DIA in JZB. It was reported that the phytoplankton community experienced a shift from DIA-dominated to FLA-dominated based on observations in an offshore shellfish-cultured area located in the northern Yellow Sea of China [60]. There was also a worldwide phenomenon that excessive discharge of DIP increased the risks of eutrophication and threatened the environmental regulations in the aquaculture waters [61]. So, the utilization efficiency of nutrients was regarded to be an important indicator, whether the aquaculture activities were sustainable or not [60]. In addition, it has been reported that aluminum pillared bentonite can decrease the P loads from aquaculture, and it can be used for the treatment of aquaculture discharge water before releasing into the coastal seas [62].

The effects of wastewater vary due to its various components in different countries. A previous study [63] reported that the low N/P ratio of wastewater inputs resulted in a DIN-limited and DIP-saturated system in the estuaries of Chesapeake Bay, which is totally different from the situation in JZB. Reduced loads of nutrients from the Luggage Point sewage treatment plants into the Moreton Bay resulted in a decrease in DIN concentration and phytoplankton biomass [64]. In the Scheldt estuary of the Netherlands, upgrading the wastewater treatment plant dramatically reduced DIN and DIP delivery to the coastal seas. In the absence of decreases in SIL export, the reduction of DIN and DIP fluxes and the improvement of the DIN/DIP ratio compared to SIL favored the growth of DIA [65]. After decreases in N and P inputs through a wastewater-improvement project, the trophic status shifted from hyper-eutrophic to eutrophic-mesotrophic in Boston Harbor. As a result, the N limitation relative to P limitation increased and the primary production and phytoplankton biomass decreased [66]. Nitrogen reduction reduced the winter spring diatom bloom and winter Chl-a levels were lower significantly after nutrient reduction in Narragansett Bay [67].

The model results showed that the atmospheric deposition had little influence on the variations in phytoplankton biomass and composition in JZB. The atmospheric deposition provided a similar amount of DIN to JZB compared with the aquaculture. However, it brought less DIP and SIL, which acted as the limitation nutrients in JZB in the 2000s. It was widely reported that atmospheric deposition was an important source of nutrients and promoted primary production in open oceans [68,69]. The atmospheric deposition of nitrogen contributed about 10% of the primary production in the western Pacific Ocean [70]. Even in regional seas, the atmospheric deposition is also an important source of nutrients. Taking YS for example, the estimated annual amount of nitrogen from atmospheric deposition was about five times that from the rivers. However, their contributions to the annual primary production of the whole Yellow Sea were not very high, because there was already surplus DIN in the seawaters. Only in certain areas and certain periods can their contributions not be neglected. The nitrogen from atmospheric deposition promoted the spring phytoplankton bloom in the central part of the Yellow Sea when DIN was used up by the fast growth of phytoplankton [71]. During heavy-polluting weathers, such as dust or haze events, the amount of nitrogen from atmospheric deposition increased to twice as much as that in normal days [72], which aggravated the DIP and SIL limitation situation and then affected the composition of the phytoplankton community [73].

It was also indicated by the nutrient budgets that the external sources of nutrients have subsequent effects on the adjacent seas besides the bay itself. When considering the effects of climate change, the remote influences of the oceans on the marginal seas always gain more attraction. However, it cannot be neglected that anthropogenic activities can also affect the remote oceans through coastal waters and marginal seas. For example, the material exchange between the East China Sea shelf (marginal sea) and Kuroshio affects not only the ecosystem of the East China Sea but also that of the Northwestern Pacific Ocean [37,74,75].

It should be noted that the influence discussed above primarily focuses on the entire JZB. The box ecosystem model was effective when focusing on the temporal variation of nutrient composition and succession of the phytoplankton community of the entire JZB under the influence of anthropogenic activities. However, it could not differentiate the spatial variations in JZB due to inhomogeneity of the nutrient sources, e.g., the aquaculture activities and wastewater have more significant impacts on shore areas. In addition, the nutrients from benthic fluxes were not taken into consideration due to the limited observational data. This is, however, probably important to the shore area because of presence of large tidal flat in JZB. Therefore, a fully three-dimensional model with a moving boundary at shore area is probably necessary in the future to further understand the spatial variation in the bay.

5. Conclusions

The simulated annual cycles of nutrient concentrations and phytoplankton biomasses in the 2000s suggest that the growth of DIA is limited by SIL, while that of FLA is limited by DIP. Changes in the external nutrient sources lead to long-term variations (the 1960s, the 1980s, and the 2000s) in nutrient limitation and the phytoplankton community. Removing the nutrients from the river input enhances SIL limitation for DIA and decreases the ratio of DIA to FLA from 2.3 to 2.0. Aquaculture activities are the main source of DIP; therefore, the FLA biomass decreases when removing this external source. The wastewater has a relatively weaker impact on the change of phytoplankton communities than aquaculture, and the influence of atmospheric deposition is likely negligible.