WAP-1D-VAR v1.0: development and evaluation of a one-dimensional variational data assimilation model for the marine ecosystem along the West Antarctic Peninsula

The West Antarctic Peninsula (WAP) is a rapidly warming region, with substantial ecological and biogeochemical responses to the observed change and variability for the past decades, revealed by multi-decadal observations from the Palmer Antarctica Long-Term Ecological Research (LTER) program. The wealth of these long-term observations provides an important resource for ecosystem modeling, but there has been a lack of focus on the development of numerical models that simulate time-evolving plankton dynamics over the austral growth season along the coastal WAP. Here, we introduce a one-dimensional variational data assimilation planktonic ecosystem model (i.e., the WAP-1D-VAR v1.0 model) equipped with a model parameter optimization scheme. We first demonstrate the modified and newly added model schemes to the pre-existing food web and biogeochemical components of the other ecosystem models that WAP-1D-VAR model was adapted from, including diagnostic sea-ice forcing and trophic interactions specific to the WAP region. We then present the results from model experiments where we assimilate 11 different data types from an example Palmer LTER growth season (October 2002–March 2003) directly related to corresponding model state variables and flows between these variables. The iterative data assimilation procedure reduces the misfits between observations and model results by 58 %, compared to before optimization, via an optimized set of 12 parameters out of a total of 72 free parameters. The optimized model results capture key WAP ecological features, such as blooms during seasonal sea-ice retreat, the lack of macronutrient limitation, and modeled variables and flows comparable to other studies in the WAP region, as well as several important ecosystem metrics. One exception is that the model slightly underestimates particle export flux, for which we discuss potential underlying reasons. The data assimilation scheme of the WAP-1D-VAR model enables the available observational data to constrain previously poorly understood processes, including the partitioning of primary production by different phytoplankton groups, the optimal chlorophyll-to-carbon ratio of the WAP phytoplankton community, and the partitioning of dissolved organic carbon pools with different lability. The WAP-1D-VAR model can be successfully employed to link the snapshots collected by the available data sets together to explain and understand the observed dynamics along the coastal WAP.

1 Introduction

The West Antarctic Peninsula (WAP) has experienced significant atmospheric and surface ocean warming since the 1950s, resulting in decreased winter sea-ice duration, the retreat of glaciers, and changes in upper-ocean dynamics (Clarke et al., 2009; Cook et al., 2005; Henley et al., 2019; King, 1994; Meredith and King, 2005; Stammerjohn et al., 2008; Vaughan et al., 2003, 2006; Whitehouse et al., 2008). These climate-driven changes propagate through marine food webs by affecting physiology of individual organisms and the whole communities (Ducklow et al., 2007). Long-term observational efforts through the Palmer Antarctica Long-Term Ecological Research program (LTER; since 1991) have demonstrated a range of ecological and biogeochemical responses to changing environments, including phytoplankton (Montes-Hugo et al., 2009; Saba et al., 2014; Schofield et al., 2017), marine heterotrophic bacteria (Bowman and Ducklow, 2015; Ducklow et al., 2012; Kim and Ducklow, 2016; Luria et al., 2014, 2017), nutrient drawdown (Kim et al., 2016), and micro- and macrozooplankton (Garzio and Steinberg, 2013; Steinberg et al., 2015; Thibodeau et al., 2019).

The wealth of Palmer LTER time-series observations provides an important resource for ecological and biogeochemical modeling, and different types of modeling approaches have been developed to explore the WAP responses to climate change and variability. For instance, an inverse modeling study estimated the steady-state dynamics of the WAP food web by deriving snapshots of flows among different plankton functional types and higher trophic levels (Sailley et al., 2013). However, there has been less focus on numerical ecosystem models that simulate time-evolving plankton dynamics over the full austral growth season along the coastal WAP. Numerical ecosystem models provide estimates of key rate processes for which observations have been less frequently or seldom made compared to frequently measured stocks and rates. Despite its strengths, constructing an ecosystem model is a challenge due to the lack of a priori knowledge on model parameter values and incomplete understanding of ecological processes that should be explicitly presented in the model structure (Ducklow et al., 2008; Murphy et al., 2012). Due to many observational studies, a more robust, yet still incomplete, data-based picture is emerging of WAP food-web interactions and ecosystem dynamics, which could guide a development of the WAP-specific numerical ecosystem model.

Here, we introduce a one-dimensional (1-D) variational data assimilation model specific to the coastal WAP (i.e., the WAP-1D-VAR v1.0 model) that we develop by adapting an existing biogeochemical–planktonic model of different ocean basins (Friedrichs, 2001; Friedrichs et al., 2006, 2007; Luo et al., 2010). The WAP-1D-VAR model is compared against the roughly semi-weekly biophysical observations over the austral growth season near Palmer Station on Anvers Island, Antarctica (64.77^∘ S, 64.05^∘ W). The field data record the seasonal variations in the initiation, peak, and termination of phytoplankton blooms and other biogeochemical processes modulated by variations in surface light, mixed layer depth, and sea-ice cover. In the present study, we (1) describe the structure and schemes of the WAP-1D-VAR model in great detail, (2) evaluate the model performance and robustness using a variety of quantitative metrics, and (3) discuss the model applicability with regard to capturing the key WAP ecological and biogeochemical features using the data from an example growth season.

Figure 1 The WAP-1D-VAR v1.0 model is forced by five different physical forcing, denoted as a horizontal row across the top of the schematic of the ecosystem model. The ecosystem component incorporates 11 different prognostic state variables. Higher level and refractory dissolved organic matter (RDOM) are represented implicitly.

2 Model development and implementation

2.1 Model state variables

The WAP-1D-VAR v1.0 model (Fig. 1) is originally derived and modified from data-assimilative, ocean regional test-bed models of the Arabian Sea, the equatorial Pacific, and the Hawaii Ocean Time-series Station ALOHA (A Long-Term Oligotrophic Habitat Assessment) (Friedrichs, 2001; Friedrichs et al., 2006, 2007; Luo et al., 2010). The WAP-1D-VAR model simulates stocks and flows of C, N, and P through 11 different model prognostic state variables. The two size-fractionated phytoplankton compartments (i.e., diatoms and cryptophytes) and the two different zooplankton compartments (i.e., microzooplankton and krill) are separately simulated following the plankton functional types as in Sailley et al. (2013) and the observations of the phytoplankton community structure along the coastal WAP. Typically, the coastal WAP is associated with large phytoplankton accumulations dominated by large (>20 µm) diatoms, but nanoflagellates (<20 µm) or cryptophytes are also an important component of the food web (Schofield et al., 2017). Mixed flagellates, prasinophytes, and type-4 haptophytes are also found in the region, but we choose to model only diatoms and cryptophytes, in order to avoid too many free (optimizable) parameters associated with each phytoplankton group. The third most dominant species is mixed flagellate but little is known about these taxa in the region and these taxa generally exhibit low interannual variability (Schofield et al., 2017). Functional grazing relationships are defined in which diatoms are consumed by both krill (Euphausia superba) and microzooplankton (mostly ciliates and other protozoa), cryptophytes are consumed by microzooplankton, and microzooplankton are grazed by krill. Other abundant zooplankton taxa in the WAP, such as salps, pteropods, and copepods (Steinberg et al., 2015), are not explicitly simulated in the WAP-1D-VAR model, in part to limit the model complexity and in part because of the limited data constraints on these groups, especially feeding. Higher trophic levels are implicitly represented to close the model. The WAP-1D-VAR model allows for the partitioning between labile dissolved organic matter (LDOM) and semi-labile DOM (SDOM) such that the entire LDOM pool is available but only a limited portion of the SDOM is available for bacterial utilization to account for lower lability of SDOM. Refractory DOM (RDOM) is not explicitly modeled due to its much longer turnover time than labile and semi-labile pools, but some mass flows are included to RDOM from other prognostic model compartments, such as bacteria, krill, and SDOM, to account for loss terms for those state variables. Detritus represents an average particulate organic matter (POM) pool after removing living phytoplankton and bacterial biomass, and sinking of the detritus pool contributes to particle export flux. The WAP-1D-VAR model explicitly simulates NO₃, NH₄, and PO₄ for inorganic (macro)nutrient compartments, but there is not a separate Fe model compartment or Fe uptake processes, given that even during the peak of the blooms iron is still measurable and iron limitation is absent or occurs only minimally and seasonally in the nearshore Palmer Station area (Carvalho et al., 2016; Sherrell et al., 2018).

2.2 Model equations

Here, we demonstrate key model processes that are either based on the existing schemes or built as new schemes for the coastal WAP region. The original model schemes are detailed in Supplementary Material of Luo et al. (2010). The WAP-1D-VAR v1.0 model simulates biological–physical model processes for a 1-D vertical water column, solving numerically for a discretized version of the time rate of change for each model state variable. For a generic tracer variable C, the time rate of change equation takes the form (Glover et al., 2011)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(wC\right)+{\displaystyle \frac{\partial }{\partial z}}\left(K_z{\displaystyle \frac{\partial C}{\partial z}}\right)+J_{\mathrm{C}},\end{array}$$

where z is the depth, w is the vertical velocity (the sum of water motion and gravitational particle sinking), K_z is the turbulent eddy diffusivity, and J_C is the biological and biogeochemical net source and sink term for C (Appendix A Equations; Eqs. A42–A45, A82–A85, A138–A140, A164–A166, A193–A195, A199–A201, A205–A210, and A212–A214). The physical advection and mixing terms are discussed below in Sect. 2.3 and applied sequentially following the computation of the biological and biogeochemical terms J_C using a constant time step of 1 h. The contributions of the source sink terms J_C to the full time rate of change equations are constructed as a series of coupled ordinary differential equations, detailed in Appendix A (Sects. A1–A9), and solved using a second-order Runge–Kutta numerical integration scheme. The WAP-1D-VAR model simulates the dynamics of C, N, and P, but here we only focus on the presentation of the model C dynamics. The cellular molar (e.g., N$/$C, P$/$C) quota parameter values of most state variables are fixed (Table 1) and not submitted to the optimization and data assimilation procedure. To first order, most model physiological processes are affected by water temperature, including the maximum growth rates of phytoplankton, bacteria, and zooplankton and basal respiration rates of bacteria and zooplankton. The Arrhenius function is implemented to change these physiological rates as a function of water temperature (Eq. A1).

Table 1 Summary of the model parameter symbol and definition, initial parameter values (p₀) and optimized values (p_f) for optimizable parameters, the cost function gradient with regard to the optimized parameter ($\partial J/\partial p$), and prescribed values for fixed model parameters over the course of simulations. The parameter with “n/a” in the parentheses is an optimized parameter with a large relative uncertainty, while the parameter with values in the parentheses is a constrained parameter (optimized with a low relative uncertainty) with its upper and lower bounds. The uncertainties for these upper and lower bounds are calculated as $p_{\mathrm{f}}\times e^{\pm {\mathit{\sigma }}_{\mathrm{f}}}$, where p_f is the value of the constrained parameter and σ_f is the square roots of diagonal elements of the inverse of the Hessian matrix. The cost function gradient with regard to the optimized parameter ($\partial J/\partial p$) after data assimilation is defined as $\mathrm{\Delta }J/e^{\mathrm{\Delta }p}$ where e^Δp≈Δp for an infinitely small Δp.

The net change of phytoplankton (both diatoms and cryptophytes) C biomass is driven by gross growth, DOC excretion, particulate organic carbon (POC) production via aggregation, respiration, and grazing (Eqs. A42 and A82), the net change of their N and P biomass by gross growth, dissolved organic nitrogen (phosphorus) excretion, particulate organic nitrogen (phosphorus) production, and grazing (Eqs. A43–A44 and A83–A84). The net change of their chlorophyll a (Chl a) by gross growth, DOM excretion, and grazing (Eqs. A45 and A85). The WAP-1D-VAR model adapts a phytoplankton growth scheme with flexible stoichiometry, in which phytoplankton cells are allowed to accumulate and store more nutrients under light stress (Bertilsson et al., 2003; Droop, 1974, 1983; McCarthy, 1980). The phytoplankton C growth rate is limited by their cellular nutrient quota (Eqs. A2–A3 and A46–A47). Modified from Geider et al. (1998), phytoplankton nitrogen uptake decreases when their cellular N$/$C quota is higher than their reference (Redfield) ratio but not limited when lower than their reference ratio (Eqs. A5, A9, A49, and A53). The nitrogen consumption completely ceases when the phytoplankton cellular quota reaches their maximum allowable ratios and is additionally limited by the ambient NO₃ and NH₄ concentrations with a Monod function (Eqs. A11–A12 and A55–A56). NH₄ inhibition on NO₃ uptake of phytoplankton is modeled by assigning lower $k^{{\mathrm{NH}}_{\mathrm{4}}}$ compared to $k^{{\mathrm{NO}}_{\mathrm{3}}}$ (Table 1). The inhibition term does not exist for PO₄. The uptake scheme is similar for PO₄ (Eqs. A14 and A58), but PO₄ can be consumed in great excess of current needs (Armstrong, 2006). Such luxury uptake is modeled by assigning smaller maximum and minimum P quota, which acts to alleviate P limitation. The maximum photosynthesis rate decreases when the phytoplankton cellular quota is lower than their reference ratio and approaches zero near their minimum ratio (Eqs. A7 and A51). The Chl production decreases with lowering photosynthetic active radiation (PAR) and completely ceases in dark (Eqs. A15 and A59). Phytoplankton release LDOM via passive diffusion of the low molecular weights DOM (e.g., neutral sugars and dissolved free amino acids) with the same cellular elemental ratio as that of phytoplankton (Fogg, 1966, Bjørnsen, 1988, Biddanda and Benner, 1997; Eqs. A17–A19 and A61–A63). Phytoplankton also release L- and SDOM actively, in the form of carbohydrate, as 75 % of the labile (Eqs. A20, A24, A64, and A68) and 25 % of the semi-labile pools (Eqs. A21, A27, A65, and A71). This active DOM production enables phytoplankton to adjust their stoichiometry to approach their reference ratio. If cellular organic C is in excess, DOC is released on a timescale of 2 d, and if excess nitrogen (phosphorus), DON (DOP) is released on a timescale of 8 d (Eqs. A22–A23 and A66–A67). Diatoms are grazed by both microzooplankton and krill (Eqs. A34–A41), while cryptophytes are only grazed by microzooplankton (Eqs. A78–A81). Microzooplankton grazing functions are altered by assigning grazing limitation terms (ϵ) to provide a limit on diatom grazing and route more cryptophytes to microzooplankton (Eqs. A34 and A78), based on initial modeling attempts where elevated diatom Chl was not simulated due to their much stronger removal by microzooplankton than cryptophytes. In principle, optimization should be able to capture the elevated diatom Chl by adjusting free parameters unless (1) the right parameters are not adjusted and/or the baseline (non-optimized) parameters need significant adjusting, and/or (2) the model equations are not adequate even with the optimized parameters. In our initial modeling attempts, the model failed to simulate the elevated diatom Chl with varying sets of the model initial parameter values assigned to decouple diatoms from their grazers. Thus, grazing limitation terms (ϵ) are instead assigned to limit microzooplankton grazing on diatoms for modeling purposes, the implementation of which is not strictly based on the ecological evidence of prey switching or of zooplankton mortality thresholds.

The net change of bacterial biomass is driven by their gross growth (via L- and SDOM uptake; Eqs. A97–A99, A100–A101, and A106–A107), respiration (Eq. A110), S- and RDOM excretion (Eqs. A111–A113 and A120–A128), grazing (Eqs. A129–A131), and mortality due to viral attack (Eqs. A132–A134). The WAP-1D-VAR model allows both L- and SDOM as the substrate sources for bacteria, and bacterial nutrient quota lets the lability of SDOM variable for their selective utilization. All the LDOM pool is available, while only a limited portion of the SDOM pool is allowed for bacterial utilization, the degree of which is controlled by an optimizable parameter controlling the relative utilization of SDOM to LDOM, or SDOM lability (i.e., r_SDOC, Eq. A96, Table 1). Bacterial C growth is determined by their cellular quota and available L- and SDOC concentration (Eqs. A97–A98), in which the growth would be limited if bacterial cellular nitrogen (phosphorus) quota is smaller than their reference ratios (Eqs. A93–A94). Bacteria take up LDOM in the way that the ratio of labile dissolved organic nitrogen (LDON) (LDOP) to LDOC uptake is the same as the bulk N$/$C (P$/$C) ratio of the LDOM (Eqs. A100 and A106). Bacteria take up SDOM with higher N$/$P ratios to reflect that SDOM with higher N$/$P ratios is more labile (Eq. A98). The ratio of semi-labile dissolved organic nitrogen (SDON) to SDOC uptake by bacteria would vary between the bulk N$/$C of SDOM and the bacterial reference cellular quota (Eqs. A101 and A107). Bacteria are modeled to either take up or release NH₄ and PO₄ to maintain their stable and consistent stoichiometry (Kirchman, 2000). Bacteria take up NO₃ only if their cellular N$/$C ratio is smaller than their reference ratio (i.e., when bacteria are short of nitrogen), in order to reflect higher energetic cost of NO₃ uptake than NH₄, but the amount of NO₃ uptake is modeled to be no more than 10 % of N-specific bulk L- and SDOM uptake, and the sum of NO₃ and NH₄ uptake is modeled to be no more than N-specific bulk L- and SDOM uptake (Eqs. A102–A105). These limit the maximum NO₃ uptake rate and set the inhibition of NH₄ uptake on NO₃ uptake. Bacteria excrete RDOM by transforming LDOM to RDOM (Eqs. A111–A113). Bacteria also adjust their cellular stoichiometry by remineralizing NH₄ and PO₄ if carbon is in short (i.e., N and P in excess; Eqs. A114–A115) and by excreting SDOC if carbon is in excess (i.e., nitrogen and phosphorus are in short; Eqs. A123–A128). Bacteria are grazed by microzooplankton (Eqs. A129–A131), and a certain percentage of bacteria gets lost to LDOC pool due to viral attack (Eqs. A132–A134).

