Constraining the response of phytoplankton to zooplankton grazing and photo-acclimation in a temperate shelf sea with a 1-D model towards S2P3 v8.0

. An established 1-dimensional model of Shelf Sea Physics and Primary Production (S2P3) has been developed into three different new models: S2P3-NPZ which includes a Nutrient-Phytoplankton-Zooplankton (NPZ) framework, where the grazing rate is no longer ﬁxed, but instead varies over time depending on different functions chosen to represent the predator-prey relationship between zooplankton and phytoplankton; S2P3-Photoacclim which includes a representation of the process of photo-acclimation and ﬂexible stoichiometry in phytoplankton; and S2P3 v8.0 which combines the NPZ framework and 5 the variable stoichiometry of phytoplankton at the same time. These model formulations are compared to buoy and CTD observations, as well as zooplankton biomass and in situ phytoplankton physiological parameters obtained in the Central Celtic Sea (CCS). Models were calibrated by comparison to observations of the timing and magnitude of the spring phytoplankton bloom, magnitude of the spring zooplankton bloom, and phytoplankton physiological parameters obtained throughout the water column. A sensitivity study was also performed for each model to understand the effects of individual parameters on 10 model dynamics. Results demonstrate that better agreement with biological observations can be obtained through the addition of representations of photo-acclimation, ﬂexible stoichiometry, and grazing provided these can be adequately constrained.