The net change of zooplankton (both microzooplankton and krill) biomass is driven by their gross growth (via grazing on prey; Eqs. A143–A145 and A169–A171), L- and SDOM excretion (Eqs. A146–A154 and A172–A180), respiration (Eqs. A157 and A183), POM production (Eqs. A158–A160 and A184–A186), and grazing (Eqs. A161–A163 and A190–A192). Microzooplankton C growth is supported by consuming cryptophytes and bacteria (Eqs. A143–A145), while krill carbon growth is supported by consuming diatoms and microzooplankton (Eqs. A169–A171). Both zooplankton compartments follow the Holling type-2 density-dependent grazing function with a preferential selection on different prey species (Eqs. A34, A38, A78, A129, and A161). Both zooplankton groups release a portion of the organic matter that they ingest as DOM via sloppy feeding and excretion (Eqs. A146–A148, A149–A151, A172–A174, and A175–A177) such that the ratio of the released DON (DOP) to LDOC is equivalent to the N$/$C (P$/$C) ratio of zooplankton. The amount of SDOC excretion is a function of the total carbon growth (Eqs. A149 and A175), while the amount of SDON (SDOP) excretion is also a function of the zooplankton cellular N$/$C (P$/$C) ratio relative to their reference ratio (Eqs. A150–A151 and A176–A177). Zooplankton adjust their body cellular quota by either releasing SDOM if carbon is in excess or by regenerating NH₄ or PO₄ if nitrogen or phosphorus is in excess (Eqs. A152–A156 and A178–A182), similar to the bacterial scheme. Respiration is formulated such that basal respiration is based on a portion of zooplankton biomass, while active respiration is based on a portion of their grazed C (Eqs. A157 and A183). Both zooplankton egest fecal matter as POM (Eqs. A158–A160 and A184–A186), but only krill additionally excrete RDOM with N$/$C and P$/$C similar to bacteria (Eqs. A187–A189). Microzooplankton are grazed by krill (Eqs. A161–A163), while krill are removed by implicit higher trophic levels (Eqs. A190–A192), similarly calculated as a bacterial mortality term rather than as an explicit grazing process.

The net change of detritus is driven by POM production by all phyto- and zooplankton compartments that is routed to detrital pool (Eqs. A30–A33, A74–A77, A158–A160, and A184–A186) and its dissolution (Eqs. A196–A198). An optimizable vertical sinking speed is assigned to detritus to derive particle export flux (i.e., particle export flux is equal to detrital concentration multiplied by particle sinking velocity, wnsv, Table 1). The detritus that is lost due to dissolution is incorporated to SDOM pool when it sinks (Eqs. A208–A210) before regenerated to inorganic nutrients, rather than directly regenerated from as the particulate form. The net change of LDOM is driven by LDOM excretion by all phyto- and zooplankton compartments (Eqs. A17–A20, A61–A64, A146–A148, and A172–A174) and the amount of bacterial mortality that is incorporated to LDOM due to viral attack (Eqs. A132–A134) and its uptake by bacteria (Eq. A97). The net change of SDOM is driven by SDOM excretion by all organisms (Eqs. A21–A23, A65–A67, A149–A154, and A175–A180), the amount of detrital dissolution (Eqs. A196–A198), uptake by bacteria (Eqs. A98–A99), and conversion to RDOM pool (Eqs. A202–A204). The conversion of SDOM to RDOM pool is a function of the stoichiometry of SDOM, in which the conversion process is slower for higher N$/$C and P$/$C of SDOM, to reflect that nitrogen- and phosphorus-enriched SDOM is more likely labile. A certain percentage of NH₄ is converted to NO₃ on a daily basis to represent a simple nitrification process in the model (Eq. A211).

2.3 Physical forcing

The WAP-1D-VAR v1.0 model is forced by mixed layer depth (MLD), PAR at the ocean surface, sea-ice concentration, water-column temperature, vertical velocity, and vertical eddy diffusivity, at a temporal resolution of 1 d. Temperature, sea ice, and vertical eddy diffusivity are set up at every vertical grid (depth) point.

MLD is determined based on a finite difference density criterion with a threshold value of $\mathrm{\Delta }{\mathit{\sigma }}_{\mathit{\theta }}=\mathrm{0.03}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (Montégut et al., 2004) after calculating potential density of water mass from temperature and salinity conductivity–temperature–depth (CTD) observations. Vertical velocity, w, is assigned as zero because it is very weak in the surface waters of the study site and materials are transported vertically mostly by diffusion. The vertical eddy diffusivity scheme treats the rapid vertical mixing in the surface boundary layer by homogenizing model state variables instantaneously in the mixed layer (i.e., by averaging at every time step). Thus, K_z value above MLD is not required, and only K_z below MLD is calculated as follows:

$$\begin{array}{}\text{(2)}& K_z\left(z\right)=K_{z_{\mathrm{0}}}\times \exp \left\{-\mathit{\alpha }\times (z-\text{MLD})\right\},\end{array}$$

where z is depth (m) below MLD and $K_{z_{\mathrm{0}}}$ is the vertical eddy diffusivity at the bottom of the mixed layer ($\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) (Klinck, 1998; Smith et al., 1999), and α is 0.01 (m⁻¹).

Daily surface downward solar radiation flux (National Centers for Environmental Prediction reanalysis daily averages) is used to calculate sea surface PAR. PAR is estimated as 46 % of the total solar radiation (Pinker and Laszlo, 1992, Kirk, 1994). The attenuation of PAR as a function of depth is calculated as follows:

$$\begin{array}{}\text{(3)}& \text{PAR}\left(z\right)={\text{PAR}}_{\mathrm{0}}\times \exp \left\{-(k_{\mathrm{w}}+k_{\mathrm{c}}\times \text{CHL})\times z\right\},\end{array}$$

where z is depth (m), PAR₀ is PAR level at sea surface (W m⁻²), k_w is the attenuation coefficient for seawater (m⁻¹), k_c is the attenuation coefficient for Chl ((mg Chl)⁻¹ m²), and CHL is the Chl concentration (mg Chl m⁻³).

Sea-ice conditions in the coastal WAP do not necessarily represent solely local temperature and climate conditions, given that sea ice can be impacted by temperature, mixed layer, heat fluxes, regional winds, and other physical processes (Saenz et al., in review). We implement a sea-ice model scheme to account for light transmission through sea ice (5 % of incident irradiance, as a typical transmittance value used in the Community Earth System Model) and non-linearities in the photosynthesis–irradiance (P–I) response under partial ice concentration (Long et al., 2015) using percent daily sea-ice concentration data (GSFC Bootstrap versions 2/3, derived from SMMR/SSMI satellite temperature brightness data binned into 25 by 25 km grid cells). In many previous models, the light-limitation term ℒ(I) is calculated as a function of mean irradiance $\stackrel{\mathrm{‾}}{I}$ averaged over both ice-covered and open-water conditions, so $\mathcal{L}\left(\stackrel{\mathrm{‾}}{I}\right)$; instead, we compute the mean of light-limitation term ($\stackrel{\mathrm{‾}}{\mathcal{L}\left(I\right)})$ as a function of fractional sea ice and open-water and incident irradiance:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathcal{L}\left(I\right)=P^{\mathrm{C}}/P_{\text{MAX}}^{\mathrm{C}}=\mathrm{1}-\exp (-I/I_{\mathrm{k}})\text{(5)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathcal{L}\left(I\right)}=f_{\mathrm{i}}\times \mathcal{L}\left(I_{\mathrm{i}}\right)+f_{\mathrm{o}}\times \mathcal{L}\left(I_{\mathrm{o}}\right),\end{array}$$

where P^C is the C-specific photosynthetic rate (d⁻¹), $P_{\text{MAX}}^{\mathrm{C}}$ is the maximum photosynthetic rate (d⁻¹), I_k is the parameter describing the light-saturation behavior of the P–I curve (W m⁻²), I_o is the open-water irradiance, I_i is the under-ice irradiance (i.e., $I_{\mathrm{i}}=\mathrm{0.05}\times I_{\mathrm{o}}$), f_i is the fraction of area covered with sea ice, and f_o is the fraction of open water (i.e., $f_{\mathrm{o}}=\mathrm{1}-f_{\mathrm{i}}$).

2.4 Variational data assimilation

The WAP-1D-VAR v1.0 model is equipped with a built-in data assimilation scheme based on a variational adjoint method (Lawson et al., 1995). This method generates optimal model solutions that minimize the difference between model results and observations by objectively optimizing model parameter values (Friedrichs, 2001; Spitz et al., 2001; Ward et al., 2010). In detail, the assimilation scheme (Fig. 2) consists of four steps (Glover et al., 2011): (1) starting with initial values of the model parameters (see below), the model is integrated forward in time from specified initial conditions to calculate the difference between the model simulation and the field data, or the model–observation misfit (i.e., cost function; Sect. 2.5, Eq. 6); (2) an adjoint model constructed using the Tangent linear and Adjoint Model Compiler (TAPENADE) is integrated backward in time to compute the gradients of the total cost with respect to the model parameters; (3) the computed gradients are then passed to a limited-memory quasi-Newton optimization software M1QN3 3.1 (Gilbert and Lemaréchal, 1989) to determine the direction and optimal step size by which the selected model parameters (see below) need to be modified in order to reduce the total cost; and (4) a new forward-in-time simulation is conducted using the new set of modified (optimized) parameter values. These four-step procedures are conducted in an iterative manner until the preset convergence criteria (i.e., low gradients of the total cost function with respect to optimized parameters and positive eigenvalues of the Hessian matrix) are satisfied to ensure that the optimized parameters converge and the total cost function reaches a local minimum.

Figure 2 A variational adjoint scheme is employed for the parameter optimization and data assimilation processes (adapted from Glover et al., 2011). Gradient: the sensitivity of the total cost function with respect to model parameter from optimization.

Initial values of the model parameters (total of 72 free or optimizable parameters, Table 1) are assigned based on literature values (Caron et al., 2000, Luo et al., 2010, Garzio et al., 2013) without examining the effects of the initial parameter values on the model results prior to optimization. As is typical for many types of ecosystem models, a collection of what appear to be reasonable initial parameter estimates can result in relatively poor overall system behavior because of system-level interactions of different model components. In most marine ecosystem models, these initial parameter values are subjectively and manually adjusted to improve the simulation, and the simulations with the initial, unadjusted parameter values are rarely shown. However, here with a more objective optimization approach that we conduct, the initial and optimized solutions can be explicitly compared (Sect. 4). Optimization starts by submitting a subset of the 72 free model parameters rather than submitting all of them at once. This initial parameter subset consists of 10 different model parameters, the change of which yields the largest decrease in the total cost function (i.e., which also happens to be usually one per each state variable). These include α_DA (initial slope of photosynthesis vs. irradiance curve of diatoms, $\mathrm{mol}\mathrm{C}\left(\mathrm{g}\mathrm{Chl}\mathit{\text{a}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}\right(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$), α_CR (initial slope of photosynthesis vs. irradiance curve of cryptophytes, $\mathrm{mol}\mathrm{C}\left(\mathrm{g}\mathrm{Chl}\mathit{\text{a}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}\right(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$), Θ (maximum Chl$/$N ratio, g Chl a (mol N)⁻¹), μ_BAC (maximum bacterial growth rate, d⁻¹), $r_{\text{max,BAC}}^{\mathrm{A}}$ (maximum bacterial active respiration rate, d⁻¹), g_BAC (half-saturation density of bacteria in microzooplankton grazing, mmol C m⁻³), μ_MZ (maximum microzooplankton growth rate, d⁻¹), μ_KR (maximum krill growth rate, d⁻¹), and remv_KR (krill removal rate by higher trophic levels ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$; Table 1).

When computed at the minimum of the cost function value, the inverse of the Hessian matrix provides the uncertainties of optimized parameters, cross-correlations among parameters, and sensitivities of the total cost function to each parameter (Matear, 1996; Tziperman and Thacker, 1989). High off-diagonal values in the inversed Hessian matrix indicate highly cross-correlated model parameters, so one of the highly cross-correlated parameters is removed from the optimization. The square root of a diagonal element in the inversed Hessian matrix is the logarithm of the relative uncertainty (σ_f) of the corresponding optimized parameter. The absolute uncertainty of the constrained parameter is calculated as $p_{\mathrm{f}}\times e^{\pm {\mathit{\sigma }}_{\mathrm{f}}}$ where p_f is the value of the optimized parameter (Table 1). If optimized to ecologically unrealistic values, parameters are kept back to their respective initial values and removed from the next optimization cycle. Optimized parameters with σ_f larger than 50 % are updated but removed from the next optimization cycle (i.e., defined as “optimized” parameters), while optimized parameters with σ_f smaller than 50 % are updated and kept for the next optimization cycle (i.e., defined as “constrained parameters”). Constrained parameters are reported with the uncertainties, while optimized parameters are reported without the uncertainties (Table 1) because both changed parameters consist of an optimized model parameter set, but the parameters reported with the uncertainty ranges are the ones optimized with relatively small uncertainties and considered constrained. This way, a part of the initial parameter subset forms a final optimized parameter set. The gradients of the total cost function with respect to all 72 parameters are then evaluated, the parameters with large gradients (e.g., >5) are re-submitted to optimization to further reduce the total cost, the gradients are evaluated again, and these cycles repeat until the termination of optimization. Optimization terminates when the gradients are reasonably low (e.g., $<{\mathrm{10}}^{-\mathrm{2}}$ for constrained parameters, <5 for optimized parameters, and <10 for unoptimized parameters). This final optimized model parameter set forms the basis of the results presented throughout this study (Sect. 4). Additionally, in order to assess the sensitivity of the model optimization results with regard to the initial parameter choice, we perturb by ±50 % a subset of the initial parameter values used in the reference (original) optimization experiments to form different initial parameter sets (a total of 15 sets consisting of partially or fully perturbed 18 parameters, Tables B1–2) and conduct new optimization experiments from each set (Sect. 4.1).

2.5 Cost function

To represent a misfit between observations and model output, a total cost function is calculated as follows (Luo et al., 2010):

$$\begin{array}{}\text{(6)}& J=\sum _{m=\mathrm{1}}^M{\displaystyle \frac{\mathrm{1}}{N_m}}\sum _{n=\mathrm{1}}^{N_m}{\left({\displaystyle \frac{{\widehat{a}}_{m,n}-a_{m,n}}{{\mathit{\sigma }}_m}}\right)}^{\mathrm{2}},\end{array}$$

where m and n represent assimilated data types and data points, respectively, M and N_m are the total number of assimilated data types and data points for data type m, respectively, σ_m is the target error for data type m, a_m,n is observations, and ${\widehat{a}}_{m,n}$ is model output. Given the high biological productivity of the WAP waters and the approximate log-normal distribution of many marine biological variables, the base-10 logarithms of Chl and primary production (PP) are used in the cost function calculation to capture phytoplankton dynamics (Campbell, 1995; Glover et al., 2018). The target error is calculated for each data type as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_m=\stackrel{\mathrm{‾}}{a_{m,n}}\times {\text{CV}}_m,\end{array}$$

where $\stackrel{\mathrm{‾}}{a_{m,n}}$ is the climatological mean (over the select nine growth seasons; see below) of the observations and CV_m is the averaged coefficient of variation (CV) of the observations of each data type in the mixed layer (due to observational error and seasonal and interannual variations) calculated using all of the observational data over nine growth-season periods between 2002–2003 and 2011–2012, except the 2007–2008 growth season due to its missing data. These nine growth seasons are chosen, instead of the multi-decadal observations available from Palmer LTER (since 1991), due to the relatively more complete data coverage in those seasons. The standard deviations are used as target errors of the log-converted data types. The CV of the log-converted data type is estimated as the average of ±1 standard deviation in log space converted back into normal space (Doney et al., 2003; Glover et al., 2018). Hereafter, we present the total cost normalized by M (J equivalent to $J/M$ hereafter) as it indicates the model–observation misfit equivalent to a reduced chi-square estimate of model goodness of fit. We report the normalized total cost J along with normalized costs of individual data types throughout this study. J=1 indicates a good fit, J≫1 indicates a poor fit or underestimation of the error variance, and J≪1 indicates an overfitting of the data, fitting the noise, or overestimation of the error variance.

3 Model experiments

3.1 Modeling framework

To examine the applicability of the WAP-1D-VAR v1.0 model to the coastal WAP region, we select a nearshore Palmer LTER water-column time-series station, Station E (64.77^∘ S, 64.05^∘ W), as the modeling site that is ∼200 m deep and situated approximately 3 km south of Palmer Station and 6.5 km northeast of the head of Palmer Deep (Sherrell et al., 2018). Physical forcing (Fig. 3) and data types assimilated are derived from roughly semi-weekly physical, chemical, and biological profiles collected from small boat via a profiling CTD and discrete water samples at Station E. When weather and ice conditions permit, water-column sampling at the station has been conducted twice a week over the growth season. Seven upper-ocean layer depths (2.5, 10, 20, 30, 40, 50, and 60 m) are chosen for the model vertical grids. The model depth can be extended to as deep as needed, but this study is focused to upper 60 m water column to fully take advantage of the large data availability. Also, conceptually, the application of the 1-D model framework makes the most sense for the upper water column dominated by local seasonal processes, and extension of the model into deeper water well below the maximum seasonal mixed layer becomes more problematic because of the growing importance of lateral advective processes that are not well captured in the 1-D model framework. The vertical structure of the water column can be affected by growing sea ice due to reduced wind-driven turbulence and brine rejection during winter, but this is what a prognostic, coupled ocean–ice 1-D model can offer to simulate, not our diagnostic-forcing-based model that was used in this study. Also, because our model simulates only the spring–summer growth season, the impact of winter sea ice on ecosystem dynamics is less of a concern.

Figure 3 Physical forcing used in the WAP-1D-VAR v1.0 model, including surface photosynthetically active radiation (PAR) (a), sea-ice concentration (b), water temperature (c), and vertical eddy diffusivity (d) overlaid with mixed layer depth (MLD; dashed line) in the modeled growth season of 2002–2003.

Given the routine observations of Palmer LTER available over the growth season (October–March), we simulate one example growth season with the most complete data coverage, from October 2002 to March 2003 (2002–2003 growth season hereafter), instead of a series of different growth seasons in a continuous manner. The example growth season simulations utilize this year's specific observed physical forcing fields and assimilated biological and biogeochemical observations. Each Palmer LTER growth season should be modeled to have its own unique optimized parameter set, as well as initial conditions and physical forcing that together determine the model solution for that year; however, only the 2002–2003 growth season simulations are modeled in this study for model analysis and evaluation.

3.2 Initial and boundary conditions

Model initial conditions are prescribed 135 d before the model start date for the growth season (15 October 2002) on 1 June 2002. This 135 d spinup is conducted to minimize the impact of initial conditions on the model output over the growth season. Initial conditions are prepared by optimizing the full growth seasonal cycle forced by climatological physics and assimilated with climatological observations and with the same bottom boundary conditions used in the optimization of the 2002–2003 growth season (i.e., climatological model; using climatological physics and observations averaged over nine growth-season periods between 2002–2003 and 2011–2012 except the 2007–2008 growth season due to its missing data). For the first climatological model simulation, initial conditions are prepared by adjusting manually following literature values (e.g., Luo et al., 2010). Due to strong interannual variability in the phytoplankton bloom phenology at Palmer Station, averaging across all these 9 years does not reflect distinct seasonal phytoplankton peaks, leading to underestimated phytoplankton values (not shown). To capture this non-linear aspect of the coastal WAP system, we construct the climatological year by applying a single time shift to all variables so that a seasonal PP peak of each year lines up with an average date of seasonal PP peaks from all years. Most biological initial conditions on 1 June are close to zero given the lack of active physiological processes in the very low light and the presence of sea ice during wintertime before the model growth season starts. All the data types are set to zero at the lower boundary (bottom) except for NO₃, PO₄, SDOC, SDON, and SDOP in which the climatological values at 65 m are used for lower boundary values (25.9, 1.9, 6.5, 0.6, and 0.03 mmol m⁻³, respectively).

3.3 Assimilated data

We include the data types directly related to corresponding model outputs, including a mix of ecosystem stocks or state variables – NO₃, PO₄, Chl for diatoms and cryptophytes, bacterial biomass, microzooplankton biomass, SDOC, POC, and particulate organic nitrogen (PON), as well as carbon flows among model stocks – bulk net PP and bacterial production (BP). These data sets have been sampled semi-weekly at Palmer Station E (64.77^∘ S, 64.05^∘ W), the same location where our model is set up, and are available from the Palmer LTER data website (see data availability). The distinction between diatoms and cryptophytes is established by assimilating phytoplankton taxonomic-specific Chl data for diatoms and non-diatom species derived from a high-performance liquid chromatography (HPLC) and CHEMTAX analysis (Schofield et al., 2017), but given cryptophytes being the second dominant species in the water samples at the study site, cryptophytes are assumed to represent all non-diatom species for modeling purposes. Given that POC (PON) from bottle filtration may capture both living biomass and detrital material, we adjust the observed POC (PON) by subtracting phytoplankton and bacterial C (N) biomass to estimate the detrital pool, in order to only include non-living particles to detrital pool. When phytoplankton or bacterial biomass data are not available, we assign climatological (2002–2003 to 2011–2012) fractions of POC (PON) to detrital pool. Phytoplankton and bacterial biomass accounts for 26 % of total POC and 29 % of total PON. In converting Chl to phytoplankton carbon (nitrogen) biomass, the maximum Chl$/$C (Chl$/$N) ratio submitted for optimization is used along with other reference ratios (Table 1). Microzooplankton biomass data are not available for the full time series, so their data from grazing experiments at Palmer Station (Garzio et al., 2013) are assimilated to at least provide constraints on bacterial and cryptophyte grazing processes. However, due to the discrepancy in the timing and location from model simulations of this study, the microzooplankton model–observation misfits are not analyzed in the present study. Krill biomass data are not assimilated due to the strong patchiness of their distribution that may hinder proper model optimization. The vertical profiles of most of the data types are assimilated, whereas average NO₃ and PO₄ concentrations in the mixed layer are assimilated due to the difficulty of simulating depth-dependent nutrient concentrations and the fact that net PP is mostly determined by surface nutrient concentrations (Luo et al., 2010). BP ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$) is derived from the ³H-leucine incorporation rate ($\mathrm{pmol}{\mathrm{l}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) data using the conversion factor of 1.5 kg C (mol leucine)⁻¹ incorporated (Ducklow, 2000). Bacterial biomass (mmol C m⁻³) is estimated from bacterial abundance measured by flow cytometry with the conversion factor of 10 fg C cell⁻¹ (Fukuda et al., 1998). SDOC is calculated by subtracting the background concentration (41.2 mmol m⁻³ for the modeling site) from total DOC concentration.

3.4 Uncertainty analysis

Uncertainties of the optimized parameters are computed from a finite difference approximation of the complete Hessian matrix (i.e., second derivatives of the cost function with respect to the model parameters) during the iterative optimization process (details in Sect. 2.4). We then conduct Monte Carlo experiments to calculate the impact of the optimized parameter uncertainties on the model results. We first create an ensemble of parameter sets (n=1000) by randomly sampling values within the uncertainty ranges of the constrained parameters and then perform a model simulation using each parameter set. A total of 1000 Monte Carlo experiments were shown to be adequate from a series of tests with different numbers of Monte Carlo sampling (n=500–2000), where standard deviations of model-simulated values converged after >1000 Monte Carlo sampling (not shown). All uncertainty estimates are calculated following standard error propagation rules and presented as ±1 standard deviation in the study.

4 Results and discussion

4.1 Model skill assessment

In the case of the example growth season (2002–2003) modeled in this study, the iterative data assimilation–parameter optimization procedure reduced by 58 % the misfits between observations and model output compared to the misfits obtained from the initial parameter values (Table 2). The optimized model solution satisfied the preset convergence criteria, with the low gradients of the total cost with respect to the optimized parameters and positive eigenvalues of the Hessian matrix. Notably, this was achieved by optimizing a subset of 12 (9 constrained and 3 optimized) parameters among the total of 72 optimizable parameters (Table 1, Sect. 4.2). To examine the sensitivity of the optimized model solution to the initial parameter choice, a series of new optimization experiments (n=15) was conducted with a varying subset of the initial parameter values perturbed by ±50 % of those used in the original optimization experiment (Table B1). These experiments showed that the optimized model results (i.e., the reference case; Table 1) were not sensitive to the initial choice of the parameters. The 15 different initial parameter sets resulted in a range of initial model–observation misfits, some substantially larger than the reference case (14.25–28.24 vs. 14.85 for the reference case). However, the total normalized optimized cost values of the 15 sensitivity experiments (5.79–7.19) were similar to that of the reference case of 6.42. In sensitivity experiment no. 12, the initial model–observation misfit was ∼2 times larger than that of the reference case, and there was up to 76 % of the reduction in the model–observational misfit (vs. 58 % of the reduction in the reference case; Table B1, Table 1). These results suggest that no matter where in parameter space the optimization process starts from, the optimization scheme takes the model cost function to similar local minima. Importantly, this was achieved by similar subsets of the optimized parameters: μ_DA, μ_CR, $r_{\text{max,BAC}}^{\mathrm{A}}$, and μ_MZ were optimized in all cases, while α_DA, α_CR, Θ, μ_BAC, g_CR, $g_{\text{DA}}^{\prime }$, g_BAC, μ_KR, g_MZ, and remv_KR were optimized except for a few cases (Table B1). The uncertainties of the optimized parameters were also similar among different optimizations, with most of the relative errors <0.5. Constrained parameter values and their uncertainty ranges averaged over the sensitivity experiments (Table B2) were comparable to those in the original optimization experiments (Table 1).

Table 2 The observed mean ($\stackrel{\mathrm{‾}}{a}$), coefficient of variation (CV), and target error (σ) of each assimilated data type used for calculating the cost function before and after optimization. J₀ is the normalized cost function before optimization and J_f is the normalized cost function after optimization (Eq. 6). Data type units: mmol m⁻³ for NO₃, PO₄, bacterial biomass, SDOC, and POC; mg m⁻³ for diatom Chl, cryptophyte Chl, mmol N m⁻³ for PON; and $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$ for PP and BP. The average error (ε_bias) of each data type (for non-transformed or raw Chl and PP) is calculated from Stow et al. (2009) before and after optimization where a positive value indicates the model overestimation of the observation, and vice versa.

Overall, there was a good model–data fit with the largely decreased cost value for each data type after optimization (Table 2). Optimization yielded J_f close to 1 for all data types, compared to the initial model solution where three data types – diatom Chl, crypto Chl, and bacterial biomass – had particularly poor model fits to observations and underestimated error variances (J≫1). Compared to the initial (unoptimized) model results, the average errors (ε_bias, Doney et al., 2009; Stow et al., 2009) in the optimized model results indicated that diatom Chl, cryptophyte Chl, bacterial biomass, BP, and POC had reduced model biases, while NO₃, PO₄, PP, SDOC, and PON had increased model biases (for both positive and negative biases, defined as ε_bias>0 and ε_bias<0, for model overestimation and underestimation of the observation, respectively). Optimization resulted in the negative model bias for PO₄, compared to the positive model bias in the initial model results. The point-to-point comparison plots showed that there were seasonally consistent, negative model biases for PP, POC, and PON (Fig. B1). Model skill was further evaluated with a Taylor diagram (Taylor, 2001) summarizing the statistics of the correlation coefficient between model output and observations, normalized standard deviation (by the standard deviation of the observations), and centered (bias removed) root-mean-square difference (RMSD) for each data type, in which a better model skill is characterized by a higher correlation, a normalized standard deviation close to 1, and a lower RMSD (Fig. 4). Optimization resulted in better model skills for cryptophyte Chl, PP, BP, and bacterial biomass via increased correlation coefficients and lowered RMSD (Fig. 4b), compared to those in the unoptimized model results (Fig. 4a). After optimization, the normalized standard deviations of PP, BP, bacterial biomass, phosphate, POC, and PON were closer to 1 (Fig. 4b). Direct comparisons with the observational data showed that the optimized model parameter set captured better the increases in diatom biomass early in the season, cryptophyte biomass in January, and bacterial biomass in mid-February, compared to the unoptimized model parameter set (Figs. 5a and b and B2).

Figure 4 Model skill assessment: the Taylor diagrams using a polar-coordinate system summarizing the model–observation correspondence for each model stock and flow for the modeled growth season of 2002–2003 before (a) and after optimization (b). The angular coordinate for the Pearson correlation coefficient (r), the distance from the origin for the standard deviation normalized by the standard deviation of the observation, and the distance from point (1,0), marked as REF on x axis, for the centered (bias removed) root-mean-square difference (RMSD) between model results and observations.

Figure 5 The model state variables for the modeled growth season of 2002–2003 (x axis; month/day) for before (a) and after optimization (b) with init. as the initial (unoptimized) model results and opt. as the optimized model results. The error (standard deviation) of each model state variable from the Monte Carlo experiments (n=1000) is available in Fig. B3.

4.2 Optimized parameters

The number of the optimized parameters in this study is small and comparable to those from other data-assimilative models focused on different marine environments (Friedrichs, 2001; Friedrichs et al., 2006, 2007; Luo et al., 2010). This is consistent with the general behavior of marine plankton ecosystem models, in which well-posed model solutions would be found with only a subset of independent model parameters due to many cross-correlated parameters inherent in non-linear model equations (Fennel et al., 2001; Harmon and Challenor, 1997; Matear, 1996; Prunet et al., 1996a, b). Ecosystem models with a relatively large number of unconstrained parameters (i.e., equivalent to the optimized parameters with high uncertainties in the present study) might reduce total costs to a greater extent but could possess low predictive skill as a result of being overtuned to fit noise in the observations (Friedrichs et al., 2007). Also, there are several field- and lab-based studies at the study site or in a similar polar environment that reported the values of the model parameters used in the WAP-1D-VAR v1.0 model, including the bacterial growth rate of 0.82 d⁻¹, total phytoplankton (including large cells like diatoms) growth rate of 0.33–0.55 d⁻¹, nanophytoplankton (corresponding to cryptophytes) growth rate of 0.52–0.99 d⁻¹ (Garzio et al., 2013), and the microzooplankton growth rate of up to 1.0 d⁻¹ (Caron et al., 2000). The optimized values of the maximum bacterial, diatom, cryptophyte, and microzooplankton growth rates in this study were 1.06 d⁻¹ (0.93–1.20 d⁻¹), 0.77 d⁻¹ (0.68–0.88 d⁻¹), 0.72 d⁻¹ (0.61–0.85 d⁻¹), and 1.18 d⁻¹ (1.10-1.26 d⁻¹), falling in the ranges of those measured bacterial, total phytoplankton, nanophytoplankton, and microzooplankton growth rates, respectively.

The ensemble of the model state variables and flows obtained from the Monte Carlo experiments had generally small standard variations at each model time step and grid, suggesting the robustness of the modeled fields against the variations in the optimized parameter values (Figs. B3 and B4). To examine the translation of the optimized parameter values to altered functioning of the WAP biogeochemical processes, we compared two different sets of the model simulation results – one based on the initial parameter values (Figs. 5a, 6a, and 7a) and the other based on the optimized parameter values (Figs. 5b, 6b, and 7b). However, due to the non-linearities in the model, it is not straightforward to identify what causes the parameter variations, except for a few cases in which the changes in the parameter values are clearly linked to the difference in the model state variables and flows. The first case is the relation of the increased g_BAC value (bacterial half-saturation concentration in microzooplankton grazing, mmol C m⁻³) to the elevated bacterial accumulations after optimization (Table 1, Figs. 5 and 7). The second case is the link between Θ (maximum Chl$/$N ratio, g Chl a (mol N)⁻¹) and the relative dominance of cryptophytes in total phytoplankton accumulations. It has been demonstrated that the variations of Θ are driven by an imbalance between the rate of light absorption and energy demands for photosynthesis and biosynthesis in phytoplankton cells (Geider et al., 1997). Θ can also change because of the variations in phytoplankton photo-acclimation or physiological differences across phytoplankton groups, from a lower Θ value for smaller species to a higher Θ value for larger diatom cells (Geider, 1987). Θ was optimized to a 30 % lower value than the initial parameter value (Table 1), in order to simulate the relatively larger proportion of cryptophytes in total phytoplankton accumulations in the optimized model results compared to the unoptimized model results (Figs. 5 and 7). By contrast, the remaining cases are not as clear because the first-order impact of parameter variations on the model results is less direct and more nuanced. Compared to the unoptimized results, the decreases in μ_DA (diatom C-specific maximum growth rate, d⁻¹), μ_CR (cryptophytes C-specific maximum growth rate, d⁻¹), α_DA (initial slope of P–I curve of diatoms, $\mathrm{mol}\mathrm{C}\left(\mathrm{g}\mathrm{Chl}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}\right(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$), and α_CR (initial slope of P–I curve of cryptophytes, $\mathrm{mol}\mathrm{C}\left(\mathrm{g}\mathrm{Chl}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}\right(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}$) did not lead to decreased diatom and cryptophyte accumulations, presumably due to decreased g_MZ (microzooplankton half-saturation concentration in krill grazing, mmol C m⁻³) and increased remv_KR (krill removal rate by higher trophic levels, $(\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$) after optimization (Table 1, Figs. 5 and 7). Similarly, the decreased μ_BAC (maximum bacterial growth rate, d⁻¹) and the increased $r_{\text{MAX,BAC}}^{\mathrm{A}}$ (bacterial maximum active respiration rate, d⁻¹) did not lead to decreased bacterial accumulations, presumably due to the increased g_BAC (bacterial half-saturation concentration in microzooplankton grazing, mmol C m⁻³) and the decreased g_MZ (microzooplankton half-saturation concentration in krill grazing, mmol C m⁻³).

Figure 6 Model ecosystem indices before (a) and after optimization (b). The key ecosystem indices for the modeled growth season of 2002–2003 (x axis; month/day), with init. as the initial (unoptimized) model results and opt. as the optimized model results. NPP: net primary production, NCP: net community production, C export flux: particulate organic carbon (POC) export flux, and BP: bacterial production. The error (standard deviation) of each rate process from the Monte Carlo experiments (n=1000) is available in Fig. B4.

Figure 7 Depth-integrated (0–60 m) carbon stocks (mmol C m⁻³), flows ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$), and POC export flux at 60 m ($\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) averaged over the modeled austral growth season (October 2002–March 2003) before (a) and after optimization (b). Values in parentheses as uncertainties from season averaging (a) and as uncertainties propagated from season averaging and depth integration of the Monte Carlo errors for panel (b).

4.3 Ecosystem indices

We calculated key ecosystem indices for the modeled growth season, including NPP (directly comparable to ¹⁴C–PP observations), net community production (NCP; i.e., NCP=NPP – bacterial, microzooplankton, and krill respiration), BP, and POC export (sinking) flux (Fig. 6). Setting an upper limit for lateral or vertical carbon export from the euphotic zone (Dugdale and Goering, 1967), over appropriate timescales and space scales, NCP is quantitatively equivalent to new production that is supported via external sources of nitrogen (Ducklow and Doney, 2013). In both optimized and unoptimized model results, NPP increased after complete sea-ice retreat, but a brief ice-edge bloom was simulated under sea ice at the beginning of the growth season (Figs. 3 and 6). Seasonal patterns of NCP resembled those of NPP and occasionally fell below zero (i.e., the net heterotrophy) in subsurface waters for both optimized and unoptimized cases (Fig. 6). The POC export flux increased over time and reached the maximum value at the end of the growth season in both model results, but there were two major POC flux events separated by weaker, in-between flux events in December in the optimized results that the initial model results did not capture (Fig. 6b). After optimization, the correlation coefficients adjusted from 0.88 to 0.36, 0.89 to 0.68, and 0.45 to 0.73 for the NPP vs. NCP pair, the NPP vs. BP pair, and the NCP vs. POC export flux (lagged by 30 d) pair (all p<0.001). In the optimized model results, the growth-season mean of the depth-integrated NPP, NCP, and BP in the 60 m water column, and the 30 d lagged POC export flux at 60 m were 19±8 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, 10±3 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, 1±1 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and 2±0.3 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (uncertainties propagated from season averaging in Fig. 6b and Monte Carlo uncertainties in Fig. B4), compared to 28±6 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, 13±3 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, 3±1 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and 2±0.2 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the unoptimized model results (uncertainties from season averaging in Fig. 6a).

The mean e ratio, defined as the growth-season mean of the 30 d lagged POC export flux divided by the growth-season mean NPP, was 0.11±0.05 (uncertainties propagated from season averaging in Fig. 6b and Monte Carlo uncertainties in Fig. B4) in the optimized model results, compared to 0.07±0.02 (uncertainties from season averaging in Fig. 6a) in the unoptimized model results. The mean f ratio, defined as the amount of NO₃ uptake divided by the amount of NO₃ and NH₄ uptake both, was 0.88±1.52 in the optimized model results, compared to 0.84±0.19 in the unoptimized model results (not shown). The higher mean f ratio relative to the mean e ratio in this study implies an imbalance between production and export at the study site, at least during the modeled period. Excess new production relative to export production (as derived from sediment traps and ²³⁴Th disequilibrium; Ducklow et al., 2018) was previously observed in the WAP, presumably due to diel vertical migration, DOM export, lateral export, and diffusive loss of PON via diapycnal mixing (Stukel et al., 2015). Stukel et al. (2015) reported up to 5 times larger new production via NO₃ uptake than export production via Th-based N export along the coastal WAP. Several additional mechanisms might be responsible for driving the discrepancy between production and export. First, given that the assimilated pool of suspended POC in the model formulation is not a good indicator of a rapidly sinking detrital pool dominating particle export, the WAP-1D-VAR v1.0 model does not capture large, short-lived particle flux events (e.g., fecal pellets produced by a large swarm of krill), underestimating POC export flux. Second, the WAP-1D-VAR v1.0 model export scheme does not consider DOC export that would lower the production–export discrepancy. Finally, RDOC is not explicitly modeled in the model, due to its much longer timescale than the model timescale, so accumulated and not-exportable RDOC pool would contribute to the deviation of the modeled e ratio from the modeled f ratio. Indeed, the modeled mean e ratios in our study, for both optimized and unoptimized cases, are situated at the lower end of the range of the e ratios measured or estimated in the WAP waters (Ducklow et al., 2018; Sailley et al., 2013; Stukel et al., 2015; Weston et al., 2013), but optimization increased the e ratio by 60 % and thus made it closer to the literature values.

The mean BP$/$NPP ratio was 0.05±0.06 (uncertainties propagated from season averaging in Fig. 6b and Monte Carlo uncertainties in Fig. B4) in the optimized model results, compared to 0.11±0.04 (uncertainties from season averaging in Fig. 6a) in the unoptimized model results. The modeled mean BP$/$NPP ratio for both optimized and unoptimized cases corresponds well to the estimates from other measurement- and observation-based studies (Ducklow et al., 2012; Kim and Ducklow, 2016). Relatively low bacterial activity in productive Antarctic waters, typically reflected as a low BP$/$PP ratio, has been attributed to low LDOM availability for bacterial growth (Kirchman et al., 2009), low temperature (Pomeroy and Wiebe, 2001), or top-down control via grazing and viral lysis (Bird and Karl, 1999).

4.4 Mean carbon stocks and flows

We summarized the growth-season means of the carbon stocks and flows in the entire food web (Fig. 7). The WAP-1D-VAR v1.0 model captured several key WAP ecological and biogeochemical features, including the lack of macronutrient limitation (NO₃ and PO₄ drawdown by phytoplankton utilization but remaining well above their half-saturation constants, Table 1) and comparable values of the assimilated and non-assimilated model state variables (Ducklow et al., 2007, 2012, 2018; Kim et al., 2016; Moline et al., 2008; Smith et al., 2008), providing confidence in the model simulations. For instance, growth-season measurements in 2017–2018 at Palmer Station showed a strongly patchy krill distribution, with the mean biomass of 0.12±0.04 mmol C m⁻³ and the maximum biomass of 0.57 mmol C m⁻³ when krill were present (unpublished data provided by D. Steinberg), falling in the range of the modeled krill biomass values (0.13±0.03 mmol C m⁻³; calculated from Fig. 7b). The WAP-1D-VAR v1.0 model also simulated several important ecosystem metrics comparable to other statistical modeling studies. For instance, the modeled phytoplankton seasonal patterns in the present study are consistent with physicochemical attributes revealed by a distinct ecological seascape pattern in the coastal WAP (Bowman et al., 2018), including low Chl and high nutrients in the first half of the growth season followed by high Chl and low nutrients in the second half of the growth season. A steady-state-solution-based inverse modeling study quantified different food-web states using ecosystem network indices from Palmer LTER annual summer cruises along the WAP shelf region (Sailley et al., 2013). Their network indices include the ratio of C export to total PP (i.e., equivalent to e ratio in our study) and the ratio of recycling (the sum of flows into respiration and DOC pool) to total PP, where more (less) recycling favorable microbial food webs are characterized by greater (smaller) ratios of recycling to total PP and smaller (greater) ratios of total C export to total PP (Legendre and Rassoulzadegan, 1996). As discussed above, the modeled mean e ratio in the present study is smaller than the estimates in the inverse modeling study for the WAP shelf region (Sailley et al., 2013) but consistent with their conclusion on the food-web status of the modeled growth season (2002–2003) positioned on the microbial food-web side. The discrepancy in the e-ratio values between the present study and Sailley et al. (2013) may be attributed to fundamentally different model formulation (i.e., time-evolving modeling for the WAP-1D-VAR v1.0 model vs. steady-state modeling) and optimization approach, or due to relatively strong microbial food-web activity at our coastal site compared to the shelf region. Microbial food-web activity can be approximated by quantifying the amount of fixed carbon flowing through heterotrophic bacteria (Carlson et al., 1999; del Giorgio and Cole, 1998; Ducklow, 2000; Ducklow et al., 2012). According to this approach, microbial food-web activity from the optimized model results was around 38 %±16 %, calculated as the ratio of bacterial L- and SDOC uptake to PP (i.e., (arrow 13 + arrow 14)$/$arrow 1 in Fig. 1, mean ± uncertainties from season averaging and Monte Carlo uncertainties in Fig. 7b). On average, SDOC supported 1 %±2 % of the total bacterial C uptake or C demand (i.e., arrow 14/(arrow 13 + arrow 14) in Fig. 1, mean ± uncertainties from season averaging and Monte Carlo uncertainties in Fig. 7b) but could be an important bacterial C source when LDOC became scarce as the growth season progressed (Fig. 5b). Indeed, several observational studies speculated that the WAP bacteria utilize SDOM in short of LDOM (Ducklow et al., 2011; Kim and Ducklow, 2016; Luria et al., 2017).

5 Summary

We developed the WAP-1D-VAR v1.0 model, a one-dimensional variational data assimilation model specific to the coastal WAP region, evaluated the model performance and robustness using a variety of quantitative metrics, and discussed the model applicability with regard to capturing the key WAP ecological and biogeochemical features using the data from the example growth season. The data assimilation scheme significantly reduced the model–observation misfits via the optimized model parameter set that adjusted the simulation to better match the Palmer LTER observations in 2002–2003. We also explored the nuanced question of how the observations influenced the data assimilation process, drove the variations in optimized parameter values relative to their corresponding initial parameter values, and affected the resulting model simulations. The WAP-1D-VAR v1.0 model successfully simulated the variables and flows not included in data assimilation, with the values consistent and comparable with other measurement- and observation-based studies in the coastal WAP. Importantly, the data assimilation scheme enabled the available observational data to constrain processes that were poorly understood, including the partitioning of NPP by different phytoplankton groups, the optimal Chl$/$C ratio of the WAP phytoplankton community, and the partitioning of DOC pools with different lability. Up to this point, a range of observational studies has provided snapshots of ecosystem and biogeochemical processes in the WAP. Yet, we have little understanding of the driving processes that underlie the connections between each component in complex food-web interactions. We used data-assimilative modeling to glue these snapshots together to better explain the observed dynamics and further understand the previously poorly constrained processes in the coastal WAP system.

Appendix A

A1 Temperature effect

$$\begin{array}{}\text{(A1)}& T_{\mathrm{f}}=\exp \left\{-A_{\mathrm{E}}\times (\mathrm{1}/T-\mathrm{1}/T_{\text{ref}})\right\}\end{array}$$

A2 Diatom processes

• Cellular quota (ratio):

  $$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,DA}}^{\mathrm{C}}=N_{\text{DA}}/C_{\text{DA}}\text{(A3)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,DA}}^{\mathrm{C}}=P_{\text{DA}}/C_{\text{DA}}\text{(A4)}& {\displaystyle }& {\displaystyle }{Q}_{\text{CHL,DA}}^{\mathrm{C}}={\text{CHL}}_{\text{DA}}/C_{\text{DA}}.\end{array}$$

• N and P limitation function:

  $$\begin{array}{ll}{\displaystyle }{N}_{\text{f,DA}}=& {\displaystyle }(Q_{\text{N,DA}}^{\mathrm{C}}-q_{\text{N,MIN,DA}}^{\mathrm{C}})/(q_{\text{N,RDF,DA}}^{\mathrm{C}}-q_{\text{N,MIN,DA}}^{\mathrm{C}})\\ \text{(A5)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le N_{\text{f,DA}}\le \mathrm{1}{\displaystyle }{P}_{\text{f,DA}}=& {\displaystyle }(Q_{\text{P,DA}}^{\mathrm{C}}-q_{\text{P,MIN,DA}}^{\mathrm{C}})/(q_{\text{P,RDF,DA}}^{\mathrm{C}}-q_{\text{P,MIN,DA}}^{\mathrm{C}})\\ \text{(A6)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le P_{\text{f,DA}}\le \mathrm{1}.\end{array}$$

• Maximum primary production rate:

  $$\begin{array}{}\text{(A7)}& P_{\text{MAX}}^{\mathrm{C}}={\mathit{\mu }}_{\text{DA}}\times T_{\mathrm{f}}\times min(N_{\text{f,DA}},P_{\text{f,DA}}).\end{array}$$

• C-specific gross primary production:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{DA}}^{\mathrm{C}}=& {\displaystyle }{C}_{\text{DA}}\times P_{\text{MAX}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }\times \mathit{\{}\mathrm{1}-\exp (-{\mathit{\alpha }}_{\text{DA}}\times Q_{\text{CHL}}^{\mathrm{C}}\times \text{PAR})/P_{\text{MAX}}^{\mathrm{C}})\mathit{\}}\\ \text{(A8)}& {\displaystyle }& {\displaystyle }\times \exp (-{\mathit{\beta }}_{\text{DA}}\times \text{PAR}).\end{array}$$

• Limitation on N and P uptake:

  $$\begin{array}{ll}{\displaystyle }{V}_{\text{MAX}}^{\mathrm{N}}=& {\displaystyle }(q_{\text{N,MAX}}^{\mathrm{C}}-Q_{\text{N,DA}}^{\mathrm{C}})/(q_{\text{N,MAX}}^{\mathrm{C}}-q_{\text{N,RDF}}^{\mathrm{C}})\\ \text{(A9)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le V_{\text{MAX}}^{\mathrm{N}}\le \mathrm{1}{\displaystyle }{V}_{\text{MAX}}^{\mathrm{P}}=& {\displaystyle }(q_{\text{P,MAX}}^{\mathrm{C}}-Q_{\text{P,DA}}^{\mathrm{C}})/(q_{\text{P,MAX}}^{\mathrm{C}}-q_{\text{P,RDF}}^{\mathrm{C}})\\ \text{(A10)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le V_{\text{MAX}}^{\mathrm{P}}\le \mathrm{1}.\end{array}$$

• N assimilation:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{DA}}^{{\mathrm{NH}}_{\mathrm{4}}}=& {\displaystyle }{C}_{\text{DA}}\times V_{\text{REF}}^{\mathrm{N}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{N}}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{NH}}_{\mathrm{4}}/({\mathrm{NH}}_{\mathrm{4}}+k^{{\mathrm{NH}}_{\mathrm{4}}}+{\mathrm{NO}}_{\mathrm{3}}\times k^{{\mathrm{NH}}_{\mathrm{4}}}/k^{{\mathrm{NO}}_{\mathrm{3}}}\left)\mathit{\right\}}{\displaystyle }{G}_{\text{DA}}^{{\mathrm{NO}}_{\mathrm{3}}}=& {\displaystyle }{C}_{\text{DA}}\times V_{\text{REF}}^{\mathrm{N}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{N}}\\ \text{(A12)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{NO}}_{\mathrm{3}}/({\mathrm{NO}}_{\mathrm{3}}+k^{{\mathrm{NO}}_{\mathrm{3}}}+{\mathrm{NH}}_{\mathrm{4}}\times k^{{\mathrm{NO}}_{\mathrm{3}}}/k^{{\mathrm{NH}}_{\mathrm{4}}}\left)\mathit{\right\}}\text{(A13)}& {\displaystyle }{G}_{\text{DA}}^{\mathrm{N}}=& {\displaystyle }{G}_{\text{DA}}^{{\mathrm{NH}}_{\mathrm{4}}}+G_{\text{DA}}^{{\mathrm{NO}}_{\mathrm{3}}}.\end{array}$$