Introduction
Shelf seas are ocean regions where water depth is less than a few hundred metres (∼ 200 m) and represent only ∼ 10% by area of the global ocean. However, these systems have a disproportionate importance because of their exceptionally high biological 15 productivity (Holt and Proctor, 2008), being responsible for 15 to 30% of the total oceanic primary production (PP) (Wollast, 1998;Muller-Karger et al., 2005;Davis et al., 2014). Research vessels and remote, autonomous vehicles have been used to study shelf sea regions such as in the SSB programme (https://www.uk-ssb.org/), whose aim was to increase the understanding of how physical, chemical and biological processes interact on UK and European shelf seas, collecting observations throughout 2014 and 2015 in different regions of the UK shelf sea, although this data is not synoptic (i.e. it is not sampled at different 20 locations simultaneously). To complement the available data from research vessels, ocean models have been used to study and understand marine biogeochemistry, including a variety of high spatial resolution models to represent the biogeochemistry of shelf seas with high complexity and horizontal spatial resolution (Sharples, 1999(Sharples, , 2008Edwards et al., 2012;Marsh et al., In contrast, simpler models like the Shelf Sea Physics and Primary Production (S2P3) model (Sharples, 1999;Simpson and Sharples, 2012) have been used to study the dynamics of shelf seas and to simulate seasonal stratification with greater computational efficiency by using a 1-D Nutrient-Phytoplankton (NP) model to represent physical and biogeochemical processes in the water column. In temperate shelf seas, away from advective sources such as the shelf break or plumes from rivers, horizontal processes can be neglected in comparison to vertical processes, thus, although S2P3 does not consider advective fluxes it can 40 make a good representation of the dynamics in the water column for temperate shelf seas Sharples, 2008;Marsh et al., 2015).
The S2P3 v7.0 model was introduced and developed in the work of Marsh et al. (2015), where it was outlined that further development of that model would include resolving phytoplankton physiology. In this study, the S2P3 v7.0 model is devel-45 oped by allowing variations in light intensity to produce phenotypic adjustments in the phytoplankton cells by changing the chlorophyll content of the phytoplankton and, therefore, the cellular absorption cross-section (Macintyre et al., 2002). This phenotypic change in response to variations in the photon flux density is called photo-acclimation (Falkowski and Laroche, 1991;Moore et al., 2006). The main property of photo-acclimation is the reduction of photosynthetic pigment content in response to increased irradiance (Falkowski and Laroche, 1991). Moreover, changes in nutrient availability can further alter cellular 50 chlorophyll and nitrogen quotas (Droop, 1983;Geider et al., 1998), incorporating a combined representation of these two processes (Geider et al., 1998), this new version of S2P3 v7.0 we term S2P3-Photoacclim and relates phytoplankton growth rates to cell quota (Droop, 1983), describing the light-, nutrient-. and temperature-dependencies of phytoplankton growth rate to varying ratios of N : C : Chl (Geider et al., 1998). On the other hand, simpler models using an NPZ or NPZD framework, with the use of nutrients, phytoplankton, zooplankton, and detritus as the main model components (e.g., Steele, 1974;Wroblewski 55 et al., 1988;Anderson, 2005) have shown good agreement with observations in terms of chlorophyll and PP, by simulating the timing and magnitude of the spring phytoplankton bloom in different regions of the ocean. Despite their relative simplicity, NPZ models can be a better option to approach an understanding of the physics and biology of an ecosystem, which lead to a further development of the S2P3 v7.0 model where the simplest assumption of a fixed proportion of phytoplankton being grazed and remineralised into the DIN pool (grazing rate) is developed into an NPZ framework (S2P3-NPZ model), using a 60 Holling Type 2 or Ivlev grazing functional response of zooplankton grazing on phytoplankton (Franks, 2002), which shows a saturating response to increasing food.
A combination of photo-acclimation, flexible stoichiometry, and NPZ framework is then performed to produce a newly developed model called S2P3 v8.0. This paper presents a thorough analysis in terms of sensitivity to biological parameter values in 65 each new developed model, resulting in differences of the model structure. A comparison between each model demonstrates how structural differences influence the representation of the spring phytoplankton bloom and annual PP. The aim of this paper is to provide a better understanding about the predator-prey relationship between zooplankton and phytoplankton, the effects of photo-acclimation and flexible stoichiometry and thus the importance of these processes for accurate representation of phytoplankton dynamics in shelf sea modelling.

2 Study region and model setup
This study is focused on the Central Celtic Sea (CCS), a region located in the North-Western (NW) European Shelf, which is characterised by its tidally dynamic environment and summer stratification (Pingree et al., 1978;Sharples and Holligan, 2006;Hickman et al., 2012). Daily meteorological data, provided by the National Centers for Environmental Predictions (NCEP) 75 Reanalysis data (http://www.esrl.noaa.gov/psd) was used to force the model at the CCS site, located at 49.4 • N , 8.6 • W . Wind speed (m s −1 ), cloud coverage (%), air temperature ( • C), and relative humidity (%) variables from this dataset were all used to force each model version.
The following description of model setup is applied to each model structure developed in this work. Tidal components consist 80 of the u-component (semi-major axis) and the v-component (semi-minor axis) for the M 2 , S 2 , and N 2 tidal constituents. This data was obtained from a fine mesh (12km resolution) covering the UK shelf. Tidal currents are predicted using the Proudman Oceanographic Laboratory Coastal Ocean Modelling Systems (POLCOMS) 3-D shelf model Wakelin et al., 2009) with an output extracted for the CCS location. Moreover, each model is initialised on 1 st January of the first year of simulation with a temperature of 10.10 • C at all depths, and water column presumed mixed throughout, including a vertical 85 resolution of 1m (i.e. 140 vertical levels). Initial values of physical variables are consistent with former studies (Sharples, 1999(Sharples, , 2008Marsh et al., 2015), whereas initial values of biological variables are based on observations of zooplankton biomass at the CCS location over winter months of 0.02 mmol N m −3 (Giering et al., 2018); phytoplankton chlorophyll correspond to a typical winter value for the CCS location of 0.2 mg Chl m −3 , and the DIN initial value is 7 mmol N m −3 . The initialised variables are only set up at the start of each simulation and do not reset in between years.

90
The S2P3 v7.0 model can be divided into two different components: a physical part and a biological part. For this research, the model is an improved version of the original described in  to be compiled and executed in a Unix environment (Marsh et al., 2015), allowing a 1-D representation of physical and biological processes in shelf seas by simulating 95 the seasonal cycle of phytoplankton, water column stratification, and PP at a selected location defined by water depth and tidal current amplitude. The physical part of the model has been greatly described in many other studies (Sharples, 1999;Sharples, 2008;Simpson and Sharples, 2012;Marsh et al., 2015) and is not described in this section again. This model uses the turbulence closure scheme based on Canuto et al. (2001). Likewise, the biological part of the S2P3 v7.0 model is described in Marsh et al. (2015) (Figure 1a).

100
In order to explicitly account for the influence of zooplankton grazing and, hence, predator-prey dynamics, the S2P3 v7.0 model was developed into an NPZ framework. This new version of the model (S2P3-NPZ) includes zooplankton as a state variable, contrary to the S2P3 v7.0 model where grazing was calculated as a fixed seasonal cycle which was represented as a sink term in the phytoplankton tendency equation. Addition of an explicit grazer also allows comparison to zooplankton  The S2P3-Photoacclim model is a new version of the S2P3 v7.0 model, incorporating a more complete representation of phytoplankton physiology (Geider et al., 1998). This model allows phytoplankton to acclimate to changes in light and nutrients, therefore, the ratios of N : C : Chl and characteristics of phytoplankton physiology can vary, allowing direct comparison to physiological data. The biological part of the S2P3-Photoacclim model uses three currencies of phytoplankton biomass: carbon (P hyto C ; mg C m −3 ), nitrogen (P hyto N ; mg N m −3 ), and chlorophyll (P hyto chl ; mg Chl m −3 ). The S2P3-Photoacclim 115 model calculates phytoplankton growth as a function of both nitrogen assimilation and carbon fixation (i.e. variable Chl : N and Chl : C ratios). It is assumed that respiration (R) is equal for all cellular components as a function of temperature: synthesis (ρ chl ), which reflects the ratio of energy assimilated to energy absorbed (Geider et al., 1996), and with phytoplankton carbon regulated through a cost associated with biosynthesis (uζ) (Vries et al., 1974;Geider, 1992;Geider et al., 1998), respiration, and grazing (G); (d) S2P3 v8.0 model, including two different varying quotas for P hytoN : P hytoc and P hytoN : P hyto chl .

Validation of the models: observations
To calibrate or tune each model, they were adjusted on a trial-and-error basis until disagreement with the in situ observations was minimised, allowing investigation of the sensitivity of each model to changes in the parameters listed in Table 1. from long-term mooring deployments including Smartbuoy (Mills et al., 2003) were collected. Moreover, CTD (stainless and titanium) casts were performed in different locations of the NW European shelf from the CCS location to the shelf break, with discrete samples of temperature, DIN, and chlorophyll-a. At the CCS location, the CTD samples were collected from pre-dawn to midday with a 1m vertical resolution over the whole water column (140m depth). the S2P3 v8.0 models, only mesozooplankton biomass were considered, with a community composition that included: amphipods, appendicularian, chaetogratha, copepods, euphausiacea, polychaeta, and others (e.g. cladocerans, dinoflagellates, echinoderm, eggs, foraminifera, gymnosomata, unidentified larvae, nauplii, ostracods and radiolarian, all of which contributed < 3% in all samples). The complete data set can be obtained from the British Oceanographic Data Centre (BODC), (http: //www.bodc.ac.uk/data). Photosynthesis versus irradiance (P vs E) experiments were conducted in short-term incubations (2 -4h) using a photosyn-155 thetron (Moore et al., 2006). From these P vs E experiments chlorophyll-a normalised PP was derived from 14 C uptake to obtain the chlorophyll-a specific maximum light-saturated photosynthesis rate P Chl max (mg C (mg Chl − a) −1 h −1 ) and the maximum light utilisation coefficient, (Jassby and Platt, 1976;Hickman et al., 2012).
Multiple samples were collected in the surface and in deeper layers to obtain different phytoplankton populations throughout the photoperiod. Values of α chl and the light saturation parameter, E k (µE m −2 s −1 ) (given by E k = P Chl max /α chl ) were spec-160 trally corrected to the in situ irradiance at the sample depth according to the phytoplankton light absorption (Moore et al., 2006). For this work, observations from both cruises were used considering only the stations from the seasonally stratified sites 6 https://doi.org/10.5194/gmd-2019-345 Preprint. Discussion started: 30 January 2020 c Author(s) 2020. CC BY 4.0 License.

165
(µE m −2 s −1 ) −1 ; std = 4.38 × 10 −6 ), and therefore, the values of P Chl max and E k were used as variables for comparison with equivalent modelled values.
Loss rate of zooplankton due to predation and physiological death 0.05 0.02 Maximum value of the carbon-specific rate of photosynthesis 2.0 3.5 Maximum value of the chlorophyll : phytoplankton nitrogen ratio 0.3 0.15 Respiration rates 0 0.02

Model calibration
For this study approximately two years of phytoplankton chlorophyll data were available for the CCS location, while in the Moreover, the S2P3-Photoacclim and the S2P3 v8.0 models were also evaluated with phytoplankton physiological variability observations ( Figure 5), both showing a good agreement with the observations, with plausible values of P chl max and E k found through the water column particularly in terms of the vertical gradients. Calibration of these models against physiological 205 data is relatively novel as such comparisons remain rare. Greater complexity allowed the S2P3 v8.0 model to resolve a more diverse range of biogeochemical dynamics, explicitly accounting for zooplankton biomass and for the dynamics of internal quotas of phytoplanktonic cells, with phytoplankton biomass being in carbon, nitrogen, and chlorophyll currencies, allowing the decoupling of nutrient uptake from carbon fixation (Klausmeier et al., 2004;Flynn, 2008;Bougaran et al., 2010;Bernard, 2011;Mairet et al., 2011;Ayata et al., 2013). Including additional parameters in models can add more unconstrained degrees 210 of freedom (Ward et al., 2010), but also allows for more parameter combinations and, therefore, more flexibility to constrain S2P3 v8.0 in order to reproduce observations. The new added parameters and variables in the S2P3 v8.0 model, were carefully chosen to allow the model to be constrained against additional data, specifically zooplankton abundances and photosynthetic physiological measurements, therefore, a higher complexity allowed better representations of the temporal dynamics. Additionally, despite having more sophisticated formulations of the ecosystem, the S2P3 v8.0 model continues to be a 1-D model, allowing multiple experiments to be run at the same time with relatively low computational cost.

Sensitivity studies
An analysis was performed to assess the sensitivity of the model to selected parameter values for the S2P3 v8.0 model. Each parameter listed in Table 1 was varied in turn in order to understand how sensitive the model is to those changes and the effect that they have on the modelled ecosystem dynamics at the CCS location. In each case parameters were varied from the best calibrated value by +50% and -50%. Sensitivity studies are important tools to improve the accuracy of shelf sea models (Chen 245 et al., 2013), but developing these analyses has to be done carefully in order to identify which processes are responsible for the observed model behaviour (Ward et al., 2013). A direct comparison was calculated in terms of the S2P3 v8.0 model attributes considering: the timing and magnitude of the spring phytoplankton bloom, and the total annual zooplankton biomass for the calibrated version of the model and each experiment, providing better insights on the effects that each parameter produces in the behaviour of each model. Table 2 only shows the experiments involving the maximum ingestion rate of phytoplankton (R m ), 250 mortality of zooplankton (m), and the maximum value of the Chl : N ratio (θ N max ), as they demonstrated to be the most significant parameters in terms of the sensitivity of the model according to the attributes calculated, with the rest of the experiments omitted for this discussion. The S2P3 v8.0 model is strongly influenced by NPZ parameters, with the maximum ingestion rate of phytoplankton (R m ) (Stegert et al., 2007) and zooplankton mortality rate (m) having the largest effect in the magnitude of the spring phytoplankton bloom (Figures A1a, A4a) and in the total annual zooplankton biomass (Figures A1b, A4b). The ing low values of surface chlorophyll-a during spring ( Figure A5). It is well known that zooplankton are key players in the biogeochemical cycling of carbon and nutrients in marine ecosystems (Beaugrand and Kirby, 2010;, influencing the export of organic matter to the deep ocean (González et al., 2009;Juul-Pedersen et al., 2010). Additionally, grazing responses comprise the dominant losses for phytoplankton in the ocean (Banse, 1994), influencing plankton stocks and primary production (Franks et al., 1986).

265
Representation of phytoplankton physiology had an important influence on the timing of the spring phytoplankton bloom (Table 2), with θ N max affecting the S2P3 v8.0 model the most in terms of this attribute of the model structure, showing less productive and delayed spring blooms when θ N max is lower. The spring zooplankton bloom coincides with the timing of the phytoplankton blooms. Ayata et al. (2013) demonstrated that taking into account photoacclimation and variable stoichiometry 270 of phytoplankton growth in marine ecosystem models, produce qualitative and quantitative differences in phytoplankton dynamics. Moreover, these quota formulations in S2P3 v8.0 were compared to the available dataset of physiological observations ( Figure A9). It is interesting to note that the sensitivity analysis of NPZ parameters produced differences in the physiological variables P chl max and E k , specially at the surface ( Figures A3, A6), suggesting that the predator-prey interactions are indirectly influencing phytoplankton physiology, presumably through feedbacks whereby zooplankton mediated nutrient cycling and 275 which subsequently have an effect on phytoplankton physiology. A key point from this study showed that the greater variety of data used in a model, the more constrained the parameter choices seem to get. 6 Comparison of overall model performance

Experiments
The four calibrated models were further tested by running each of the extended period 1965 -2015, to evaluate the statistics 280 of productivity, partitioned between spring and summer. Table 3 shows that for the S2P3 v7.0 model, on average, 69.2% of the annual phytoplankton production occurs during the spring phytoplankton bloom. On the other hand, for the S2P3-NPZ model, on average, only 37.8% of the annual production occurs during spring months, showing that the predator-prey relationship has a strong influence on the magnitude of the spring phytoplankton bloom every year. Additionally, the S2P3-Photoacclim model shows a very strong spring phytoplankton bloom, corresponding to 97.4% of the total annual NPP. Finally, for the S2P3 v8.0 285 model, only 67.4% of the annual production corresponds to the spring bloom period. It is clear that, on average, the least productive model overall was S2P3 v8.0, followed by the S2P3-NPZ, S2P3 v7.0, and S2P3-Photoacclim. This shows the impact and complexity that the predator-prey relationship has on the model dynamics, with the addition of explicit zooplankton and their grazing activity as one of the main losses of phytoplankton (Franks et al., 1986). On the other hand, the S2P3 v7.0 model has, on average, more total annual NPP than the S2P3-NPZ model, suggesting that the influence of a constant grazing rate is 290 not as strong in comparison to the one provided by the zooplankton grazing (NPZ framework), because the predator-prey rela-tionship can not be entirely represented in the S2P3 v7.0 model. Additionally, S2P3-Photoacclim had, on average, the highest total annual NPP, demonstrating that the phenotypic responses of phytoplankton to changes in the environment and the addition of variable stoichiometric ratios of N : C : Chl, allows for more efficient phytoplankton, which use the available resources to grow and fix carbon more efficiently.  Table 3. Comparison between the S2P3 v7.0, S2P3-NPZ, S2P3-Photoacclim, and S2P3 v8.0 models in terms of the total spring NPP, total summer NPP, and total annual NPP calculated from 1965 to 2015, including the mean, maximum, minimum, and STD values.

Conclusions
This study demonstrates that the combination of an NPZ framework, photo-acclimation, and flexible stoichiometry of phytoplankton in one model produces a better representation of the ecosystem based on the comparison to observations. This combined framework offers a novel and innovative improvement to the S2P3 v7.0 model, for application to the CCS location insights in studies about the effects of physical forcing through tides (Sharples, 2008), intra and inter-annual variations in me-305 terology  and other drivers on PP, and phytoplankton dynamics would be possible.
Appropriate parameterisations to represent shelf seas is a subject that should be further supported by fieldwork campaigns and future work should aim to include additional datasets with longer time-series (Friedrichs et al., 2007;Ward et al., 2010). For this study, constraining the seasonal cycle of the phytoplankton physiology is not possible due to the lack of physiological 310 observations during other periods of the year, furthermore, phytoplankton and zooplankton biomass datasets only include the years 2014 and 2015 but longer time-series would help to improve the model calibration. Many model parameters quantities are poorly constrained observationally mainly due to the fact that model state variables are highly integrated pools, which are affected by biotic and abiotic factors in the environment, making them difficult to be determined by in situ measurements (Fennel et al., 2001).

315
The tuning and sensitivity analysis performed in this work allows a better understanding of the ecosystem dynamics represented in the model and how it is influenced by each parameter. By considering the timing and magnitude of the spring phytoplankton bloom, the annual zooplankton biomass, and summer phytoplankton photosynthetic physiology as features to be calculated, a quantitative and clearer comparison between each model could be developed. Finally, the model calibration will never be 320 in perfect agreement with all the observations, particularly in this case, where only one type of phytoplankton and zooplankton were considered. Thus responses must represent some typical or average dynamics and cannot represent any effects of competition between types. Additionally no stages of zooplankton growth were taken into account, which might affect the predator-prey interactions (Wroblewski, 1982;Fennel, 2001). Despite such potential limitations, the S2P3 v8.0 model is able to reasonably represent the integrated behaviour of the mixture of species that inhabit the NW European Shelf Sea.

325
Code and data availability. S2P3 v8.0 The current version of model is available from https://doi.org/10.5281/zenodo.3600467 under the Creative Commons Attribution 4.0 International license. The exact version of the model used to produce the results used in this paper is archived on Zenodo (https://zenodo.org), as are input data and scripts to run the model and produce the plots for all the simulations pre-330 sented in this paper.

335
-/domain contains, location (latitude and longitude), tidal components, and bathymetry (total depth) for the Central Celtic Sea in the western English Channel (s12_m2_s2_n2_h_tim.dat).