• P assimilation:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{DA}}^{{\mathrm{PO}}_{\mathrm{4}}}=& {\displaystyle }{C}_{\text{DA}}\times V_{\text{REF}}^{\mathrm{P}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{P}}\\ \text{(A14)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{PO}}_{\mathrm{4}}/({\mathrm{PO}}_{\mathrm{4}}+k^{{\mathrm{PO}}_{\mathrm{4}}}\left)\mathit{\right\}}.\end{array}$$

• Chlorophyll production:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{DA}}^{\text{CHL}}=& {\displaystyle }\mathit{\theta }\times (G_{\text{DA}}^{{\mathrm{NH}}_{\mathrm{4}}}+G_{\text{DA}}^{{\mathrm{NO}}_{\mathrm{3}}})\\ {\displaystyle }& {\displaystyle }\times \mathit{\{}{G}_{\text{DA}}^{\mathrm{C}}/\mathit{\alpha }\times {\text{CHL}}_{\text{DA}}\times \text{PAR}\\ \text{(A15)}& {\displaystyle }& {\displaystyle }\times \exp (-\mathit{\beta }\times \text{PAR})\mathit{\}}.\end{array}$$

• Respiration:

  $$\begin{array}{}\text{(A16)}& R_{\text{DA}}^{\mathrm{C}}=G_{\text{DA}}^{{\mathrm{NO}}_{\mathrm{3}}}\times {\mathit{\zeta }}^{{\mathrm{NO}}_{\mathrm{3}}}.\end{array}$$

• Passive excretion of LDOM:

  $$\begin{array}{}\text{(A17)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDOC,PSV}}^{\mathrm{C}}={\mathit{\text{ex}}}_{\text{DA,PSV}}\times C_{\text{DA}}.\text{(A18)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDON,PSV}}^{\mathrm{N}}={\mathit{\text{ex}}}_{\text{DA,PSV}}\times N_{\text{DA}}.\text{(A19)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDOP,PSV}}^{\mathrm{P}}={\mathit{\text{ex}}}_{\text{DA,PSV}}\times P_{\text{DA}}.\end{array}$$

• Active excretion of LDOC:

  $$\begin{array}{}\text{(A20)}& E_{\text{DA,LDOC,ACT}}^{\mathrm{C}}={\mathit{\text{ex}}}_{\text{DA,ACT}}\times G_{\text{DA}}^{\mathrm{C}}.\end{array}$$

• Active excretion of SDOC:

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{DA,SDOC,ACT}}^{\mathrm{C}}=& {\displaystyle }\mathrm{0.5}\times C_{\text{DA}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{N,DA}}^{\mathrm{C}}/q_{\text{N,RDF,DA}}^{\mathrm{C}},\\ \text{(A21)}& {\displaystyle }& {\displaystyle }\mathrm{1}-Q_{\text{P,DA}}^{\mathrm{C}}/q_{\text{P,RDF,DA}}^{\mathrm{C}},\mathrm{0}).\end{array}$$

• Active excretion of SDON and SDOP (if ${\text{EX}}_{\text{DA,SDOC,ACT}}^{\mathrm{C}}>\mathrm{0}$, otherwise 0):

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{DA,SDON,ACT}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0.5}\times \mathrm{0.25}\times N_{\text{DA}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{P,DA}}^{\mathrm{N}}/(q_{\text{P,RDF,DA}}^{\mathrm{C}}\\ \text{(A22)}& {\displaystyle }& {\displaystyle }/q_{\text{N,RDF,DA}}^{\mathrm{C}}),\mathrm{0}){\displaystyle }{E}_{\text{DA,SDOP,ACT}}^{\mathrm{P}}=& {\displaystyle }\mathrm{0.5}\times \mathrm{0.25}\times P_{\text{DA}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{N,DA}}^{\mathrm{P}}/(q_{\text{N,RDF,DA}}^{\mathrm{C}}\\ \text{(A23)}& {\displaystyle }& {\displaystyle }/q_{\text{P,RDF,DA}}^{\mathrm{C}}),\mathrm{0}).\end{array}$$

• Partitioning between LDOM and SDOM:

  $$\begin{array}{}\text{(A24)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDOC}}^{\mathrm{C}}=E_{\text{DA,LDOC,PSV}}^{\mathrm{C}}+\mathrm{0.75}\times E_{\text{DA,LDOC,ACT}}^{\mathrm{C}}\text{(A25)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDON}}^{\mathrm{N}}=E_{\text{DA,LDON,PSV}}^{\mathrm{N}}\text{(A26)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,LDOP}}^{\mathrm{P}}=E_{\text{DA,LDOP,PSV}}^{\mathrm{P}}\text{(A27)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,SDOC}}^{\mathrm{C}}=E_{\text{DA,SDOC,ACT}}^{\mathrm{C}}+\mathrm{0.25}\times E_{\text{DA,LDOC,ACT}}^{\mathrm{C}}\text{(A28)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,SDON}}^{\mathrm{N}}=E_{\text{DA,SDON,ACT}}^{\mathrm{N}}\text{(A29)}& {\displaystyle }& {\displaystyle }{E}_{\text{DA,SDOP}}^{\mathrm{P}}=E_{\text{DA,SDOP,ACT}}^{\mathrm{P}}.\end{array}$$

• POM production by aggregation:

  $$\begin{array}{}\text{(A30)}& {\displaystyle }& {\displaystyle }{D}_{\text{DA}}^{\mathrm{C}}={\text{pom}}_{\text{DA}}\times C_{\text{DA}}\times C_{\text{DA}}\text{(A31)}& {\displaystyle }& {\displaystyle }{D}_{\text{DA}}^{\mathrm{N}}=Q_{\text{N,DA}}^{\mathrm{C}}\times D_{\text{DA}}^{\mathrm{C}}\text{(A32)}& {\displaystyle }& {\displaystyle }{D}_{\text{DA}}^{\mathrm{P}}=Q_{\text{P,DA}}^{\mathrm{C}}\times D_{\text{DA}}^{\mathrm{C}}\text{(A33)}& {\displaystyle }& {\displaystyle }{D}_{\text{DA}}^{\text{CHL}}=Q_{\text{CHL,DA}}^{\mathrm{C}}\times D_{\text{DA}}^{\mathrm{C}}.\end{array}$$

• Grazing by microzooplankton:

  $$\begin{array}{ll}{\displaystyle }G{Z}_{\text{DA,MZ}}^{\mathrm{C}}=& {\displaystyle }{T}_{\mathrm{f}}\times {\mathit{\mu }}_{\text{MZ}}\times C_{\text{MZ}}\\ {\displaystyle }& {\displaystyle }\times \left[\right(C_{\text{DA}}-{\mathit{ϵ}}_{\text{DA}}{)}^{\mathrm{2}}/\mathit{\left\{}\right(C_{\text{DA}}-{\mathit{ϵ}}_{\text{DA}}{)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+g_{\text{DA}}^{\mathrm{2}}+(C_{\text{CRYPTO}}\times g_{\text{DA}}/g_{\text{CR}}{)}^{\mathrm{2}}\\ \text{(A34)}& {\displaystyle }& {\displaystyle }+(C_{\text{BAC}}\times g_{\text{DA}}/g_{\text{BAC}}{)}^{\mathrm{2}}\mathit{\}}]\text{(A35)}& {\displaystyle }G{Z}_{\text{DA,MZ}}^{\mathrm{N}}=& {\displaystyle }{Q}_{\text{N,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,MZ}}^{\mathrm{C}}\text{(A36)}& {\displaystyle }G{Z}_{\text{DA,MZ}}^{\mathrm{P}}=& {\displaystyle }{Q}_{\text{P,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,MZ}}^{\mathrm{C}}\text{(A37)}& {\displaystyle }G{Z}_{\text{DA,MZ}}^{\text{CHL}}=& {\displaystyle }{Q}_{\text{CHL,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,MZ}}^{\mathrm{C}}.\end{array}$$

• Grazing by krill:

  $$\begin{array}{ll}{\displaystyle }G{Z}_{\text{DA,KR}}^{\mathrm{C}}=& {\displaystyle }{T}_{\mathrm{f}}\times {\mathit{\mu }}_{\text{KR}}\times C_{\text{KR}}\\ {\displaystyle }& {\displaystyle }\times [C_{\text{DA}}^{\mathrm{2}}/\mathit{\{}{C}_{\text{DA}}^{\mathrm{2}}+g_{{\text{DA}}^{\mathrm{2}}}^{\prime }\\ \text{(A38)}& {\displaystyle }& {\displaystyle }+(C_{\text{MZ}}\times g_{\text{DA}}^{\prime }/g_{\text{MZ}}{)}^{\mathrm{2}}\mathit{\}}]\text{(A39)}& {\displaystyle }G{Z}_{\text{DA,KR}}^{\mathrm{N}}=& {\displaystyle }{Q}_{\text{N,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,KR}}^{\mathrm{C}}\text{(A40)}& {\displaystyle }G{Z}_{\text{DA,KR}}^{\mathrm{P}}=& {\displaystyle }{Q}_{\text{P,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,KR}}^{\mathrm{C}}\text{(A41)}& {\displaystyle }G{Z}_{\text{DA,KR}}^{\text{CHL}}=& {\displaystyle }{Q}_{\text{CHL,DA}}^{\mathrm{C}}\times G{Z}_{\text{DA,KR}}^{\mathrm{C}}.\end{array}$$

• The net growth rate equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{DA}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{DA}}^{\mathrm{C}}-E_{\text{DA,LDOC}}^{\mathrm{C}}-E_{\text{DA,SDOC}}^{\mathrm{C}}-D_{\text{DA}}^{\mathrm{C}}\\ \text{(A42)}& {\displaystyle }& {\displaystyle }-R_{\text{DA}}^{\mathrm{C}}-G{Z}_{\text{DA,MZ}}^{\mathrm{C}}-G{Z}_{\text{DA,KR}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{DA}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{DA}}^{\mathrm{N}}-E_{\text{DA,LDON}}^{\mathrm{N}}-E_{\text{DA,SDON}}^{\mathrm{N}}-D_{\text{DA}}^{\mathrm{N}}\\ \text{(A43)}& {\displaystyle }& {\displaystyle }-G{Z}_{\text{DA,MZ}}^{\mathrm{N}}-G{Z}_{\text{DA,KR}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{DA}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{DA}}^{\mathrm{P}}-E_{\text{DA,LDOP}}^{\mathrm{P}}-E_{\text{DA,SDOP}}^{\mathrm{P}}-D_{\text{DA}}^{\mathrm{N}}\\ \text{(A44)}& {\displaystyle }& {\displaystyle }-G{Z}_{\text{DA,MZ}}^{\mathrm{P}}-G{Z}_{\text{DA,KR}}^{\mathrm{P}}\text{(A45)}& {\displaystyle \frac{\mathrm{d}{\text{CHL}}_{\text{DA}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{DA}}^{\text{CHL}}-D_{\text{DA}}^{\text{CHL}}-G{Z}_{\text{DA,MZ}}^{\text{CHL}}-G{Z}_{\text{DA,KR}}^{\text{CHL}}.\end{array}$$

A3 Cryptophyte processes

• Cellular quota (ratio):

  $$\begin{array}{}\text{(A46)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,CR}}^{\mathrm{C}}=N_{\text{CR}}/C_{\text{CR}}\text{(A47)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,CR}}^{\mathrm{C}}=P_{\text{CR}}/C_{\text{CR}}\text{(A48)}& {\displaystyle }& {\displaystyle }{Q}_{\text{CHL,CR}}^{\mathrm{C}}={\text{CHL}}_{\text{CR}}/C_{\text{CR}}.\end{array}$$

• N and P limitation function:

  $$\begin{array}{ll}{\displaystyle }{N}_{\text{f,CR}}=& {\displaystyle }(Q_{\text{N,CR}}^{\mathrm{C}}-q_{\text{N,MIN,CR}}^{\mathrm{C}})/(q_{\text{N,RDF,CR}}^{\mathrm{C}}-q_{\text{N,MIN,CR}}^{\mathrm{C}})\\ \text{(A49)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le N_{\text{f,CR}}\le \mathrm{1}{\displaystyle }{P}_{\text{f,CR}}=& {\displaystyle }(Q_{\text{P,CR}}^{\mathrm{C}}-q_{\text{P,MIN,CR}}^{\mathrm{C}})/(q_{\text{P,RDF,CR}}^{\mathrm{C}}-q_{\text{P,MIN,CR}}^{\mathrm{C}})\\ \text{(A50)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le P_{\text{f,CR}}\le \mathrm{1}.\end{array}$$

• Maximum primary production rate:

  $$\begin{array}{}\text{(A51)}& P_{\text{MAX}}^{\mathrm{C}}={\mathit{\mu }}_{\text{CR}}\times T_{\mathrm{f}}\times min(N_{\text{f,CR}},P_{\text{f,CR}}).\end{array}$$

• C-specific gross primary production:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{CR}}^{\mathrm{C}}=& {\displaystyle }{C}_{\text{CR}}\times P_{\text{MAX}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }\times \mathit{\{}\mathrm{1}-\exp (-{\mathit{\alpha }}_{\text{CR}}\times Q_{\text{CHL}}^{\mathrm{C}}\times \text{PAR})/P_{\text{MAX}}^{\mathrm{C}})\mathit{\}}\\ \text{(A52)}& {\displaystyle }& {\displaystyle }\times \exp (-{\mathit{\beta }}_{\text{CR}}\times \text{PAR}).\end{array}$$

• Limitation on N and P uptake:

  $$\begin{array}{ll}{\displaystyle }{V}_{\text{MAX}}^{\mathrm{N}}=& {\displaystyle }(q_{\text{N,MAX}}^{\mathrm{C}}-Q_{\text{N,CR}}^{\mathrm{C}})/(q_{\text{N,MAX}}^{\mathrm{C}}-q_{\text{N,RDF}}^{\mathrm{C}})\\ \text{(A53)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le V_{\text{MAX}}^{\mathrm{N}}\le \mathrm{1}{\displaystyle }{V}_{\text{MAX}}^{\mathrm{P}}=& {\displaystyle }(q_{\text{P,MAX}}^{\mathrm{C}}-Q_{\text{P,CR}}^{\mathrm{C}})/(q_{\text{P,MAX}}^{\mathrm{C}}-q_{\text{P,RDF}}^{\mathrm{C}})\\ \text{(A54)}& {\displaystyle }& {\displaystyle }\mathrm{0}\le V_{\text{MAX}}^{\mathrm{P}}\le \mathrm{1}.\end{array}$$

• Nitrogen assimilation:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{CR}}^{{\mathrm{NH}}_{\mathrm{4}}}=& {\displaystyle }{C}_{\text{CR}}\times V_{\text{REF}}^{\mathrm{N}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{N}}\\ \text{(A55)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{NH}}_{\mathrm{4}}/({\mathrm{NH}}_{\mathrm{4}}+k^{{\mathrm{NH}}_{\mathrm{4}}}+{\mathrm{NO}}_{\mathrm{3}}\times k^{{\mathrm{NH}}_{\mathrm{4}}}/k^{{\mathrm{NO}}_{\mathrm{3}}}\left)\mathit{\right\}}{\displaystyle }{G}_{\text{CR}}^{{\mathrm{NO}}_{\mathrm{3}}}=& {\displaystyle }{C}_{\text{CR}}\times V_{\text{REF}}^{\mathrm{N}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{N}}\\ \text{(A56)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{NO}}_{\mathrm{3}}/({\mathrm{NO}}_{\mathrm{3}}+k^{{\mathrm{NO}}_{\mathrm{3}}}+{\mathrm{NH}}_{\mathrm{4}}\times k^{{\mathrm{NO}}_{\mathrm{3}}}/k^{{\mathrm{NH}}_{\mathrm{4}}}\left)\mathit{\right\}}\text{(A57)}& {\displaystyle }{G}_{\text{CR}}^{\mathrm{N}}=& {\displaystyle }{G}_{\text{CR}}^{{\mathrm{NH}}_{\mathrm{4}}}+G_{\text{CR}}^{{\mathrm{NO}}_{\mathrm{3}}}.\end{array}$$

• Phosphorus assimilation:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{CR}}^{{\mathrm{PO}}_{\mathrm{4}}}=& {\displaystyle }{C}_{\text{CR}}\times V_{\text{REF}}^{\mathrm{P}}\times T_{\mathrm{f}}\times V_{\text{MAX}}^{\mathrm{P}}\\ \text{(A58)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathrm{PO}}_{\mathrm{4}}/({\mathrm{PO}}_{\mathrm{4}}+k^{{\mathrm{PO}}_{\mathrm{4}}}\left)\mathit{\right\}}.\end{array}$$

• Chlorophyll production:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{CR}}^{\text{CHL}}=& {\displaystyle }\mathit{\theta }\times (G_{\text{CR}}^{{\mathrm{NH}}_{\mathrm{4}}}+G_{\text{CR}}^{{\mathrm{NO}}_{\mathrm{3}}})\\ {\displaystyle }& {\displaystyle }\times \mathit{\{}{G}_{\text{CR}}^{\mathrm{C}}/\mathit{\alpha }\times {\text{CHL}}_{\text{CR}}\times \text{PAR}\\ \text{(A59)}& {\displaystyle }& {\displaystyle }\times \exp (-\mathit{\beta }\times \text{PAR})\mathit{\}}.\end{array}$$

• Respiration:

  $$\begin{array}{}\text{(A60)}& R_{\text{CR}}^{\mathrm{C}}=G_{\text{CR}}^{{\mathrm{NO}}_{\mathrm{3}}}\times {\mathit{\zeta }}^{{\mathrm{NO}}_{\mathrm{3}}}.\end{array}$$

• Passive excretion of LDOM:

  $$\begin{array}{}\text{(A61)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDOC,PSV}}^{\mathrm{C}}={\mathit{\text{ex}}}_{\text{CR,PSV}}\times C_{\text{CR}}\text{(A62)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDON,PSV}}^{\mathrm{N}}={\mathit{\text{ex}}}_{\text{CR,PSV}}\times N_{\text{CR}}\text{(A63)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDOP,PSV}}^{\mathrm{P}}={\mathit{\text{ex}}}_{\text{CR,PSV}}\times P_{\text{CR}}.\end{array}$$

• Active excretion of LDOC:

  $$\begin{array}{}\text{(A64)}& E_{\text{CR,LDOC,ACT}}^{\mathrm{C}}={\mathit{\text{ex}}}_{\text{CR,ACT}}\times G_{\text{CR}}^{\mathrm{C}}.\end{array}$$

• Active excretion of SDOC:

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{CR,SDOC,ACT}}^{\mathrm{C}}=& {\displaystyle }\mathrm{0.5}\times C_{\text{CR}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{N,CR}}^{\mathrm{C}}/q_{\text{N,RDF,CR}}^{\mathrm{C}},\\ \text{(A65)}& {\displaystyle }& {\displaystyle }\mathrm{1}-Q_{\text{P,CR}}^{\mathrm{C}}/q_{\text{P,RDF,CR}}^{\mathrm{C}},\mathrm{0}).\end{array}$$

• Active excretion of SDON and SDOP (if ${\text{EX}}_{\text{CR,SDOC,ACT}}^{\mathrm{C}}>\mathrm{0}$, otherwise 0):

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{CR,SDON,ACT}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0.5}\times \mathrm{0.25}\times N_{\text{CR}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{P,CR}}^{\mathrm{N}}/(q_{\text{P,RDF,CR}}^{\mathrm{C}}\\ \text{(A66)}& {\displaystyle }& {\displaystyle }/q_{\text{N,RDF,CR}}^{\mathrm{C}}),\mathrm{0}){\displaystyle }{E}_{\text{CR,SDOP,ACT}}^{\mathrm{P}}=& {\displaystyle }\mathrm{0.5}\times \mathrm{0.25}\times P_{\text{CR}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{1}-Q_{\text{N,CR}}^{\mathrm{P}}/(q_{\text{N,RDF,CR}}^{\mathrm{C}}\\ \text{(A67)}& {\displaystyle }& {\displaystyle }/q_{\text{P,RDF,CR}}^{\mathrm{C}}),\mathrm{0}).\end{array}$$

• Partitioning between LDOM and SDOM:

  $$\begin{array}{}\text{(A68)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDOC}}^{\mathrm{C}}=E_{\text{CR,LDOC,PSV}}^{\mathrm{C}}+\mathrm{0.75}\times E_{\text{CR,LDOC,ACT}}^{\mathrm{C}}\text{(A69)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDON}}^{\mathrm{N}}=E_{\text{CR,LDON,PSV}}^{\mathrm{N}}\text{(A70)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,LDOP}}^{\mathrm{P}}=E_{\text{CR,LDOP,PSV}}^{\mathrm{P}}\text{(A71)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,SDOC}}^{\mathrm{C}}=E_{\text{CR,SDOC,ACT}}^{\mathrm{C}}+\mathrm{0.25}\times E_{\text{CR,LDOC,ACT}}^{\mathrm{C}}\text{(A72)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,SDON}}^{\mathrm{N}}=E_{\text{CR,SDON,ACT}}^{\mathrm{N}}\text{(A73)}& {\displaystyle }& {\displaystyle }{E}_{\text{CR,SDOP}}^{\mathrm{P}}=E_{\text{CR,SDOP,ACT}}^{\mathrm{P}}.\end{array}$$

• POM production by aggregation:

  $$\begin{array}{}\text{(A74)}& {\displaystyle }& {\displaystyle }{D}_{\text{CR}}^{\mathrm{C}}={\text{pom}}_{\text{CR}}\times C_{\text{CR}}\times C_{\text{CR}}\text{(A75)}& {\displaystyle }& {\displaystyle }{D}_{\text{CR}}^{\mathrm{N}}=Q_{\text{N,CR}}^{\mathrm{C}}\times A_{\text{CR}}^{\mathrm{C}}\text{(A76)}& {\displaystyle }& {\displaystyle }{D}_{\text{CR}}^{\mathrm{P}}=Q_{\text{P,CR}}^{\mathrm{C}}\times A_{\text{CR}}^{\mathrm{C}}\text{(A77)}& {\displaystyle }& {\displaystyle }{D}_{\text{CR}}^{\text{CHL}}=Q_{\text{CHL,CR}}^{\mathrm{C}}\times A_{\text{CR}}^{\mathrm{C}}.\end{array}$$

• Grazing by microzooplankton:

  $$\begin{array}{ll}{\displaystyle }G{Z}_{\text{CR}}^{\mathrm{C}}=& {\displaystyle }{T}_{\text{f}}\times {\mathit{\mu }}_{\text{MZ}}\times C_{\text{MZ}}\\ {\displaystyle }& {\displaystyle }\times \left[\right(C_{\text{CR}}-{\mathit{ϵ}}_{\text{CR}}{)}^{\mathrm{2}}/\mathit{\left\{}\right(C_{\text{CR}}-{\mathit{ϵ}}_{\text{CR}}{)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+g_{\text{CR}}^{\mathrm{2}}+(C_{\text{DA}}\times g_{\text{CR}}/g_{\text{DA}}{)}^{\mathrm{2}}\\ \text{(A78)}& {\displaystyle }& {\displaystyle }+(C_{\text{BAC}}\times g_{\text{CR}}/g_{\text{BAC}}{)}^{\mathrm{2}}\mathit{\}}]\text{(A79)}& {\displaystyle }G{Z}_{\text{CR}}^{\mathrm{N}}=& {\displaystyle }{Q}_{\text{N,CR}}^{\mathrm{C}}\times G{Z}_{\text{CR,MZ}}^{\mathrm{C}}\text{(A80)}& {\displaystyle }G{Z}_{\text{CR}}^{\mathrm{P}}=& {\displaystyle }{Q}_{\text{P,CR}}^{\mathrm{C}}\times G{Z}_{\text{CR,MZ}}^{\mathrm{C}}\text{(A81)}& {\displaystyle }G{Z}_{\text{CR}}^{\text{CHL}}=& {\displaystyle }{Q}_{\text{CHL,CR}}^{\mathrm{C}}\times G{Z}_{\text{CR,MZ}}^{\mathrm{C}}.\end{array}$$

• The net growth rate equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{CR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{CR}}^{\mathrm{C}}-E_{\text{CR,LDOC}}^{\mathrm{C}}-E_{\text{CR,SDOC}}^{\mathrm{C}}-D_{\text{CR}}^{\mathrm{C}}\\ \text{(A82)}& {\displaystyle }& {\displaystyle }-R_{\text{CR}}^{\mathrm{C}}-G{Z}_{\text{CR}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{CR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{CR}}^{\mathrm{N}}-E_{\text{CR,LDON}}^{\mathrm{N}}-E_{\text{CR,SDON}}^{\mathrm{N}}-D_{\text{CR}}^{\mathrm{N}}\\ \text{(A83)}& {\displaystyle }& {\displaystyle }-G{Z}_{\text{CR}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{CR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{CR}}^{\mathrm{P}}-E_{\text{CR,LDOP}}^{\mathrm{P}}-E_{\text{CR,SDOP}}^{\mathrm{P}}-D_{\text{CR}}^{\mathrm{N}}\\ \text{(A84)}& {\displaystyle }& {\displaystyle }-G{Z}_{\text{CR}}^{\mathrm{P}}\text{(A85)}& {\displaystyle \frac{\mathrm{d}{\text{CHL}}_{\text{CR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{CR}}^{\text{CHL}}-D_{\text{CR}}^{\text{CHL}}-G{Z}_{\text{CR}}^{\text{CHL}}.\end{array}$$

A4 Bacterial processes

• Cellular quota (ratio):

  $$\begin{array}{}\text{(A86)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,BAC}}^{\mathrm{C}}=N_{\text{BAC}}/C_{\text{BAC}}\text{(A87)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,BAC}}^{\mathrm{C}}=P_{\text{BAC}}/C_{\text{BAC}}\text{(A88)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,BAC}}^{\mathrm{P}}=N_{\text{BAC}}/P_{\text{BAC}}\text{(A89)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,LDOM}}^{\mathrm{C}}=N_{\text{LDOM}}/C_{\text{LDOM}}\text{(A90)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,LDOM}}^{\mathrm{C}}=P_{\text{LDOM}}/C_{\text{LDOM}}\text{(A91)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,SDOM}}^{\mathrm{C}}=N_{\text{SDOM}}/C_{\text{SDOM}}\text{(A92)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,SDOM}}^{\mathrm{C}}=P_{\text{SDOM}}/C_{\text{SDOM}}.\end{array}$$

• N and P limitation function:

  $$\begin{array}{}\text{(A93)}& {\displaystyle }& {\displaystyle }{N}_{\text{f,BAC}}=Q_{\text{N,BAC}}^{\mathrm{C}}/q_{\text{N,BAC}}^{\mathrm{C}}& {\displaystyle }\mathrm{0}\le N_{\text{f,BAC}}\le \mathrm{1}\text{(A94)}& {\displaystyle }& {\displaystyle }{P}_{\text{f,BAC}}=Q_{\text{P,BAC}}^{\mathrm{C}}/q_{\text{P,BAC}}^{\mathrm{C}}& {\displaystyle }\mathrm{0}\le P_{\text{f,BAC}}\le \mathrm{1}.\end{array}$$

• Maximum available LDOC and SDOC:

  $$\begin{array}{}\text{(A95)}& {\displaystyle }& {\displaystyle }\text{ALC}=C_{\text{LDOC}}\text{(A96)}& {\displaystyle }& {\displaystyle }\text{ASC}=r_{\text{SDOC}}\times C_{\text{SDOC}}.\end{array}$$

• Bacterial uptake of LDOC and SDOC (i.e., bacterial gross C growth):

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{BAC,LDOC}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\mu }}_{\text{BAC}}\times T_{\mathrm{f}}\times C_{\text{BAC}}\\ {\displaystyle }& {\displaystyle }\times min(N_{\text{f,BAC}},P_{\text{f,BAC}})\\ \text{(A97)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}\text{ALC}/(\text{ALC}+k^{\text{DOC}}+\text{ASC}\left)\mathit{\right\}}{\displaystyle }{G}_{\text{BAC,SDOC}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\mu }}_{\text{BAC}}\times T_{\mathrm{f}}\times C_{\text{BAC}}\\ {\displaystyle }& {\displaystyle }\times min(N_{\text{f,BAC}},P_{\text{f,BAC}})\\ \text{(A98)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}\text{ASC}/(\text{ASC}+k^{\text{DOC}}+\text{ALC}\left)\mathit{\right\}}\text{(A99)}& {\displaystyle }{G}_{\text{BAC,DOC}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDOC}}^{\mathrm{C}}+G_{\text{BAC,SDOC}}^{\mathrm{C}}.\end{array}$$

• Bacterial N uptake:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{BAC,LDON}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDOC}}^{\mathrm{C}}\times Q_{\text{N,LDOM}}^{\mathrm{C}}\\ \text{(A100)}& {\displaystyle }{G}_{\text{BAC,SDON}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,SDOC}}^{\mathrm{C}}{\displaystyle }& {\displaystyle }\times min\mathit{\{}{q}_{\text{N,BAC}}^{\mathrm{C}},Q_{\text{N,SDOM}}^{\mathrm{C}}+f_{\text{S}}/N_{\text{f,BAC}}\\ \text{(A101)}& {\displaystyle }& {\displaystyle }\times (q_{\text{N,BAC}}^{\mathrm{C}}-Q_{\text{N,SDOM}}^{\mathrm{C}})\mathit{\}}{\displaystyle }{G}_{\text{BAC},{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDON}}^{\mathrm{C}}\\ \text{(A102)}& {\displaystyle }& {\displaystyle }\times {\mathrm{NH}}_{\mathrm{4}}/N_{\text{LDOM}}/min(\mathrm{1},N_{\text{f,BAC}}),\end{array}$$

  if N_f,BAC<1,

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{BAC,}{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{C}}=& {\displaystyle }min\mathit{\{}\mathrm{0.1}\times {\mathrm{NO}}_{\mathrm{3}}\times \mathrm{1}/min(\mathrm{1},N_{\mathrm{f}}{,}_{\text{BAC}})\\ {\displaystyle }& {\displaystyle }\times (G_{\text{BAC,LDON}}^{\mathrm{C}}+G_{\text{BAC,SDON}}^{\mathrm{C}})\\ {\displaystyle }& {\displaystyle }/(N_{\text{LDOM}}+N_{\text{SDOM}}),({\mathrm{NO}}_{\mathrm{3}}+{\mathrm{NH}}_{\mathrm{4}})\\ {\displaystyle }& {\displaystyle }\times (G_{\text{BAC,LDON}}^{\mathrm{N}}+G_{\text{BAC,SDON}}^{\mathrm{N}})\\ \text{(A103)}& {\displaystyle }& {\displaystyle }/(N_{\text{LDOM}}+N_{\text{SDOM}})-G_{\text{BAC}}^{{\mathrm{NH}}_{\mathrm{4}}}\mathit{\}},\end{array}$$

  else,

  $$\begin{array}{ll}\text{(A104)}& {\displaystyle }{G}_{\text{BAC,}{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{C}}=& {\displaystyle }\mathrm{0}{\displaystyle }{G}_{\text{BAC,N}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDON}}^{\mathrm{C}}+G_{\text{BAC,SDON}}^{\mathrm{C}}+G_{\text{BAC},{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{C}}\\ \text{(A105)}& {\displaystyle }& {\displaystyle }+G_{\text{BAC},{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{C}}.\end{array}$$

• Bacterial P uptake:

  $$\begin{array}{ll}{\displaystyle }{G}_{\text{BAC,LDOP}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDOC}}^{\mathrm{C}}\times Q_{\text{P,LDOM}}^{\mathrm{C}}\\ \text{(A106)}& {\displaystyle }{G}_{\text{BAC,SDOP}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,SDOC}}^{\mathrm{C}}{\displaystyle }& {\displaystyle }\times min\mathit{\{}{q}_{\text{P,BAC}}^{\mathrm{C}},Q_{\text{P,SDOM}}^{\mathrm{C}}+f_{\mathrm{S}}/P_{\text{f,BAC}}\\ \text{(A107)}& {\displaystyle }& {\displaystyle }\times (q_{\text{P,BAC}}^{\mathrm{C}}-Q_{\text{P,SDOM}}^{\mathrm{C}})\mathit{\}}{\displaystyle }{G}_{\text{BAC,}{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDON}}^{\mathrm{C}}\\ \text{(A108)}& {\displaystyle }& {\displaystyle }\times {\mathrm{PO}}_{\mathrm{4}}/P_{\text{LDOM}}/min(\mathrm{1},P_{\text{f,BAC}})\text{(A109)}& {\displaystyle }{G}_{\text{BAC,P}}^{\mathrm{C}}=& {\displaystyle }{G}_{\text{BAC,LDOP}}^{\mathrm{C}}+G_{\text{BAC,SDOP}}^{\mathrm{C}}+G_{\text{BAC,}{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{C}}.\end{array}$$

• Respiration:

  $$\begin{array}{ll}{\displaystyle }{R}_{\text{BAC}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\zeta }}^{{\mathrm{NO}}_{\mathrm{3}}}\times G_{\text{BAC},{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{C}}+r_{\text{BAC}}^{\mathrm{B}}\times T_{\mathrm{f}}\times C_{\text{BAC}}\\ {\displaystyle }& {\displaystyle }+\mathit{\{}{r}_{\text{min,BAC}}^{\mathrm{A}}+(r_{\text{max,BAC}}^{\mathrm{A}}-r_{\text{min,BAC}}^{\mathrm{A}})\\ {\displaystyle }& {\displaystyle }\times \exp (-b_{\text{R,BAC}}\times G_{\text{BAC,DOC}}^{\mathrm{C}})\mathit{\}}\\ \text{(A110)}& {\displaystyle }& {\displaystyle }\times G_{\text{BAC,DOC}}^{\mathrm{C}}.\end{array}$$

• RDOC release:

  $$\begin{array}{}\text{(A111)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,RDOC}}^{\mathrm{C}}={\text{refr}}_{\text{BAC}}\times C_{\text{BAC}}\text{(A112)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,RDON}}^{\mathrm{N}}=E_{\text{BAC,RDOC}}^{\mathrm{C}}\times q_{\text{N,RDOM}}^{\mathrm{C}}\text{(A113)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,RDOP}}^{\mathrm{P}}=E_{\text{BAC,RDOC}}^{\mathrm{C}}\times q_{\text{P,RDOM}}^{\mathrm{C}}.\end{array}$$

• Remineralization of inorganic nutrients:
  if $Q_{\text{N,BAC}}^{\mathrm{C}}>q_{\text{N,BAC}}^{\mathrm{C}}$ and $Q_{\text{P,BAC}}^{\mathrm{C}}>q_{\text{P,BAC}}^{\mathrm{C}}$ (i.e., C in short)

  $$\begin{array}{}\text{(A114)}& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{N}}=& {\displaystyle }{\text{remi}}_{\text{BAC}}\times (N_{\text{BAC}}-C_{\text{BAC}}\times q_{\text{N,BAC}}^{\mathrm{C}})\text{(A115)}& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{P}}=& {\displaystyle }{\text{remi}}_{\text{BAC}}\times (P_{\text{BAC}}-C_{\text{BAC}}\times q_{\text{P,BAC}}^{\mathrm{C}})\end{array}$$

  elseif $Q_{\text{N,BAC}}^{\mathrm{C}}<q_{\text{N,BAC}}^{\mathrm{C}}$ and $Q_{\text{N,BAC}}^{\mathrm{P}}<q_{\text{N,BAC}}^{\mathrm{C}}/q_{\text{P,BAC}}^{\mathrm{C}}$ (i.e., N in short)

  $$\begin{array}{}\text{(A116)}& {\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{N}}=\mathrm{0}\text{(A117)}& {\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{P}}=\mathrm{0},\end{array}$$

  else (i.e., P in short)

  $$\begin{array}{}\text{(A118)}& {\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{N}}=\mathrm{0}\text{(A119)}& {\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{BAC}}^{\mathrm{P}}=\mathrm{0}.\end{array}$$

• SDOM excretion to adjust stoichiometry:
  if $Q_{\text{N,BAC}}^{\mathrm{C}}>q_{\text{N,BAC}}^{\mathrm{C}}$ and $Q_{\text{P,BAC}}^{\mathrm{C}}>q_{\text{P,BAC}}^{\mathrm{C}}$ (i.e., C in short)

  $$\begin{array}{}\text{(A120)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,SDOC}}^{\mathrm{C}}=\mathrm{0}\text{(A121)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,SDON}}^{\mathrm{N}}=\mathrm{0}\text{(A122)}& {\displaystyle }& {\displaystyle }{E}_{\text{BAC,SDOP}}^{\mathrm{P}}=\mathrm{0},\end{array}$$

  else if $Q_{\text{N,BAC}}^{\mathrm{C}}<q_{\text{N,BAC}}^{\mathrm{C}}$ and $Q_{\text{N,BAC}}^{\mathrm{P}}<q_{\text{N,BAC}}^{\mathrm{C}}/q_{\text{P,BAC}}^{\mathrm{C}}$ (i.e., N in short)

  $$\begin{array}{ll}\text{(A123)}& {\displaystyle }{E}_{\text{BAC,SDOC}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,BAC}}\times (C_{\text{BAC}}-N_{\text{BAC}}/q_{\text{N,BAC}}^{\mathrm{C}})\text{(A124)}& {\displaystyle }{E}_{\text{BAC,SDOC}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0}{\displaystyle }{E}_{\text{BAC,SDOP}}^{\mathrm{P}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,BAC}}\\ \text{(A125)}& {\displaystyle }& {\displaystyle }\times (P_{\text{BAC}}-N_{\text{BAC}}/q_{\text{N,BAC}}^{\mathrm{C}}\times q_{\text{P,BAC}}^{\mathrm{C}}),\end{array}$$

  else (i.e., P in short)

  $$\begin{array}{ll}\text{(A126)}& {\displaystyle }{E}_{\text{BAC,SDOC}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,BAC}}\times (C_{\text{BAC}}-P_{\text{BAC}}/q_{\text{P,BAC}}^{\mathrm{C}}){\displaystyle }{E}_{\text{BAC,SDON}}^{\mathrm{N}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,BAC}}\\ \text{(A127)}& {\displaystyle }& {\displaystyle }\times (N_{\text{BAC}}-P_{\text{BAC}}/q_{\text{P,BAC}}^{\mathrm{C}}\times q_{\text{N,BAC}}^{\mathrm{C}})\text{(A128)}& {\displaystyle }{E}_{\text{BAC,SDOP}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0}.\end{array}$$

• Grazing by microzooplankton:

  $$\begin{array}{ll}{\displaystyle }G{Z}_{\text{BAC}}^{\mathrm{C}}=& {\displaystyle }{T}_{\mathrm{f}}\times {\mathit{\mu }}_{\text{MZ}}\times C_{\text{MZ}}\\ {\displaystyle }& {\displaystyle }\times [C_{\text{BAC}}^{\mathrm{2}}/\mathit{\{}{C}_{\text{CR}}^{\mathrm{2}}+g_{\text{BAC}}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+(C_{\text{DA}}\times g_{\text{BAC}}/g_{\text{DA}}{)}^{\mathrm{2}}\\ \text{(A129)}& {\displaystyle }& {\displaystyle }+(C_{\text{CR}}\times g_{\text{BAC}}/g_{\text{CR}}{)}^{\mathrm{2}}\mathit{\}}]\text{(A130)}& {\displaystyle }G{Z}_{\text{BAC}}^{\mathrm{N}}=& {\displaystyle }G{Z}_{\text{BAC}}^{\mathrm{C}}\times Q_{\text{N,BAC}}^{\mathrm{C}}\text{(A131)}& {\displaystyle }G{Z}_{\text{BAC}}^{\mathrm{P}}=& {\displaystyle }G{Z}_{\text{BAC}}^{\mathrm{C}}\times Q_{\text{P,BAC}}^{\mathrm{C}}.\end{array}$$

• Viral mortality:

  $$\begin{array}{}\text{(A132)}& {\displaystyle }& {\displaystyle }{M}_{\text{BAC}}^{\mathrm{C}}=m_{\text{BAC}}\times C_{\text{BAC}}\text{(A133)}& {\displaystyle }& {\displaystyle }{M}_{\text{BAC}}^{\mathrm{N}}=m_{\text{BAC}}\times N_{\text{BAC}}\text{(A134)}& {\displaystyle }& {\displaystyle }{M}_{\text{BAC}}^{\mathrm{P}}=m_{\text{BAC}}\times P_{\text{BAC}}.\end{array}$$

• Net flux of inorganic nutrients through bacteria:

  $$\begin{array}{}\text{(A135)}& {\displaystyle }& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{NH}}_{\mathrm{4}}}={\text{REMI}}_{\text{BAC}}^{\mathrm{N}}-G_{\text{BAC,}{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{C}}\text{(A136)}& {\displaystyle }& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{NO}}_{\mathrm{3}}}=-G_{\text{BAC,}{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{C}}\text{(A137)}& {\displaystyle }& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{PO}}_{\mathrm{4}}}={\text{REMI}}_{\text{BAC}}^{\mathrm{P}}-G_{\text{BAC,}{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{C}}.\end{array}$$

• The net growth rate equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{BAC}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{BAC,DOC}}^{\mathrm{C}}-E_{\text{BAC,SDOC}}^{\mathrm{C}}-E_{\text{BAC,RDOC}}^{\mathrm{C}}\\ \text{(A138)}& {\displaystyle }& {\displaystyle }-R_{\text{BAC}}^{\mathrm{C}}-G{Z}_{\text{BAC}}^{\mathrm{C}}-M_{\text{BAC}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{BAC}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{BAC,DON}}^{\mathrm{N}}-E_{\text{BAC,SDON}}^{\mathrm{N}}-E_{\text{BAC,RDON}}^{\mathrm{N}}\\ \text{(A139)}& {\displaystyle }& {\displaystyle }-R_{\text{BAC}}^{\mathrm{N}}-G{Z}_{\text{BAC}}^{\mathrm{N}}-M_{\text{BAC}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{BAC}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{BAC,DOP}}^{\mathrm{P}}-E_{\text{BAC,SDOP}}^{\mathrm{P}}-E_{\text{BAC,RDOP}}^{\mathrm{P}}\\ \text{(A140)}& {\displaystyle }& {\displaystyle }-R_{\text{BAC}}^{\mathrm{P}}-G{Z}_{\text{BAC}}^{\mathrm{P}}-M_{\text{BAC}}^{\mathrm{P}}.\end{array}$$

A5 Microzooplankton processes

• Cellular quota (ratio):

  $$\begin{array}{}\text{(A141)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,MZ}}^{\mathrm{C}}=N_{\text{MZ}}/C_{\text{MZ}}\text{(A142)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,MZ}}^{\mathrm{C}}=C_{\text{MZ}}/P_{\text{MZ}}.\end{array}$$

• Gross growth:

  $$\begin{array}{}\text{(A143)}& {\displaystyle }& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{C}}=G{Z}_{\text{CR}}^{\mathrm{C}}+G{Z}_{\text{BAC}}^{\mathrm{C}}\text{(A144)}& {\displaystyle }& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{N}}=G{Z}_{\text{CR}}^{\mathrm{N}}+G{Z}_{\text{BAC}}^{\mathrm{N}}\text{(A145)}& {\displaystyle }& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{P}}=G{Z}_{\text{CR}}^{\mathrm{P}}+G{Z}_{\text{BAC}}^{\mathrm{P}}.\end{array}$$

• LDOM excretion:

  $$\begin{array}{}\text{(A146)}& {\displaystyle }& {\displaystyle }{E}_{\text{MZ,LDOC}}^{\mathrm{C}}=f_{\text{ex,MZ}}\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{C}}\text{(A147)}& {\displaystyle }& {\displaystyle }{E}_{\text{MZ,LDON}}^{\mathrm{N}}=f_{\text{ex,MZ}}\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{N}}\text{(A148)}& {\displaystyle }& {\displaystyle }{E}_{\text{MZ,LDOP}}^{\mathrm{P}}=f_{\text{ex,MZ}}\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{P}}.\end{array}$$

• SDOM excretion:

  $$\begin{array}{ll}\text{(A149)}& {\displaystyle }{E}_{\text{MZ,SDOC,1}}^{\mathrm{C}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,MZ}})\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{C}}{\displaystyle }{E}_{\text{MZ,SDON,1}}^{\mathrm{N}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,MZ}})\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{N}}\\ \text{(A150)}& {\displaystyle }& {\displaystyle }\times Q_{\text{N,MZ}}^{\mathrm{C}}/q_{\text{N,MZ}}^{\mathrm{C}}{\displaystyle }{E}_{\text{MZ,SDOP,1}}^{\mathrm{P}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,MZ}})\times {\mathit{\text{ex}}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{P}}\\ \text{(A151)}& {\displaystyle }& {\displaystyle }\times Q_{\text{P,MZ}}^{\mathrm{C}}/q_{\text{P,MZ}}^{\mathrm{C}}.\end{array}$$

• SDOM excretion to adjust stoichiometry:

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{MZ,SDOC,2}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,MZ}}\times C_{\text{MZ}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{0},\mathrm{1}-Q_{\text{N,MZ}}^{\mathrm{C}}/q_{\text{N,MZ}}^{\mathrm{C}},\\ \text{(A152)}& {\displaystyle }& {\displaystyle }\mathrm{1}-Q_{\text{P,MZ}}^{\mathrm{C}}/q_{\text{P,MZ}}^{\mathrm{C}})\text{(A153)}& {\displaystyle }{E}_{\text{MZ,SDON,2}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0.5}\times E_{\text{MZ,SDOC,2}}^{\mathrm{C}}\times Q_{\text{N,MZ}}^{\mathrm{C}}\text{(A154)}& {\displaystyle }{E}_{\text{MZ,SDOP,2}}^{\mathrm{P}}=& {\displaystyle }\mathrm{0.5}\times E_{\text{MZ,SDOC,2}}^{\mathrm{C}}\times Q_{\text{P,MZ}}^{\mathrm{C}}.\end{array}$$

• Remineralization of inorganic nutrients:

  $$\begin{array}{ll}{\displaystyle }{\text{REMI}}_{\text{MZ}}^{\mathrm{N}}=& {\displaystyle }{\text{remi}}_{\text{MZ}}\times max(\mathrm{0},N_{\text{MZ}}-C_{\text{MZ}}\times q_{\text{N,MZ}}^{\mathrm{C}},\\ \text{(A155)}& {\displaystyle }& {\displaystyle }{N}_{\text{MZ}}-q_{\text{N,MZ}}^{\mathrm{C}}/P_{\text{MZ}}\times q_{\text{P,MZ}}^{\mathrm{C}}){\displaystyle }{\text{REMI}}_{\text{MZ}}^{\mathrm{P}}=& {\displaystyle }{\text{remi}}_{\text{MZ}}\times max(\mathrm{0},P_{\text{MZ}}-C_{\text{MZ}}\times q_{\text{P,MZ}}^{\mathrm{C}},\\ \text{(A156)}& {\displaystyle }& {\displaystyle }{P}_{\text{MZ}}-q_{\text{P,MZ}}^{\mathrm{C}}/N_{\text{MZ}}\times q_{\text{N,MZ}}^{\mathrm{C}}).\end{array}$$

• Respiration:

  $$\begin{array}{}\text{(A157)}& R_{\text{MZ}}^{\mathrm{C}}=r_{\text{MZ}}^{\mathrm{B}}\times T_{\text{f}}\times C_{\text{MZ}}+r_{\text{MZ}}^{\mathrm{A}}\times G_{\text{MZ}}^{\mathrm{C}}.\end{array}$$

• POM production:

  $$\begin{array}{}\text{(A158)}& {\displaystyle }& {\displaystyle }{D}_{\text{MZ}}^{\mathrm{C}}={\text{pom}}_{\text{MZ}}\times G_{\text{MZ}}^{\mathrm{C}}\text{(A159)}& {\displaystyle }& {\displaystyle }{D}_{\text{MZ}}^{\mathrm{N}}=q_{\text{N,POM}}^{\mathrm{C}}\times G_{\text{MZ}}^{\mathrm{C}}\text{(A160)}& {\displaystyle }& {\displaystyle }{D}_{\text{MZ}}^{\mathrm{P}}=q_{\text{P,POM}}^{\mathrm{C}}\times G_{\text{MZ}}^{\mathrm{C}}.\end{array}$$

• Grazing by krill:

  $$\begin{array}{ll}{\displaystyle }G{Z}_{\text{MZ}}^{\mathrm{C}}=& {\displaystyle }{T}_{\text{f}}\times {\mathit{\mu }}_{\text{MZ}}\times C_{\text{KR}}\times [C_{\text{MZ}}^{\mathrm{2}}/\mathit{\{}{C}_{\text{MZ}}^{\mathrm{2}}+g_{\text{MZ}}^{\mathrm{2}}\\ \text{(A161)}& {\displaystyle }& {\displaystyle }+(C_{\text{DA}}\times g_{\text{MZ}}/g_{\text{DA}}{)}^{\mathrm{2}}\mathit{\}}]\text{(A162)}& {\displaystyle }G{Z}_{\text{MZ}}^{\mathrm{N}}=& {\displaystyle }{Q}_{\text{N,MZ}}^{\mathrm{C}}\times G{Z}_{\text{MZ}}^{\mathrm{C}}\text{(A163)}& {\displaystyle }G{Z}_{\text{MZ}}^{\mathrm{P}}=& {\displaystyle }{Q}_{\text{P,MZ}}^{\mathrm{C}}\times G{Z}_{\text{MZ}}^{\mathrm{C}}.\end{array}$$

• The net growth rate equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{MZ}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{C}}-E_{\text{MZ,LDOC}}^{\mathrm{C}}-E_{\text{MZ,SDOC,1}}^{\mathrm{C}}\\ \text{(A164)}& {\displaystyle }& {\displaystyle }-E_{\text{MZ,SDOC,2}}^{\mathrm{C}}-P_{\text{MZ}}^{\mathrm{C}}-R_{\text{MZ}}^{\mathrm{C}}-G{Z}_{\text{MZ}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{MZ}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{N}}-E_{\text{MZ,LDON}}^{\mathrm{N}}-E_{\text{MZ,SDON,1}}^{\mathrm{N}}\\ \text{(A165)}& {\displaystyle }& {\displaystyle }-E_{\text{MZ,SDON,2}}^{\mathrm{N}}-P_{\text{MZ}}^{\mathrm{N}}-R_{\text{MZ}}^{\mathrm{N}}-G{Z}_{\text{MZ}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{MZ}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{MZ}}^{\mathrm{P}}-E_{\text{MZ,LDOP}}^{\mathrm{P}}-E_{\text{MZ,SDOP,1}}^{\mathrm{P}}\\ \text{(A166)}& {\displaystyle }& {\displaystyle }-E_{\text{MZ,SDOP,2}}^{\mathrm{P}}-P_{\text{MZ}}^{\mathrm{P}}-R_{\text{MZ}}^{\mathrm{P}}-G{Z}_{\text{MZ}}^{\mathrm{P}}.\end{array}$$

A6 Krill processes

• Cellular quota (ratio):

  $$\begin{array}{}\text{(A167)}& {\displaystyle }& {\displaystyle }{Q}_{\text{N,KR}}^{\mathrm{C}}=N_{\text{KR}}/C_{\text{KR}}\text{(A168)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,KR}}^{\mathrm{C}}=C_{\text{KR}}/P_{\text{KR}}.\end{array}$$

• Gross growth:

  $$\begin{array}{}\text{(A169)}& {\displaystyle }& {\displaystyle }{G}_{\text{KR}}^{\mathrm{C}}=G{Z}_{\text{DA,KR}}^{\mathrm{C}}+G{Z}_{\text{MZ}}^{\mathrm{C}}\text{(A170)}& {\displaystyle }& {\displaystyle }{G}_{\text{KR}}^{\mathrm{N}}=G{Z}_{\text{DA,KR}}^{\mathrm{N}}+G{Z}_{\text{MZ}}^{\mathrm{N}}\text{(A171)}& {\displaystyle }& {\displaystyle }{G}_{\text{KR}}^{\mathrm{P}}=G{Z}_{\text{DA,KR}}^{\mathrm{P}}+G{Z}_{\text{MZ}}^{\mathrm{P}}.\end{array}$$

• LDOM excretion:

  $$\begin{array}{}\text{(A172)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,LDOC}}^{\mathrm{C}}=f_{\text{ex,KR}}\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{C}}\text{(A173)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,LDON}}^{\mathrm{N}}=f_{\text{ex,KR}}\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{N}}\text{(A174)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,LDOP}}^{\mathrm{P}}=f_{\text{ex,KR}}\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{P}}.\end{array}$$

• SDOM excretion:

  $$\begin{array}{ll}\text{(A175)}& {\displaystyle }{E}_{\text{KR,SDOC,1}}^{\mathrm{C}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,KR}})\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{C}}{\displaystyle }{E}_{\text{KR,SDON,1}}^{\mathrm{N}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,KR}})\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{N}}\\ \text{(A176)}& {\displaystyle }& {\displaystyle }\times Q_{\text{N,KR}}^{\mathrm{C}}/q_{\text{N,KR}}^{\mathrm{C}}{\displaystyle }{E}_{\text{KR,SDOP,1}}^{\mathrm{P}}=& {\displaystyle }(\mathrm{1}-f_{\text{ex,KR}})\times {\mathit{\text{ex}}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{P}}\\ \text{(A177)}& {\displaystyle }& {\displaystyle }\times Q_{\text{P,KR}}^{\mathrm{C}}/q_{\text{P,KR}}^{\mathrm{C}}.\end{array}$$

• SDOM excretion to adjust stoichiometry:

  $$\begin{array}{ll}{\displaystyle }{E}_{\text{KR,SDOC,2}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{ADJ,KR}}\times C_{\text{KR}}\\ {\displaystyle }& {\displaystyle }\times max(\mathrm{0},\mathrm{1}-Q_{\text{N,KR}}^{\mathrm{C}}/q_{\text{N,KR}}^{\mathrm{C}},\\ \text{(A178)}& {\displaystyle }& {\displaystyle }\mathrm{1}-Q_{\text{P,KR}}^{\mathrm{C}}/q_{\text{P,KR}}^{\mathrm{C}})\text{(A179)}& {\displaystyle }{E}_{\text{KR,SDON,2}}^{\mathrm{N}}=& {\displaystyle }\mathrm{0.5}\times E_{\text{KR,SDOC,2}}^{\mathrm{C}}\times Q_{\text{N,KR}}^{\mathrm{C}}\text{(A180)}& {\displaystyle }{E}_{\text{KR,SDOP,2}}^{\mathrm{P}}=& {\displaystyle }\mathrm{0.5}\times E_{\text{KR,SDOC,2}}^{\mathrm{C}}\times Q_{\text{P,KR}}^{\mathrm{C}}.\end{array}$$

• Remineralization of inorganic nutrients:

  $$\begin{array}{ll}{\displaystyle }{\text{REMI}}_{\text{KR}}^{\mathrm{N}}=& {\displaystyle }{\text{remi}}_{\text{KR}}\times max(\mathrm{0},N_{\text{KR}}-C_{\text{KR}}\times q_{\text{N,KR}}^{\mathrm{C}},\\ \text{(A181)}& {\displaystyle }& {\displaystyle }{N}_{\text{KR}}-q_{\text{N,KR}}^{\mathrm{C}}/P_{\text{KR}}\times q_{\text{P,KR}}^{\mathrm{C}}){\displaystyle }{\text{REMI}}_{\text{KR}}^{\mathrm{P}}=& {\displaystyle }{\text{remi}}_{\text{KR}}\times max(\mathrm{0},P_{\text{KR}}-C_{\text{KR}}\times q_{\text{P,KR}}^{\mathrm{C}},\\ \text{(A182)}& {\displaystyle }& {\displaystyle }{P}_{\text{KR}}-q_{\text{P,KR}}^{\mathrm{C}}/N_{\text{KR}}\times q_{\text{N,KR}}^{\mathrm{C}}).\end{array}$$

• Respiration:

  $$\begin{array}{}\text{(A183)}& R_{\text{KR}}^{\mathrm{C}}=r_{\text{KR}}^{\mathrm{B}}\times T_{\text{f}}\times C_{\text{KR}}+r_{\text{KR}}^{\mathrm{A}}\times G_{\text{KR}}^{\mathrm{C}}.\end{array}$$

• POM production:

  $$\begin{array}{}\text{(A184)}& {\displaystyle }& {\displaystyle }{D}_{\text{KR}}^{\mathrm{C}}={\text{pom}}_{\text{KR}}\times G_{\text{KR}}^{\mathrm{C}}\text{(A185)}& {\displaystyle }& {\displaystyle }{D}_{\text{KR}}^{\mathrm{N}}=q_{\text{N,POM}}^{\mathrm{C}}\times G_{\text{KR}}^{\mathrm{N}}\text{(A186)}& {\displaystyle }& {\displaystyle }{D}_{\text{KR}}^{\mathrm{P}}=q_{\text{P,POM}}^{\mathrm{C}}\times G_{\text{KR}}^{\mathrm{P}}.\end{array}$$

• RDOC release:

  $$\begin{array}{}\text{(A187)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,RDOC}}^{\mathrm{C}}={\text{refr}}_{\text{KR}}\times C_{\text{KR}}\text{(A188)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,RDON}}^{\mathrm{N}}=E_{\text{KR,RDOC}}^{\mathrm{C}}\times q_{\text{N,RDOM}}^{\mathrm{C}}\text{(A189)}& {\displaystyle }& {\displaystyle }{E}_{\text{KR,RDOP}}^{\mathrm{P}}=E_{\text{KR,RDOC}}^{\mathrm{C}}\times q_{\text{P,RDOM}}^{\mathrm{C}}.\end{array}$$

• Removal by higher trophic levels

  $$\begin{array}{}\text{(A190)}& {\displaystyle }& {\displaystyle }{M}_{\text{KR}}^{\mathrm{C}}={\text{mort}}_{\text{KR}}\times C_{\text{KR}}\times C_{\text{KR}}\text{(A191)}& {\displaystyle }& {\displaystyle }{M}_{\text{KR}}^{\mathrm{N}}=M_{\text{KR,RDOC}}^{\mathrm{C}}\times Q_{\text{N,KR}}^{\mathrm{C}}\text{(A192)}& {\displaystyle }& {\displaystyle }{M}_{\text{KR}}^{\mathrm{P}}=M_{\text{KR,RDOC}}^{\mathrm{C}}\times Q_{\text{P,KR}}^{\mathrm{C}}.\end{array}$$

• The net growth rate equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{KR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{KR}}^{\mathrm{C}}-E_{\text{KR,LDOC}}^{\mathrm{C}}-E_{\text{KR,SDOC,1}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }-E_{\text{KR,SDOC,2}}^{\mathrm{C}}-E_{\text{KR,RDOC}}^{\mathrm{C}}-P_{\text{KR}}^{\mathrm{C}}\\ \text{(A193)}& {\displaystyle }& {\displaystyle }-R_{\text{KR}}^{\mathrm{C}}-M_{\text{KR}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{KR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{KR}}^{\mathrm{N}}-E_{\text{KR,LDON}}^{\mathrm{N}}-E_{\text{KR,SDON,1}}^{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }-E_{\text{KR,SDON,2}}^{\mathrm{N}}-E_{\text{KR,RDON}}^{\mathrm{N}}-P_{\text{KR}}^{\mathrm{N}}\\ \text{(A194)}& {\displaystyle }& {\displaystyle }-R_{\text{KR}}^{\mathrm{N}}-M_{\text{KR}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{KR}}}{\mathrm{d}t}}=& {\displaystyle }{G}_{\text{KR}}^{\mathrm{P}}-E_{\text{KR,LDOC}}^{\mathrm{P}}-E_{\text{KR,SDOC,1}}^{\mathrm{P}}\\ {\displaystyle }& {\displaystyle }-E_{\text{KR,SDOC,2}}^{\mathrm{P}}-E_{\text{KR,RDOC}}^{\mathrm{P}}-P_{\text{KR}}^{\mathrm{P}}\\ \text{(A195)}& {\displaystyle }& {\displaystyle }-R_{\text{KR}}^{\mathrm{P}}-M_{\text{KR}}^{\mathrm{P}}.\end{array}$$

A7 Detrital processes

• Dissolution:

  $$\begin{array}{}\text{(A196)}& {\displaystyle }& {\displaystyle }{\text{DISS}}_{\text{DET}}^{\mathrm{C}}=\text{diss}\times C_{\text{DET}}\text{(A197)}& {\displaystyle }& {\displaystyle }{\text{DISS}}_{\text{DET}}^{\mathrm{N}}=\text{diss}\times {\text{prf}}_{\text{N}}\times N_{\text{DET}}\text{(A198)}& {\displaystyle }& {\displaystyle }{\text{DISS}}_{\text{DET}}^{\mathrm{P}}=\text{diss}\times {\text{prf}}_{\text{P}}\times P_{\text{DET}}.\end{array}$$

• The net change equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{DET}}}{\mathrm{d}t}}=& {\displaystyle }{D}_{\text{DA}}^{\mathrm{C}}+D_{\text{CR}}^{\mathrm{C}}+D_{\text{MZ}}^{\mathrm{C}}+D_{\text{KR}}^{\mathrm{C}}+D_{\text{HZ}}^{\mathrm{C}}\\ \text{(A199)}& {\displaystyle }& {\displaystyle }-D_{\text{DET}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{DET}}}{\mathrm{d}t}}=& {\displaystyle }{D}_{\text{DA}}^{\mathrm{N}}+D_{\text{CR}}^{\mathrm{N}}+D_{\text{MZ}}^{\mathrm{N}}+D_{\text{KR}}^{\mathrm{N}}+D_{\text{HZ}}^{\mathrm{N}}\\ \text{(A200)}& {\displaystyle }& {\displaystyle }-D_{\text{DET}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{DET}}}{\mathrm{d}t}}=& {\displaystyle }{D}_{\text{DA}}^{\mathrm{P}}+D_{\text{CR}}^{\mathrm{P}}+D_{\text{MZ}}^{\mathrm{P}}+D_{\text{KR}}^{\mathrm{P}}+D_{\text{HZ}}^{\mathrm{P}}\\ \text{(A201)}& {\displaystyle }& {\displaystyle }-D_{\text{DET}}^{\mathrm{P}},\end{array}$$

  where

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{D}_{\text{HZ}}^{\mathrm{C}}=f_{\text{POM,HZ}}\times M_{\text{KR}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }{D}_{\text{HZ}}^{\mathrm{N}}=f_{\text{POM,HZ}}\times M_{\text{KR}}^{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }{D}_{\text{HZ}}^{\mathrm{P}}=f_{\text{POM,HZ}}\times M_{\text{KR}}^{\mathrm{P}}.\end{array}$$

A8 DOM processes

• Conversion of SDOM to RDOM:

  $$\begin{array}{ll}{\displaystyle }{\text{REFR}}_{\text{SDOM}}^{\mathrm{C}}=& {\displaystyle }{\mathit{\text{ex}}}_{\text{REFR,SDOM}}\times C_{\text{SDOM}}\\ {\displaystyle }& {\displaystyle }\times \exp \mathit{\{}\mathrm{1}-min(Q_{\text{N,SDOM}}^{\mathrm{C}}/q_{\text{N,RDOM}}^{\mathrm{C}},\\ \text{(A202)}& {\displaystyle }& {\displaystyle }{Q}_{\text{P,SDOM}}^{\mathrm{C}}/q_{\text{P,RDOM}}^{\mathrm{C}}\left)\mathit{\right\}}\text{(A203)}& {\displaystyle }{\text{REFR}}_{\text{SDOM}}^{\mathrm{N}}=& {\displaystyle }{\text{REFR}}_{\text{SDOM}}^{\mathrm{C}}\times q_{\text{N,RDOM}}^{\mathrm{C}}\text{(A204)}& {\displaystyle }{\text{REFR}}_{\text{SDOM}}^{\mathrm{P}}=& {\displaystyle }{\text{REFR}}_{\text{SDOM}}^{\mathrm{C}}\times q_{\text{P,RDOM}}^{\mathrm{C}}.\end{array}$$

• The net change equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{C}_{\text{LDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,LDOC}}^{\mathrm{C}}+E_{\text{CR,LDOC}}^{\mathrm{C}}+E_{\text{MZ,LDOC}}^{\mathrm{C}}\\ \text{(A205)}& {\displaystyle }& {\displaystyle }+E_{\text{KR,LDOC}}^{\mathrm{C}}+M_{\text{BAC}}^{\mathrm{C}}-G_{\text{BAC,LDOC}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{LDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,LDON}}^{\mathrm{N}}+E_{\text{CR,LDON}}^{\mathrm{N}}+E_{\text{MZ,LDON}}^{\mathrm{N}}\\ \text{(A206)}& {\displaystyle }& {\displaystyle }+E_{\text{KR,LDON}}^{\mathrm{N}}+M_{\text{BAC}}^{\mathrm{N}}-G_{\text{BAC,LDON}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{LDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,LDOP}}^{\mathrm{P}}+E_{\text{CR,LDOP}}^{\mathrm{P}}+E_{\text{MZ,LDOP}}^{\mathrm{P}}\\ \text{(A207)}& {\displaystyle }& {\displaystyle }+E_{\text{KR,LDOP}}^{\mathrm{P}}+M_{\text{BAC}}^{\mathrm{P}}-G_{\text{BAC,LDOP}}^{\mathrm{P}}{\displaystyle \frac{\mathrm{d}{C}_{\text{SDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,SDOC}}^{\mathrm{C}}+E_{\text{CR,SDOC}}^{\mathrm{C}}+E_{\text{BAC,SDOC}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }+E_{\text{MZ,SDOC,1}}^{\mathrm{C}}+E_{\text{MZ,SDOC,2}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }+E_{\text{KR,SDOC,1}}^{\mathrm{C}}+E_{\text{KR,SDOC,2}}^{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }+\text{SDOC excretion by higher level}+{\text{DISS}}_{\text{DET}}^{\mathrm{C}}\\ \text{(A208)}& {\displaystyle }& {\displaystyle }-{\text{REFR}}_{\text{SDOM}}^{\mathrm{C}}-G_{\text{BAC,SDOC}}^{\mathrm{C}}{\displaystyle \frac{\mathrm{d}{N}_{\text{SDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,SDON}}^{\mathrm{N}}+E_{\text{CR,SDON}}^{\mathrm{N}}+E_{\text{BAC,SDON}}^{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }+E_{\text{MZ,SDON,1}}^{\mathrm{N}}+E_{\text{MZ,SDON,2}}^{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }+E_{\text{KR,SDON,1}}^{\mathrm{N}}+E_{\text{KR,SDON,2}}^{\mathrm{N}}\\ {\displaystyle }& {\displaystyle }+\text{SDON excretion by higher level}+{\text{DISS}}_{\text{DET}}^{\mathrm{N}}\\ \text{(A209)}& {\displaystyle }& {\displaystyle }-{\text{REFR}}_{\text{SDOM}}^{\mathrm{N}}-G_{\text{BAC,SDON}}^{\mathrm{N}}{\displaystyle \frac{\mathrm{d}{P}_{\text{SDOM}}}{\mathrm{d}t}}=& {\displaystyle }{E}_{\text{DA,SDOP}}^{\mathrm{P}}+E_{\text{CR,SDOP}}^{\mathrm{P}}+E_{\text{BAC,SDOP}}^{\mathrm{P}}\\ {\displaystyle }& {\displaystyle }+E_{\text{MZ,SDOP,1}}^{\mathrm{P}}+E_{\text{MZ,SDOP,2}}^{\mathrm{P}}\\ {\displaystyle }& {\displaystyle }+E_{\text{KR,SDOP,1}}^{\mathrm{P}}+E_{\text{KR,SDOP,2}}^{\mathrm{P}}\\ {\displaystyle }& {\displaystyle }+\text{SDOP excretion by higher level}+{\text{DISS}}_{\text{DET}}^{\mathrm{P}}\\ \text{(A210)}& {\displaystyle }& {\displaystyle }-{\text{REFR}}_{\text{SDOM}}^{\mathrm{P}}-G_{\text{BAC,SDOP}}^{\mathrm{P}}.\end{array}$$

A9 Dissolved inorganic nutrient processes

• Nitrification:

  $$\begin{array}{}\text{(A211)}& \text{NTRF}=r_{\text{ntrf}}\times {\mathrm{NH}}_{\mathrm{4}}.\end{array}$$

• The net change equations:

  $$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{d}t}}=& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{NH}}_{\mathrm{4}}}+{\text{REMI}}_{\text{MZ}}^{\mathrm{N}}+{\text{REMI}}_{\text{KR}}^{\mathrm{N}}\\ \text{(A212)}& {\displaystyle }& {\displaystyle }+{\text{REMI}}_{\text{HZ}}^{\mathrm{N}}-G_{\text{DA}}^{{\mathrm{NH}}_{\mathrm{4}}}-G_{\text{CR}}^{{\mathrm{NH}}_{\mathrm{4}}}-\text{NTRF}\text{(A213)}& {\displaystyle \frac{\mathrm{d}{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{d}t}}=& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{NO}}_{\mathrm{3}}}-G_{\text{DA}}^{{\mathrm{NO}}_{\mathrm{3}}}-G_{\text{CR}}^{{\mathrm{NO}}_{\mathrm{3}}}+\text{NTRF}{\displaystyle \frac{\mathrm{d}{\mathrm{PO}}_{\mathrm{4}}}{\mathrm{d}t}}=& {\displaystyle }{\text{FLUX}}_{\text{BAC}}^{{\mathrm{PO}}_{\mathrm{4}}}+{\text{REMI}}_{\text{MZ}}^{\mathrm{P}}+{\text{REMI}}_{\text{KR}}^{\mathrm{P}}\\ \text{(A214)}& {\displaystyle }& {\displaystyle }+{\text{REMI}}_{\text{HZ}}^{\mathrm{P}}-G_{\text{DA}}^{{\mathrm{PO}}_{\mathrm{4}}}-G_{\text{CR}}^{{\mathrm{PO}}_{\mathrm{4}}},\end{array}$$

  where

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{HZ}}^{\mathrm{N}}=M_{\text{KR}}^{\mathrm{N}}-D_{\text{HZ}}^{\mathrm{N}}-\text{SDON excretion by higher level}\\ {\displaystyle }& {\displaystyle }{\text{REMI}}_{\text{HZ}}^{\mathrm{P}}=M_{\text{KR}}^{\mathrm{N}}-D_{\text{HZ}}^{\mathrm{P}}-\text{SDOP excretion by higher level}.\end{array}$$

Appendix B

Figure B1 Comparison of the observations to the initial (unoptimized) and optimized model results. The dot points in the second panels represent how much larger model output value is compared to the corresponding observational data (i.e., the model value minus the observational value). Normalized observation: observations normalized by the mean of each model state variable.

Figure B2 Observations assimilated for each data type. BAC: bacterial biomass, CR: cryptophyte Chl, DA: diatom Chl.

Figure B3 The uncertainties (standard deviation) of the model state variables for the modeled growth season of 2002–2003 (x axis; month/day) from the Monte Carlo experiments (n=1000). Note the different contour scales among panels.

Figure B4 The uncertainties (standard deviation) of the ecosystem indices for the modeled growth season of 2002–2003 (x axis; month/day) from the Monte Carlo experiments (n=1000). Note the different contour scales among panels.

Table B1 Sensitivity tests to varying initial parameter values. The parameters with “n/a” in the parentheses are optimized parameters with large relative uncertainties (i.e., “optimized parameters”). The 50 % perturbations to A_E, wnsv, and remin (Table 1) are considered ecologically unrealistic and therefore excluded from these sensitivity experiments.

Table B2 Summary of model parameters from initial parameter perturbation experiments. Summary of the model parameter symbol and definition, initial parameter values (p_0,ref) and optimized values (p_f,ref) for the original (reference) experiments (Table 1), and initial parameter values (p₀) and optimized values (p_f) averaged from the sensitivity experiments (n=15). p_0,ref and the mean p₀ are the same for most parameters because of the perturbations by ±50 % their original initial parameter values (with standard deviation in parentheses), while p_0,ref and the mean p₀ are different for A_E and wnsv because of the perturbations only by −50 % their original initial parameter values (with standard deviation in parentheses). Numbers in parentheses for p_f are the uncertainty ranges (lower and upper bounds) averaged across the sensitivity experiments as follows. First, for each sensitivity experiment, lower and upper bounds of the constrained parameter are calculated as $p_{\mathrm{f}}\times e^{+{\mathit{\sigma }}_{\mathrm{f}}}$ and $p_{\mathrm{f}}\times e^{-{\mathit{\sigma }}_{\mathrm{f}}}$ (where p_f is the value of the constrained parameter and σ_f is the square roots of diagonal elements of the inverse of the Hessian matrix), respectively. Then we form the “lower (upper) bound parameter set” that only consists of the lower (higher) bounds of the constrained parameters from each experiment, and average those across the sensitivity experiments (n=15) to calculate the lower (upper) bound listed in parentheses.