A coupled pelagic–benthic–sympagic biogeochemical model for the Bering Sea: documentation and validation of the BESTNPZ model (v2019.08.23) within a high-resolution regional ocean model

The Bering Sea is a highly productive ecosystem, supporting a variety of fish, seabird, and marine mammal populations, as well as large commercial fisheries. Due to its unique shelf geometry and the presence of seasonal sea ice, the processes controlling productivity in the Bering Sea ecosystem span the pelagic water column, the benthic sea floor, and the sympagic sea ice environments. The Bering Ecosystem Study Nutrient-Phytoplankton-Zooplankton (BESTNPZ) model has been developed to simulate the lower-trophic-level processes throughout this region. Here, we present a version of this lower-trophic-level model coupled to a three-dimensional regional ocean model for the Bering Sea. We quantify the model's ability to reproduce key physical features of biological importance as well as its skill in capturing the seasonal and interannual variations in primary and secondary productivity over the past several decades. We find that the ocean model demonstrates considerable skill in replicating observed horizontal and vertical patterns of water movement, mixing, and stratification, as well as the temperature and salinity signatures of various water masses throughout the Bering Sea. Along the data-rich central portions of the southeastern Bering Sea shelf, it is also able to capture the mean seasonal cycle of primary production. However, its ability to replicate domain-wide patterns in nutrient cycling, primary production, and zooplankton community composition, particularly with respect to the interannual variations that are important when linking variation in productivity to changes in longer-lived upper-trophic-level species, remains limited. We therefore suggest that near-term application of this model should focus on the physical model outputs, while model development continues to elucidate potential mechanisms controlling nutrient cycling, bloom processes, and trophic dynamics.

1 Introduction

The Bering Sea is a highly productive ecosystem. Its broad, shallow eastern shelf reaches widths of over 500 km, with an average depth of only 70 m, leading to a long growing season and high annual primary production (Rho and Whitledge, 2007). Tidal mixing along the continental shelf break also leads to a highly productive off-shelf region, often referred to as the “Green Belt”, where the confluence of nitrate from the deep basin and iron from the shelf are mixed into the euphotic zone (Springer et al., 1996). This high primary productivity across the shelf and slope in turn supports a wide variety of pelagic and benthic predators, which support fisheries that land nearly half of the US annual catch (National Marine Fisheries Service, 2017; Fissel et al., 2017).

Because of the Bering Sea's economic and cultural importance, changes in its ecosystem have prompted a series of contemporary research efforts, including the Bering Ecosystem Study (BEST) and Bering Sea Integrated Ecosystem Research Program (BSIERP), that aimed to advance understanding of ecosystem processes and their relationship to the physical environment in the Bering Sea (Sigler et al., 2010). As part of these efforts, the Bering Ecosystem Study Nutrient-Phytoplankton-Zooplankton (BESTNPZ) biogeochemical model (Gibson and Spitz, 2011) was developed to simulate the key processes and features of the Bering Sea lower-trophic-level ecosystem, including primary and secondary production in the pelagic environment as well as benthic–pelagic and ice–pelagic interactions.

A regional ocean model that includes the BESTNPZ biological model has been used to investigate a variety of topics, including historical and future biophysical variability (Hermann et al., 2013, 2016, 2019), ecosystem status and variability (Ortiz et al., 2016), fish advection and recruitment (Wilderbuer et al., 2016), and community connectivity within crab populations (Richar et al., 2015); at least a dozen ongoing projects continue to rely on this model for retrospective and future analyses. Previous studies have examined the model's skill and sensitivity broadly. Hermann et al. (2013) demonstrated that the physical model shows skill in replicating observed patterns in temperature, salinity, and circulation, while Gibson and Spitz (2011) performed a thorough sensitivity analysis of the biogeochemical model in a one-dimensional environment representing the mid-shelf portion of the southeast Bering Sea shelf. However, a comprehensive evaluation of the biogeochemical skill of the BESTNPZ model in the three-dimensional ocean modeling context has been lacking.

Here, we present a thorough documentation of the BESTNPZ biogeochemical model in its current state. We also provide context and history for the various versions of the code that were used in the aforementioned publications and the changes that were made between publications and since. Finally, we evaluate several aspects of lower-trophic-level output of the BESTNPZ model within a high-resolution regional ocean model, focusing on whether the model properly captures key biophysical and biogeochemical processes necessary to realistically simulate primary production in the Bering Sea. This skill assessment reveals the model's strengths and weaknesses in reproducing historical patterns across the entire Bering Sea domain and also serves as a baseline to which further model improvements can be compared.

2 An overview of the BESTNPZ model

2.1 Model structure

The biogeochemical and ecosystem model underlying this study uses a traditional nutrient–phytoplankton–zooplankton structure. It tracks a total of 19 state variables: 14 pelagic variables, 2 benthic variables, and 3 sympagic (ice) variables.

The 14 pelagic state variables are resolved as tracer variables within the physical model and are therefore subject to advection and diffusion. The nutrients include nitrate, ammonium, and iron. Two size classes of phytoplankton (small and large) and five zooplankton groups (microzooplankton, small copepods, large copepods, euphausiids, and jellyfish) comprise the living state variables in the model. Both the large copepods and euphausiids groups are further subdivided into two state variables each, with parameterizations tailored to on-shelf and off-shelf populations; at present, the parameterizations for these two groups are very similar. Two detrital state variables, representing fast- and slow-sinking detritus, are also included.

The benthic submodel includes a living benthos group, encompassing all live infauna, and a single benthic detritus group. These state variables do not include any horizontal or vertical movement.

The three sympagic state variables are associated with seasonal sea ice and include nitrate, ammonium, and ice algae. These variables are assumed to occupy a thin skeletal layer on the underside of sea ice (when present); their horizontal movement is determined by the movement of the ice in which they reside, and they exchange material with the top layer of the ocean model. The exact thickness of the skeletal ice layer is specified via a user input parameter (see Sect. A3 for further details); for all simulations to date, a value of 0.02 m was used.

Nitrogen is used as the primary currency throughout the model, with all living and detrital state variables assumed to have a constant stoichiometry of 106 moles carbon per 16 moles nitrogen. While iron is also included as a state variable for primary production limitation purposes, its flux through the ecosystem is not tracked beyond its uptake during primary production, and water column iron is restored to a simple empirical distribution based on water depth on an annual timescale. Due to a quirk inherited from its predecessor model (Hinckley et al., 2009), many of the model's output variables, including phytoplankton, zooplankton, and detrital biomass variables as well as all flux rate diagnostic variables, are reported in carbon-based units; this conversion uses a constant N:C ratio across all state variables.

For a full description of all state variables, process equations, and input parameters in the BESTNPZ model, please see Appendix A.

In its three-dimensional setup, the BESTNPZ ecosystem model is coupled to an implementation of the Regional Ocean Modeling System (ROMS), a free-surface, primitive equation hydrographic model (Shchepetkin and McWilliams, 2005; Haidvogel et al., 2008). The Bering10K ROMS domain spans the Bering Sea and northern Gulf of Alaska, with 10 km horizontal resolution (Hermann et al., 2013). To date, it has been run with either 10 (previous studies) or 30 (this study) terrain-following depth levels (Fig. 1); we discuss this increase in vertical resolution in the next section. This horizontal domain is a subsection of a larger-domain ROMS model of the northeast Pacific (NEP5) (Danielson et al., 2011); both domains have been shown to capture the primary physical features important to biological processes, including circulation patterns, temperature, and salinity, as well as the seasonal advance and retreat of sea ice.

Figure 1 The Bering10K ROMS domain, including model bathymetry and a variety of place names within and near the domain. Dotted blue contour lines indicate the 50, 100, and 200 m isobaths; the solid blue lines indicate the approximate regions where these isobaths correspond to the fronts between the inner, middle, and outer domains of the southeastern shelf.

2.2 History and modifications

The BESTNPZ model has undergone several phases of development over the past several years. The code for the pelagic system originated from a Gulf of Alaska NPZ model known as GOANPZ (Hinckley et al., 2009). Model parameters and equations were tailored to the Bering Sea ecosystem during the Bering Ecosystem Study and Bering Sea Integrated Ecosystem Research Project, and the benthic and sympagic modules were added to the code during this phase. In the earliest publication of the model, Gibson and Spitz (2011) analyzed its sensitivity to input parameter uncertainty when coupled to a one-dimensional ocean model representing the Bering Sea M2 mooring location (56.87^∘ N, 164.06^∘ W). While this study confirmed that large-scale properties, such as modeled annual net primary production, were in line with observations from the eastern shelf region of the Bering Sea, it did not present any detailed skill analysis of the biological state variables against observations.

Following the Gibson and Spitz (2011) study, the BESTNPZ model was embedded within a realistic three-dimensional ocean model for the Bering Sea, referred to herein as the Bering10K ROMS domain. In a study focusing on physical and biological modes of variability with the Bering10K-BESTNPZ model output, Hermann et al. (2013) demonstrated that the model physics skillfully replicated observed patterns in temperature and salinity on the Eastern Bering Sea shelf, and that nutrients, phytoplankton, and zooplankton in this shelf region covaried with each other and with physical properties in a manner that supported the existing hypotheses of energy flow in the Eastern Bering Sea. However, Hermann et al.'s (2013) analysis focused on interannual variability and did not specifically assess the skill of the model to replicate seasonal patterns in primary and secondary productivity, nor did it address the behavior or skill of the model away from the eastern shelf. A few minor input parameter adjustments were made for the Hermann et al. (2013) study versus the earlier Gibson and Spitz (2011) study; apart from this and their differing one-dimensional versus three-dimensional physical environments, the versions of the BESTNPZ model used in these two publications are the same.

A third modification of the BESTNPZ model appears in Hermann et al. (2016), where an updated version of the Bering10K-BESTNPZ model complex seen in Hermann et al. (2013) was used to investigate long-term change in the Bering Sea under various climate change scenarios. In this version of the model, the light attenuation scheme was modified, and the light-related input parameters were adjusted. Several subsequent studies have used and continue to make use of the hindcast and/or climate forecast scenario output that was produced with this version of the Bering10K-BESTNPZ model (Ortiz et al., 2016; Hermann et al., 2019). Pilcher et al. (2019) also used this version of the model, with the addition of carbon variables, to support their analysis of ocean acidification in the Bering Sea.

Following the Hermann et al. (2016) publication, the BESTNPZ source code underwent a thorough revision. The revision process revealed a handful of small but consequential issues in the implementation of the BESTNPZ process equations that have led to a new version of the BESTNPZ model. These issues included the following:

1. Lack of conservation of biomass in large copepod biomass. An error in the code governing vertical migration of large copepods led to a slow but steady post-diapause accumulation of offshore large copepod (NCaO) biomass below their diapausing depth of 400 m. Over the course of a multi-year simulation, this deep biomass manifested itself as a “fall bloom” at depth and reached levels surpassing that of surface large copepod biomass; this artifact is visible in the depth-integrated biomass results presented in Ortiz et al. (2016). This error has now been corrected.

2. High light limitation. Preliminary validation of simulated phytoplankton production revealed that light was strongly limiting in certain regions of the model domain even under noon summer conditions. In deep water, the 10-layer physical model setup was too coarse to resolve variations in light within the mixed layer, leading to light-limited conditions year round. In the shallower regions of the inner shelf, an overly aggressive sediment attenuation term – introduced between the Hermann et al. (2013, 2016) versions of the model – also led to year-round light limitation. This high light limitation masked problems with micronutrient limitation in the deep basin and macronutrient limitation in the inner coastal domain. To remove this excessive light limitation, new runs of the Bering10K domain now use 30 vertical layers rather than 10, and a new equation for light attenuation was implemented. See Appendix B for further discussion of these changes to light attenuation.

3. Non-conservative behavior of macronutrients related to nudging. Gibson et al.'s (2011) version of the BESTNPZ model implemented empirical relaxation of iron, ammonium, and nitrate throughout its one-dimensional domain; both nitrate and iron were relaxed towards seasonal climatological profiles, while ammonium was relaxed toward zero, all on annual timescales. The relaxation was appropriate in the one-dimensional context, given the use of periodic lateral boundary conditions that did not allow for any advective transport of nutrients into or out of the domain. However, when moved to the three-dimensional Bering10K domain (as was done in all previous publications, including Hermann et al., 2013, 2016), this nudging becomes inappropriate, as all processes controlling nutrient distribution should conserve biomass (total nitrogen in the system is not constant, due to the open lateral boundary conditions and out-of-system fluxes from burial and benthic denitrification, but these changes to the nitrogen budget can be attributed to known processes and quantified). Particularly in the case of ammonium relaxation, the nudging applied to the three-dimensional domain introduced a phantom process that “scavenged” ammonium from the water column across most of the shelf regions. In the most recent versions of the model, nudging has been removed from the NO₃ and NH₄ state variables; it remains in place for Fe, since the simplified iron dynamics in the current model do not supply any explicit sources of bioavailable iron.

4. Preferential uptake of nitrate in high-nitrate, high-ammonium conditions. Under the previous governing equations for macronutrient limitation during the gross primary production calculations, the total nitrogen limitation factor (i.e., the factor applied to the maximum photosynthetic uptake rate to account for nutrient limitation, after Lomas and Glibert, 1999) could exceed 1.0 under high-nitrate, high-ammonium conditions; the code implemented a cap to counter this by reducing the ammonium limitation factor so the sum of the nitrate and ammonium limitation factors was less than 1. This approach led to reduced uptake of ammonium in favor of nitrate when concentrations of both macronutrients were high, and is counter to the assumption that ammonium uptake is usually energetically favored over nitrate uptake, and that high levels of ammonium inhibit nitrate uptake (Glibert et al., 2016, and references therein); it is also counter to the ammonium inhibition encoded in the nitrate limitation factor itself in this model. While this quirk in nutrient uptake in unlikely to have made a large difference in uptake dynamics given that it was only relevant during nutrient-replete conditions, it may have exacerbated the accumulation of ammonium on the shelf seen when nudging was removed. In the most recent version of the code, nitrate and ammonium limitation factors were modified to use an equation (after Frost and Franzen, 1992) that constrains the total nitrogen limitation factor to a range of 0–1 without the need for additional caps.

5. Euphausiid prey preferences. In previous versions of the BESTNPZ model, the euphausiid groups preyed upon large phytoplankton, ice algae, microzooplankton, and small copepods. However, the modeled euphausiid populations tended to drop precipitously in winter months when these prey were scarce, in contrast to fish diet data that indicate the continued presence of euphausiids in fish diets on the southeastern shelf during the winter (Livingston et al., 2017). This observation is important to replicate when using the BESTNPZ model to look at the broader flow of energy to higher-trophic-level species (e.g., Ortiz et al., 2016). In an attempt to increase overwintering success of the on-shelf euphausiid group, feeding links to the two detrital groups and to small phytoplankton were added. The addition of these feeding links for the onshore euphausiid group distinguishes the tendency of on-shelf Bering Sea euphausiids, dominated by Thysanoessa raschii, to rely on detrital feeding for overwinter survival; in contrast, the shelf-edge population, dominated by Thysanoessa inermis, accumulates higher lipid stores to support winter survival (Hunt et al., 2016).

In addition the changes to the biogeochemical model listed above, significant biases in the ice fields produced by earlier versions of Bering10K and related models have been identified and addressed by ROMS colleagues (Durski and Kurapov, 2019, Kate Hedstrom and Alexander Kurapov, personal communications, May 2016). Late melting of ice was noted by Danielson et al. (2011) and Cheng et al. (2014) for a northeast Pacific model that utilizes nearly the same ice code as Bering10K, and by Ortiz et al. (2016) and Hermann et al. (2013, 2016) for the Bering10K model itself. Corrections to the longwave radiation terms of ice thermodynamics have been implemented in the latest version of the Bering10K code to address the late ice melting bias. Additional improvements to the mechanical ice dynamics by Hedstrom (Kate Hedstrom,personal communication, 2016) corrected for a previous bias towards excessive ice thickness in some areas. Specifically, this was corrected through the addition of a quadratic ice strength term that resists ice ridging, based on the work of Overland and Pease (1988).

3 Model validation: methods

3.1 Model configuration and forcing

We ran two simulations of the Bering10K-BESTNPZ model for this study. For both simulations, the model is driven by surface atmospheric forcing from either the Common Ocean Reference Experiment (CORE) (Large and Yeager, 2009) (1970–1994), the Climate Forecast System Reanalysis (CFSR) (Saha et al., 2010) (1995–March 2011), or the Climate Forecast System Operational Analysis (CFSv2-OA) (April 2011–September 2018); bulk formulae were used to relate winds, air temperature, specific humidity, surface pressure, and shortwave and longwave radiation from these datasets to surface stress, freshwater, and heat fluxes (Fairall et al., 1996). Comparison between overlapping years from the CORE and CFSR datasets revealed small differences in values in the radiation variables; the CORE shortwave and longwave radiation values were therefore divided by factors of 0.9 and 0.97, respectively, to align with the CFSR data (note that this is the opposite of the adjustment performed in previous studies, e.g., Hermann et al. (2013, 2016), where the CFSR values were adjusted downward by 10 % and 3 %, respectively). Lateral boundary conditions for the open southern and eastern boundaries of the model domain use a hybrid nudging–radiation scheme (Marchesiello et al., 2001). During the CORE period, these boundary conditions were derived from a simulation of the larger northeast Pacific grid (Danielson et al., 2011), which relied on the Simple Ocean Data Assimilation (SODA) dataset (Carton and Giese, 2008) for its own lateral boundary conditions; the CFS periods use CFS ocean variable values. The northern boundary transport through the Bering Strait is relaxed to the observed value of 0.8 Sv (Woodgate and Aagaard, 2005); earlier sensitivity studies tested whether a seasonally varying open boundary condition could better replicate the flow patterns in the northern portion of the domain, but the simple relaxation condition was found to perform equally well. Freshwater runoff due to river input was reconstructed from observed river discharge from Alaskan and Russian rivers (Kearney, 2019); river freshwater input is distributed across model grid points near the coast as a surface freshwater flux based on river mouth location, with an e-folding scale of 20 km.

The first model simulation, referred to hereafter as the spinup simulation, looped interannually invariant surface and lateral boundary conditions over a 30-year period. We chose to use boundary condition values from 2001, a year with close to average sea ice cover. Physical variables for this simulation were initialized from the 1 January 2001 values of a previous hindcast of the Bering10K domain. Nitrate was initialized to a constant value of 40 mmol N m⁻³ below 300 m, transitioning to 0 mmol N m⁻³ at 100 m. Iron was initialized to the same empirical profile used for annual nudging within the model, which sets surface and deep iron values based on bottom depth (see Appendix A). All living biological state variables (i.e., phytoplankton, zooplankton, and benthic infauna) were initialized using a tiny seed value to allow future growth, while all other state variables were initialized at zero. The purpose of this simulation was to allow the model to reach an internally regulated state, and to verify that any accumulation of nitrogen outside the deep basin was a result of internal dynamics rather than overestimated initial conditions.

The second simulation, referred to as the hindcast simulation, was initialized using values from the final time step of the spinup simulation. The model was then run from 1970 to 2018, forced with the full time series of surface and lateral boundary conditions from the combined CORE/CFS-derived dataset described above.

3.2 Key features of biological importance

To systematically validate the Bering10K-BESTNPZ model complex, we focus on a few key features of the Bering Sea. We begin by looking at physical processes that are known to influence the primary production and then compare our modeled patterns of primary production, phytoplankton biomass, and phytoplankton and zooplankton community composition to a variety of measurements. We primarily focus this evaluation on the Eastern Bering Sea shelf due to the availability of data in this region but also qualitatively evaluate the patterns seen in the central basin and the northern shelf regions. While our physical model domain extends into the northern Gulf of Alaska, the biological model was never intended to simulate this region, and for this validation we assume that all regions south of the Aleutian Islands or east of the Alaska Peninsula are outside the biological domain of the BESTNPZ model.

3.2.1 Sea ice

Sea ice plays a key role in shaping the ecosystems of the Bering Sea. Ice advances southwestward through the Bering Strait into the Bering Sea, driven both by winds from the northeast and local ice formation, with much of the eastern shelf at least partially covered by ice beginning in early winter (October to November) through early spring (March to April). Variability in the timing of ice onset and retreat and extent of sea ice can be significant year to year, influenced by winds and air temperature related to the position and strength of the Aleutian Low and Siberian High pressure systems, as well as ocean conditions (Stabeno et al., 2001).

To analyze the extent and timing of the updated seasonal sea ice, we collected estimates of sea ice concentration derived from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data (Cavalieri et al., 1996) for the period of 1980–2018. For comparison with model output, satellite-derived fraction ice cover was interpolated from its native 6.25–20 km resolution polar stereographic grid to the Bering10K model grid via a nearest-neighbor method. For both model and observations, we calculate the location of the ice edge along 170^∘ W (the approximate longitudinal center of the seasonal sea ice in the Bering Sea), defining the edge as the southern extent of a continuous block of grid cells where all cells have at least 15 % ice cover.

3.2.2 The cold pool

As ice advances southward along the Bering Sea shelf, the freezing process and resulting brine rejection leads to the formation of cold, salty, dense bottom water underneath the ice (Stabeno et al., 2001); continuous freezing in the vicinity of the St. Lawrence polynya further intensifies the formation of this bottom water mass (Danielson et al., 2006). The resulting cold bottom water is referred to as the “cold pool” and is typically defined as waters colder than either 0 or 2 ^∘C. In the spring, warming of the surface waters coupled with melting of sea ice sharply stratifies the water column over the middle shelf region, isolating the cold pool waters from surface heating and mixing. As a result, this signature water mass can persist well into the summer months (Stabeno et al., 2001). Due to the isolation of the cold pool, this water can serve as a nutrient reservoir to the mid-shelf region when mixed with nutrient-depleted surface waters during storm events following the initial spring bloom (Sambrotto et al., 1986; Aguilar-Islas et al., 2007). In addition, the location of the cold pool influences the spatial distribution (Mueter and Litzow, 2008; Stabeno et al., 2012) and recruitment of higher-trophic-level predators (Mueter et al., 2011; Duffy-Anderson et al., 2016) through various mechanisms.

Measurements of bottom temperature are collected each summer as part of the Bering Sea Groundfish Bottom Trawl Survey. Net trawls are conducted at fixed survey stations located across the Eastern Bering Sea shelf at 20 nautical mile resolution. From 1982 to 1989, temperature data were collected via expendable bathythermographs (XBTs). More recent surveys use digital bathythermograph recorders attached to the headrope of the bottom trawl net (BRANCKER RBR XL-200 Micro BTs recorded at 6 s intervals for the 1993–2001 surveys, and a Sea-Bird SBE-39 bathythermograph continuous data recorder at 3 s intervals for 2002–present). Bottom temperature is then averaged over the on-bottom portion of the trawl to produce a single value per station per year (see Buckley et al., 2009; Lauth et al., 2011, for details of data collection and postprocessing). Here, we use bottom temperature measurements to verify that our model properly captures the characteristics of the Bering Sea cold pool. The cold pool is quantified by indices that represent the fraction of the survey area with bottom water less than 0, 1, or 2 ^∘C water. In the model, we define the Eastern Bering Sea shelf as all grid cells falling within the Eastern Bering Sea stratified sampling regions 10–62 (Fig. 2). We calculate the model cold pool index using two methods. First, we calculate the index value on 1 July of each year; choosing a fixed date allows us to compare a consistent summer snapshot of the cold pool from year to year. The second index replicates the sampling scheme used in the groundfish survey, choosing bottom temperatures from the closest grid cell and time slice to each observation point; this index allows a better comparison to the observations, given the temporal spread in the observations between the first (typically southeast) and last (northwest) sampled station. For comparison, we also analyze bottom temperature extracted from the Climate Forecast System ocean model; this lower-resolution climate model is coupled to the same atmospheric forcings we use in our hindcast simulation for this time period and therefore allows comparison between the original climate model and our dynamically downscaled representation.

Figure 2 Bering Sea groundfish survey stratified sampling polygons. The area considered for calculation of the cold pool index in the model is shaded. The grid cells for each trawl, including bottom temperature measurement samples, are indicated by the light gray lines (circles on the corner of grid cells indicate crab resampling locations; data points from these locations were used in the spatial interpolations shown in Fig. 4 but were removed from cold pool index calculations to avoid unequal weighting).

3.2.3 Inner, middle, and outer shelf domains

The wide, shallow Eastern Bering Sea shelf is divided into three domains (Fig. 1), each characterized by distinctive patterns of stratification and mixing (Coachman, 1986; Kachel et al., 2002; Stabeno et al., 2012). The inner domain stretches from the coast to approximately the 50 m isobath and is well mixed year round by both tidal and wind energy. The middle domain reaches from the 50–100 m isobath; this domain is well mixed during winter months but thermally stratified during the spring and summer, with a tidally mixed bottom layer isolated from the surface waters. The outer shelf domain, reaching from the 100 m isobath to the shelf break (approximately the 200 m isobath), more closely reflects the seasonal stratification of an oceanic system, with a tidally mixed bottom layer that is less sharply separated from the surface layer than in the middle domain (Coachman and Charnell, 1978). Tidal mixing dominates the energy across the entire shelf, with a very small net transport northward (Coachman, 1986).

We estimate stratification in the model by calculating the potential energy required to mix the water column (SI, in units of J m⁻²), after Simpson et al. (1977):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{SI}={\displaystyle \frac{\mathrm{1}}{\mathit{\zeta }+h}}\underset{-h}{\overset{\mathit{\zeta }}{\int }}g(\mathit{\zeta }-z)(\mathit{\rho }-\stackrel{\mathrm{‾}}{\mathit{\rho }})\mathrm{d}z,\text{(2)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\rho }}={\displaystyle \frac{\mathrm{1}}{\mathit{\zeta }+h}}\underset{-h}{\overset{\mathit{\zeta }}{\int }}\mathit{\rho }\mathrm{d}z,\end{array}$$

where ρ is density, h is the depth of the water column, ζ is the height of the free surface, g is gravitational acceleration, and z is depth relative to mean surface height (i.e., ζ=0).

We also calculated the location of the structural front separating the well-mixed inner domain from the thermally stratified middle domain, defining the front as the 0.5 ^∘C m⁻¹ contour in maximum vertical temperature gradient (after Schumacher et al., 1979; Kachel et al., 2002). We apply this calculation to the years 2000–2010 in the hindcast simulation; this period spans multi-year cold and warm periods, and therefore encompasses a good deal of the variability one might expect to see in this property.

3.2.4 Spatial and temporal patterns in primary production

The physical geometry of the Bering Sea, along with the seasonal presence of sea ice, leads to a diverse set of controls on primary production, varying in both space and time. For validation purposes, we focused on a few features that best reflect this complex balance of macronutrient, micronutrient, light, and temperature control of bottom-up processes in this ecosystem.

The highest sustained productivity in the Bering Sea is seen near the edge of the shelf break. This region, referred to as the “Green Belt”, coincides spatially with both the shelf break and the Bering Slope Current that carries water northward along the shelf slope (Springer et al., 1996). The shelf-break front and eddies drive high primary productivity both by supplying macronutrients (i.e., nitrate) from the deep basin and micronutrients (i.e., iron) from the shelf and by physically entraining phytoplankton (Okkonen et al., 2004). The Pribilof Islands provide an additional source of dissolved iron to the Green Belt region (Aguilar-Islas et al., 2007). Variability in mesoscale eddies in the Bering Slope Current is a primary driver of productivity variability in the Green Belt (Okkonen et al., 2004; Mizobata and Saitoh, 2004). The initial spring bloom here is dominated by diatoms but transitions to smaller species in the late summer (Springer et al., 1996).

Another hotspot of production in this region is the Pribilof Islands. Their location disrupts flow along the 100 m isobath, which leads to enhanced tidal mixing and introduces nutrients from the deep basin to this area, leading to high summer productivity. Productivity can be lower during years when mixing is decreased due to salinity-related frontal structures propagating from the inner shelf (Stabeno et al., 2008).

While the Bering Sea is generally characterized as being very productive, this production is almost entirely driven by the on-shelf regions and the Green Belt. The deep basin, in contrast, is a high-nutrient, low-chlorophyll (HNLC) region, with low iron levels limiting primary productivity year round (Aguilar-Islas et al., 2007; Suzuki et al., 2002; Leblanc et al., 2005).

On the wide eastern shelf, primary productivity is mainly controlled by the balance of stratification-induced changes in light and nitrogen limitation. An initial spring bloom, dominated by diatoms, rapidly depletes the surface macronutrients along much of the shelf (Sambrotto et al., 1986). In the marginal ice zone, ice-edge blooms can occur, accounting for a large fraction of the total spring bloom (Niebauer et al., 1995). This ice-edge bloom occurs during years where ice lingers later over the shelf, protecting the underlying water from wind mixing and setting up stronger stratification; earlier-melting ice leads to more wind mixing and a later spring bloom. Variations in spring bloom timing, its correlation or lack thereof with ice melt date, and the impact of this timing on community composition and energy transfer to higher trophic levels form the backbone of most prevailing theories of ecosystem function in the southeastern Bering Sea (Hunt et al., 2010; Sigler et al., 2016). Later summer and fall productivity can be driven by wind mixing events that mix nutrient-rich bottom water into the surface layer. Measurements along the middle and outer shelf regions indicate that nitrate drawdown accounts for 30 %–50 % of observed carbon uptake, with the remaining 50 %–70 % driven by ammonium (Sambrotto et al., 1986; Whitledge et al., 1986).

To check for these patterns in primary production, we performed a visual comparison of modeled phytoplankton biomass patterns with satellite chlorophyll estimates. The satellite chlorophyll values used were climatological monthly averages from MODIS Aqua Ocean Color Index (OCI) algorithm 4 km resolution product, spanning the period of July 2002–July 2015 (NASA Ocean Biology Processing Group, 2017). We compared these chlorophyll patterns to optically weighted, depth-integrated phytoplankton chlorophyll from the model over the same time period, assuming an attenuation length scale of 45 m. While chlorophyll is not a direct measurement of biomass, particularly when derived from satellite color in high-latitude locations, and our model-derived estimate of satellite-visible chlorophyll is a rough one, it is in this case sufficient to allow large-scale comparison of general spatial patterns in biomass between the various biophysical regions of the Bering Sea (e.g., basin versus shelf, north versus south).

We also looked at chlorophyll data measured via fluorometer at the long-term biophysical mooring at station M2 (56.87^∘ N, 164.05^∘ W) (Stabeno et al., 2012). This mooring provides a more detailed look at both surface and subsurface chlorophyll over several years, including during times of ice cover.

3.3 Plankton community composition

To evaluate plankton community composition, we focused on a few patterns of relative biomass seen in the Bering Sea. The spring bloom is typically dominated by diatoms, with small phytoplankton numbers increasing in summer and fall (Springer et al., 1996). Microzooplankton, consisting primarily of protists such as ciliates and dinoflagellates, are the primary grazers on both large and small phytoplankton size classes; measurements of the biomass of microzooplankton vary on order of 10–100 mg C m⁻³ (Olson and Strom, 2002; Howell-Kübler et al., 1996). The mesozooplankton community is dominated by larger zooplankton. Though numerically abundant, the small copepod species typically compose less than 10 % of the zooplankton biomass (Vidal and Smith, 1986). Amongst the larger zooplankton groups, the dominant species vary on and off the shelf. The offshore community is dominated by oceanic copepod and euphausiid species, such as Neocalanus spp. and Thysanoessa inermis (NCaO and EupO in the model, respectively), while on the shelf region Calanus marshallae and Thysanoessa raschii (NCaS and EupS, respectively) compose the majority of the mesozooplankton population. Biomass measurements for these larger mesozooplankton groups vary on the order of 1–10 mg C m⁻³ (Campbell et al., 2016; Hunt et al., 2016). The offshore and onshore euphausiid groups are distinguished by differing survival strategies. The shelf-edge T. inermis build up high lipid stores used for overwinter survival and spawn early in the spring, without the need to feed on the spring bloom for spawning, while the shelf-dwelling T. raschii rely more heavily on detrital feeding during the winter and spawn later, after feeding on the spring bloom (Hunt et al., 2016, and references therein).

4 Model validation: results

4.1 Sea ice

The improved sea ice model in this version of the Bering10K model demonstrates high skill in reproducing the advance of sea ice across the domain and in capturing the interannual variability of the location of maximum ice extent (Fig. 3). Over the entire 1980–2018 time series covered by satellite observations, the ice-edge location along 170^∘ W shows a small southerly bias of 0.3^∘ latitude (33.4 km) compared to the location measured via satellite; this is improved from the previous 0.56^∘ (62.3 km) southerly bias seen in the Hermann et al. (2016) version of the ice model. However, despite improvements compared to previous versions of the model, ice retreat still lags observations by approximately 2 weeks in the early spring.

Figure 3 Skill assessment of the location of the sea ice edge along 170^∘ W. Panel (a) depicts the entire time series in satellite observations, the Hermann et al. (2016) version of the Bering10K model, and the current version of the Bering10K model. Panels (b) and (c) indicate the seasonal and interannual mismatch between the location of the ice edge in the model and the observations. Panel (d) shows the seasonal median (line) and 25th–75th percentile ranges (shaded region) of ice-edge extent across the three datasets.

4.2 The cold pool

The Bering10K bottom temperature values clearly reproduce both spatial and temporal variability in the location and extent of the cold pool (Fig. 4). During low-ice years, the cold pool is primarily located in the northwest portion of the eastern shelf, while in colder, high-ice years it extends throughout much of the middle domain and into Bristol Bay.

Figure 4 Comparison of observed, Bering10K-simulated, and CFS-simulated bottom temperature on the eastern shelf under a variety of conditions, including 2004 (a warm year), 2008 (a cold year), 2016 (an atypical warm year, where anomalously warm water from the Gulf of Alaska pushed the cold pool further northwest than usual), and 2017 (one of the few survey years when samples were collected from the northern Bering Sea). Discrete samples from the groundfish survey and the survey-replicated model datasets (indicated by black dots) were interpolated to the model grid using natural neighbor interpolation; the black contour line indicates the 2 ^∘C edge of the cold pool.

For a more quantitative assessment of skill, we looked at the correlation between annual time series of mean bottom temperature and cold pool index values in the groundfish survey dataset versus the models (Table 1, Fig. 5). The Bering10K model values correlate very strongly with the observed values; correlation is highest when calculating the cold pool extent as defined by the 0 ^∘C threshold (R²=0.940), and R² values remain above 0.87 for all other metrics. To summarize model skill, we use a model efficiency metric after Stow et al. (2009), where a model efficiency value greater than 0 indicates more skill than a simple average of the observed time series, and a value of one indicates a perfect match to the observed time series. The model efficiency for the Bering10K model is high across all metrics, ranging from 0.714 to 0.8. This is in sharp contrast to the mean bottom temperature and cold pool index metrics estimated by the coarser-resolution Climate Forecast System model. The cold pool produced by the CFS model lacks the detailed structure of colder waters seen in the observations and fails to simulate bottom waters below the 1 ^∘C threshold. The model efficiency metric (MEF) suggests that the CFS model has much less skill in predicting mean interannual bottom temperature (MEF of 0.406), with only a marginal ability to capture the 2 ^∘C cold pool and no skill at all with respect to waters less than 1 ^∘C. This indicates that the dynamic downscaling offered by the higher-resolution Bering10K model is a necessary component in reproducing this feature.

Figure 5 Annual indices of observed versus modeled bottom temperature, including average bottom temperature in the eastern shelf survey area, fraction of the survey area less than 0 ^∘C, and fraction of the survey area less than 2 ^∘C.

Table 1 Correlation, bias, root mean squared difference (RMSD), and model efficiency for each model estimate of mean bottom temperature or cold pool index compared to the observed temperature/index. Statistics are applied to annual time series.

4.3 Inner, middle, and outer shelf domains

Before analyzing whether the model properly reproduces the variations in vertical structure and mixing in the three shelf domains, it is important to note that the location of isobaths in our model are offset slightly compared to the real shelf bathymetry. Sigma-coordinate models such as the ROMS model used in this study are subject to computational errors in the horizontal pressure gradient along regions where topography is steep or the vertical gradient in a property is large (Shchepetkin, 2003); this issue often necessitates applying a smoothing filter to the bathymetry (Danielson et al., 2011). As a result, our modeled outer domain, as defined purely by location of isobaths, is generally narrower than that seen in any observations (Fig. 6), particularly in the vicinity of the Pribilof Islands.

Figure 6 Difference in location between the real (black) and modeled (green) locations of the 50, 100, and 200 m isobaths in the model domain. Colors indicate the adjustments to model bottom bathymetry (m) compared to the ETOPO5 dataset from which it was derived.

The simulated patterns in vertical stratification follow those expected across the three domains (Fig. 7). During the winter, the majority of the shelf is well mixed vertically. Stratification first appears in early spring along the outer domain and by summer is also seen throughout the middle domain. Given the relatively coarse vertical resolution of our model, the distinction in the vertical structure between the middle and outer domains is not well resolved. However, the model does reproduce the structural front, also known as a tidal front, expected between the unstratified inner domain and thermally stratified middle domain during the summer months. The exact location of this front varies both seasonally and from one year to the next but is generally located just inside the 50 m isobath (Fig. 8). The front location is relatively consistent across the southern shelf region, though possibly further inshore than would be predicted by the 50 m isobath near Cape Newenham at the north end of Bristol Bay; it moves further offshore and its location becomes more temporally variable further north. The clear structural front begins to disappear north of Nunivak Island. Stratification in the northern Bering Sea and Norton Sound area is much more strongly influenced by salinity, especially near the large outflows of the Yukon River, and in this region the clear demarcation between inner, middle, and outer domains disappears.

Figure 7 Monthly climatologically averaged simulated stratification across the Eastern Bering Sea shelf. Dark red lines indicate the primary 50, 100, and 200 m contour lines based on the model shelf bathymetry.

Figure 8 Location of structural fronts, defined as the 0.5 ^∘C m⁻¹ contour line of maximum vertical temperature gradient. Each blue line indicates a front location estimated from a single weekly averaged time point taken from the August and September values between 2000 and 2010. Contours indicate model domain bathymetry. Note that fronts over water deeper than 200 m may reflect artifacts of the coarsening vertical resolution rather than true changes in vertical gradients.

Horizontal movement in the model, as expected, is dominated by tidal frequencies across the shelf domain, with low annually averaged net velocities. There is a small net counterclockwise flow along the southern edge of the eastern shelf and then northward within the inner domain, with a small net transport from off-shelf to inner shelf waters (Fig. 9). Water entering the Bering Sea from the Gulf of Alaska through Unimak Pass moves alongside and onto the eastern shelf and travels northward; it takes approximately 7–8 months to reach the northern shelf region (i.e., 60^∘ N) along the 100 m isobath, in line with drifter-derived measurement of this flow (Stabeno et al., 2016). Further north, modeled velocities are slightly slower than those seen in the observations, with water taking approximately 13–15 months to reach the Bering Strait from the Unimak Pass in the model compared to 9–13 months in the observations; this may reflect a weak Anadyr Current in this region or alternatively be the result of missing flow from off-shelf through submarine canyons that are not well resolved by the modeled bathymetry. Overall, flow within the modeled Bering Sea reproduces the important circulation patterns within this region. However, cyclical circulation patterns seen near the southern and eastern boundaries of the model domain are likely an artifact of the boundary conditions.

Figure 9 Mean currents at 40 m depth (or in the bottom layer, for locations shallower than 40 m) across the model domain, averaged over 2000–2010. Blue lines indicate bottom contours of 50, 100, and 200 m. Shading indicates the magnitude of the flow. Light gray lines indicate flow streamlines.

4.4 Spatial and temporal patterns in primary production

Satellite ocean color measurements suggest that phytoplankton blooms in the Bering Sea first reach observable levels of chlorophyll in late February to early March, primarily on the eastern shelf in regions recently vacated of ice (Fig. 10). As light levels and temperature increase throughout the domain in summer, chlorophyll levels increase both on the shelf and along the shelf slope but remain low over the western side of the basin, where iron is limiting. The bloom peaks in late May to early June, then steadily decreases through September. A late fall bloom, smaller in magnitude than the earlier spring bloom, can be seen in October along the eastern shelf.

Figure 10 Monthly mean, optical-depth-averaged chlorophyll in the Bering10K-BESTNPZ model versus satellite-estimated values.

While the BESTNPZ model produces annual cycles of primary productivity of approximately the correct magnitude compared to these observations, it does not capture many of the nuances in spatiotemporal variability (Fig. 10). In early spring, the model does not appear to capture the early ice-associated growth along the eastern shelf. Within the pelagic ecosystem, low growth rates governed by low temperature-mediated maximum production rates coupled with strong light limitation prevent any significant accumulation of phytoplankton. Measurements at the M2 mooring location suggest that peak spring bloom date varies widely, from mid-April to early June (Sigler et al., 2014); in the model, peak bloom timing is constricted to a much narrower window from early to late May. While concentrations within the ice algae state variable can reach approximately 70 mg m⁻³ (monthly climatological average) within the thin skeletal ice layer, this biomass does not contribute significantly to the pelagic large phytoplankton concentration once ice melts due to dilution coupled with unfavorable growth conditions in the underlying water. We do not include these ice algal numbers in the optically weighted chlorophyll numbers used in our satellite comparison because these satellite measurements do not typically capture the spectral signals of ice algal pigments (Wang et al., 2018).

The spring bloom in the model begins once light and temperature levels increase in April. The first stages of the bloom resemble observations, with concentrations highest along the shelf slope and along the western shelf. However, rather than producing a short-lived bloom that drops once macronutrients are exhausted, the model allows for a sustained summer bloom. This bloom is driven by regenerated production; ammonium is produced from phytoplankton and zooplankton respiration, as well as quick remineralization of egested detritus, especially of the slow-sinking detritus group fed by small phytoplankton non-predatory mortality and microzooplankton egestion and non-predatory mortality. This pattern is seen both on the eastern shelf and throughout much of the deep basin, where iron does not appear to play its expected role in limiting primary production. In the basin, production levels fall off as macronutrients are exhausted in early July; on the eastern shelf, however, high fluxes of ammonium from the benthos drive sustained production throughout the summer and into early fall. In late fall, modeled chlorophyll levels appear more similar to the satellite patterns, with production primarily limited to the middle domain of the eastern shelf.

The late bias in the spring phytoplankton bloom is also evident when comparing model output to mooring measurements at the M2 mooring (Fig. 11). The model captures the predominant bloom characteristics: the bloom begins with a large, diatom-dominated bloom starting in the near-surface waters and then migrating deeper as surface nutrients are depleted, then decreases during the summer months, with some subsurface chlorophyll remaining at the bottom of the mixed layer. However, the modeled bloom begins in mid-April, on average, later than the mid-March to early April start seen in the mooring data. Fall blooms in September and October are spurred by increased mixing and are short and localized in the observations; the modeled fall bloom matches the timing in the observations well.

Figure 11 (a) Fluorescence-derived chlorophyll measurements at the M2 mooring location. Black dots indicate individual measurements, while the shaded values indicate estimates using a spring model interpolant (D'Errico, 2016). (b) Modeled chlorophyll extracted at the M2 mooring location over the same time period as the mooring data (2004–2017). (c) Surface-only climatological time series at the M2 location from mooring-based measurements, via satellite (within a 1^∘ box around the M2 location derived from MODIS Aqua OCI algorithm 4 km 8 d average images) and in the BESTNPZ model.

4.5 Plankton community composition

The phytoplankton community composition in the model reflects the expected balance between the small and large functional groups (Fig. 12). Throughout the majority of the domain, the spring bloom is heavily dominated by the large phytoplankton group. The one exception to this is along the shallow, well-mixed inner domain, where low macronutrient levels favor the small phytoplankton group for most of the year, with a small contribution of large phytoplankton in early summer when the bloom begins there. In the inner part of the middle shelf, the large phytoplankton levels decrease following the initial bloom but small phytoplankton biomass continues to rise through the fall. Moving further outward along the shelf, late summer and early fall biomass is low across both functional groups.

Figure 12 Depth-integrated biomass of phytoplankton and zooplankton groups in the southeastern Bering Sea shelf during the final year of the spinup simulation. Colors indicate the biomass value from each functional group, while the dashed black line indicates the total across all phytoplankton (a) or zooplankton (b) groups. Values were extracted from the model at the five locations indicated by gold dots in Fig. 1 (including the M2 mooring location). All axes use the same scale, labeled in nitrogen units on the left and carbon units on the right.

Within the zooplankton community, we see little variation between the relative dominance of the functional groups across the shelf transect (indicated by gold circles and labels in Fig. 1). In all locations, microzooplankton are the dominant group. However, their biomass is often only slightly higher than that of the summed mesozooplankton groups. Within the mesozooplankton groups, very little variation in their relative contribution to the biomass pool is seen either spatially or temporally. The only big change to zooplankton community coincides with the hardcoded diapause movement of the two large copepod groups; because these groups cease grazing when they enter diapause, their populations quickly drop during the diapause period. Offshore large copepods (NCaO) die if they encounter the ocean floor prior to reaching their prescribed overwintering depth of 400 m, and this effectively keeps this functional group constrained to deep-water locations. Lacking any similar depth-based restrictions on their process rates, the remaining large zooplankton groups can be found throughout the domain, regardless of whether they are designated as onshore or offshore in name.

The model does capture a gradient in timing of the zooplankton population, with offshore populations being established early in the spring, while on-shelf populations do not appear until early summer. However, limited observations suggest that early spring offshore zooplankton increases precede the spring phytoplankton bloom offshore (Hunt et al., 2016; Harvey et al., 2012). This timing difference is not captured by the model; instead, the overwintering population in the model is reduced to negligible amounts and only begins to increase again once the primary productivity in the region reaches sufficient levels.

5 Discussion

The Bering10K model correctly reproduces a variety of physical processes known to influence primary production in the Bering Sea. Overall patterns of sea ice cover, including interannual variations in the maximum extent of sea ice and the date of ice retreat, are well captured by the model. The exact date of sea ice retreat tends to lag observations by approximately 2 weeks. This lag could, theoretically, lead to subsequent problems in the timing of phytoplankton blooms; in the model's current state, ice algae and ice-associated blooms are poorly resolved regardless of ice melt accuracy or lack thereof, but improving the ice retreat timing should remain a focus of future improvements to the physical model.

The location and extent of the cold pool in the Bering Sea is often used as an index of biophysical variability across the Bering Sea shelf (e.g., Siddon and Zador, 2018), and therefore capturing both spatial and temporal variation in cold pool extent is key for a model to be useful in this region. The Bering10K model performs well on this point, with very high correlation to observations and very small biases. It is encouraging to note that the location of the simulated cold pool in 2016 is offset to the northwest compared to other warm years, replicating the position seen in the observations. During this year, anomalously warm water from the northeast Pacific (Bond et al., 2015) prevented the cold pool from extending as far southeast as it typically does during the summer months. The ability of the model to capture these anomalous conditions lends promise to its capability to simulate novel conditions that may arise when simulating future conditions.

The Bering10K model shows good replication of both cross-shelf and along-shelf differences in stratification. In the south, we see a distinct well-mixed inner domain, with a sharp transition to a stratified middle/outer domain occurring near the 50 m isobath during summer. The model's distinction between the middle and outer domains is less defined than in observations, likely due to the limited 30 vertical layers used in the model combined with bathymetric smoothing. In the northern portions of the eastern shelf, these thermal stratification domains disappear. Salinity is more variable in the north than in the south, driven by both sea ice formation and melt, as well as the large freshwater contribution of the Yukon River. As a result, stratification in the north is driven more strongly by salinity than in the south.

While the physical dynamics of the model perform well within the Bering Sea itself, the modeling domain south of the Aleutian Islands should be treated with caution, as it is close to the model's open boundary and artifacts associated with boundary conditions are expected.

Despite its strong skill in replicating the underlying physical features of the Bering Sea that are thought to influence primary production, the Bering10K-BESTNPZ model has limited skills in reproducing observed spatial and temporal patterns of primary production. The degraded performance in the biology realm is due to several interacting deficiencies in model process equations and parameterizations.

Throughout the deep basin, according to observations, iron levels should be low and limit primary production. The model includes only a simple representation of iron, using continuous relaxation to an empirical depth profile. While the low surface concentration prescribed in our model's basin is consistent with observations, this mechanism of replenishment through relaxation is not one that reflects the true complexity of iron cycling in the ocean. In their observations of phytoplankton growth rates in the Green Belt along the slope, Aguilar-Islas et al. (2007) noted that even in this highly productive location, the diatoms showed signs of iron stress, and dissolved iron levels remained low here compared to the shelf. They hypothesized that production in this region was maintained at its observed level due to a small but persistent source of iron being mixed from deep water along the shelf break, rather than a large iron source that fully alleviated iron limitation. The climatological nudging used in our model provides exactly this – a small but persistent source of continuous dissolved iron – throughout the domain, rather than only along the shelf slope. In order to properly capture the HNLC characteristics of the deep basin, a more mechanistic model for iron, with an explicit source near the sediments only rather than throughout the water column, is likely necessary.

Across the eastern shelf, in terms of mean seasonal cycles, modeled phytoplankton biomass and primary production levels are more in line with observations and reflect the dominant seasonal pattern observed on the southeast middle shelf of a strong spring bloom, followed by low summer biomass, then a smaller late fall bloom. However, many of the prevailing hypotheses of energy flow in the Bering Sea focus not on the mean state of the phytoplankton bloom but rather on interannual variability, particularly related to the interplay between temperature, stratification, and nutrient availability during the initial stages of the spring bloom (Stabeno et al., 2001; Coyle et al., 2008; Hunt et al., 2010). In particular, the timing of the spring bloom, and its correlation or lack thereof with ice retreat timing, forms the basis for many theories of energy transfer within the eastern Bering Sea shelf ecosystem. Given the key role that phenological variability plays in the predominant theories of energy transfer, shortcomings in the model's ability to capture the processes leading to such variability raise concerns about its potential ability to predict either current or future changes in primary and secondary production.

The first issue is that the spring bloom begins nearly a month later than it should. This is particularly apparent in the north (Fig. 10), where observations indicate that phytoplankton blooms should occur both on and under the sea ice. In our model, ice algae biomass is insignificant compared to pelagic phytoplankton biomass, and pelagic production is too strongly limited by both temperature and light for any significant growth to begin during March or April as it should. Currently, the lack of early spring growth also leads to a failure of the model to react to interannual differences in ice extent and retreat timing.

The limit on ice algae biomass in our model is primarily a limitation of the conceptual framework imposed on this state variable. By assuming that the ice algae are confined within the very thin skeletal ice layer, while at the same time allowing a continuous convective exchange (see Sect. A3.8) between this layer and the surface layer of the water column, modeled ice algae can never grow to much higher concentrations than would be found in an equivalent thin layer of water. This framework does not account for other aggregation types often seen in ice algal communities, such as nets or strands on the underside of the ice surface (Ambrose et al., 2005). Once the ice melts and the ice algae are released into the water, the model framework immediately transfers this pool of biomass to the large phytoplankton group, where it is then subject to the same controls on growth rate as the pelagic-originating phytoplankton.

For pelagic phytoplankton groups, the inability to capture early ice-associated blooms is primarily an inadequacy of the equations and parameters chosen to represent photosynthetic processes. The parameters that define each group's photosynthesis–irradiance curve, as well as those setting the maximum potential light- and nutrient-replete growth rates (a function of temperature), originated from a comparable model for the Gulf of Alaska (Hinckley et al., 2009). However, the phytoplankton community of the Bering Sea includes many Arctic species that are physiologically adapted to grow in both lower temperatures and under a wider range of light levels experienced near and under the sea ice. For example, ice-associated blooms occur in very thin upper layers of the water column left behind by ice melt; these layers are typically only 1–2 m thick and characterized by temperatures of around −1.7 ^∘C (Hunt et al., 2010). At this temperature, the model equations require approximately 2 W m⁻² surface irradiance to balance respiration and non-predatory mortality costs, even in the absence of any nutrient limitation or grazing losses. But modeled under-ice surface irradiance typically remains below this level when ice of more than approximately 0.5 m thickness is present. Therefore, the chosen set of equations do not appear appropriate to reproduce the dynamics of under-ice and ice-edge blooms. We also note that while the parameters and equations controlling the maximum potential growth rate of phytoplankton (Eq. A16) produce reasonable rates within the temperature ranges seen in the hindcast period in this geographical domain, they increase exponentially above this range, well outpacing respiration rate increases; a temperature increase of 5–10 ^∘C would push these rates well beyond the physiological limits of phytoplankton division rates.

Another key problem seen for phytoplankton across the shelf is the absence of any strong macronutrient limitation following the initial spring bloom. Observations indicate that ammonium can reach high quantities (up to 15 mg m⁻³) following the initial bloom due to both phytoplankton decomposition and benthic processes (Whitledge et al., 1986; Aguilar-Islas et al., 2007). However, this ammonium is typically concentrated in the deeper shelf water in the observations, while the model tends to accumulate material in a subsurface layer, with continuous high turnover of ammonium in surface waters of the shelf from spring through fall. This likely indicates that decomposition processes in the model are proceeding too quickly, particularly in the slow-sinking detrital group whose coupled remineralization rates and sinking rates result in this material only reaching about 50 m in depth before being completely converted to ammonium. Small phytoplankton mortality and excretion as well as excretion and egestion by microzooplankton are the primary sources for the slow-sinking detritus group. The model also produces a very robust zooplankton community that persists from the start of the spring bloom until well into the winter months (mid-December), and whose byproducts of respiration, excretion, and egestion feed into this rapid, continuous cycle of regenerated primary production (see the Supplement). In general, the process equations concerning the sinking and remineralization of organic material, both in the water column and on the seafloor, are very simplistic (see Sect. A3.7 for details). A single timescale for remineralization is used for all detrital groups, with no distinction between the lability of different pools of organic material. There is also no mechanism available to account for aggregation of material, which could lead to faster sinking speeds and lower remineralization rates of organic matter. Because of the broad, shallow nature of the Eastern Bering Sea shelf, remineralization rates play a strong role in determining the concentrations of macronutrients and the strength of benthic–pelagic coupling. These processes should be a focus of future development in this model.

In contrast to the under-resolved detrital pools, the mesozooplankton groups included in the BESTNPZ model appear to be over-resolved in terms of functional differences capable of being differentiated from each other with this type of biomass box model. The prey preferences encoded into the feeding behaviors of each group produce only very small variations in their relative contribution to the overall biomass pool. Instead, the five mesozooplankton groups (small copepods, on- and offshore large copepods, and on- and offshore euphausiids) effectively function as a single herbivorous–microzooplanktivorous zooplankton functional group. The only deviation from this synchrony is seen in the large copepod groups, whose populations fall when they enter diapause and cease grazing. The structure of the model also leads to plankton dynamics that are constrained primarily by the balance of instantaneous production and loss rates. While appropriate for simulating phytoplankton bloom dynamics, this style of model is not well suited for capturing the dynamics of zooplankton life stages (for example, the necessity of overwinter survival in order to spawn a new generation in the spring) or the nuanced differences between the survival strategies of various species. Mesozooplankton populations drop to a negligible level during the winter due to the absence of sufficient prey to balance continued respiratory and non-predatory losses, and overwintering success (or lack thereof) has almost no effect on the resulting zooplankton populations in summer. In order to truly capture the gradients in relative success of different copepod and euphausiid groups throughout the domain, a model that better captures winter survival strategies (e.g., a life-stage-resolving model or another means of introducing latency between feeding, respiration, and non-predatory mortality) is likely necessary.

Given the deficiencies identified in this evaluation, future work will comprehensively reevaluate each component of the existing model. More accurate simulation of under-ice and near-ice phytoplankton blooms may be addressed by allowing seasonal plasticity in the parameters defining the photosynthesis–irradiance curve for each phytoplankton group; when used in a simple nutrients–phytoplankton–zooplankton–detritus (NPZD)-style model, this type of equation has been shown to better capture the magnitude and timing of Bering Sea blooms than constant parameters (Sloughter et al., 2019). For sea ice algae, Tedesco and Vichi (2014) note that models using a fixed-thickness skeletal ice layer tend to underestimate production in first-year ice; they suggest that varying the width of the sea ice layer in which algae is found as a function of sea ice permeability can help overcome this issue with minimal additional model complexity required. Issues related to excessive regenerated production on the eastern shelf may be addressed by more closely examining the detrital functional groups within the model and the remineralization timescales associated with each; the use of a single remineralization timescale for all detrital groups is out of step with most modern biogeochemical models (e.g., Moore et al., 2002; Aumont and Bopp, 2006; Dunne et al., 2012), and allowing for parameters that reflect the varying lability of different detrital pools may better capture the nutrient dynamics both on and off the shelf. Improving the EBS nutrient budget may also require a more complex representation of the benthic component of the ecosystem; the benthic module from a mature shelf model such as ERSEM (Butenschön et al., 2016) may offer a blueprint for future development related to benthic functional groups. Finally, we intend to reconsider the number of functional groups used to represent the planktonic consumers within this ecosystem. Banas et al. (2016) demonstrated that a much simpler six-box model was capable of capturing spring bloom dynamics representative of the M2 mooring location. However, Friedrichs et al. (2007) cautioned that though simple models are typically able to be tuned to better simulate the ecosystem dynamics of a single location, their portability is limited compared to their more complex counterparts. Given the rapidly changing conditions in the Bering Sea, and the wide range of applications for which this model was designed (ranging from hindcast-based process studies to long-term climate change forecasts), we must carefully consider the tradeoffs of parsimony versus complexity.

6 Conclusions

Overall, the BESTNPZ model coupled to the Bering10K regional ocean model demonstrates considerable skill in replicating observed horizontal and vertical patterns of water movement, mixing, and stratification, as well as the temperature and salinity signatures of various water masses throughout the Bering Sea. However, its ability to replicate large-scale patterns in nutrient cycling, primary production, and zooplankton community composition, particularly with respect to the interannual variations that are important in a fisheries management context, is limited. In its current form, the Bering10K model can offer key insights into the physical processes that may affect higher-trophic-level species directly. In particular, it offers a useful supplement to examine physical features in areas and at times of the year that are difficult or impossible to survey due to sea ice cover or harsh weather. However, we caution that the use of the biological state variable output should be limited until the model is better able to capture observed characteristics of the Bering Sea phytoplankton and zooplankton communities.

Appendix A: Documentation for the Bering Sea Ecosystem Study Nutrient Phytoplankton Zooplankton Model (BESTNPZ)

A1 Summary and notation

This section provides a mathematical overview of processes through which biological state variables exchange material with each other in the BESTNPZ model.

The BESTNPZ model assumes a model geometry that includes N water column layers, a single benthic layer of unspecified depth, and a skeletal ice layer with a constant thickness $h_{s_{\mathrm{ice}}}$. The skeletal ice layer refers to the base of an ice sheet; this very thin layer is characterized by a looser crystal structure than the more solid ice overlying it and is the site of the most rapid algal growth in ice (Arrigo et al., 1993; Jin et al., 2006). Within this geometry, BESTNPZ tracks the concentration of 19 biological state variables (Table A1, Fig. A1).

Exchange of material between these state variables, and across vertical layers within a single state variable, results from a variety of processes. In the code, and in this description, these processes are divided into three types.

The first type of process, described in Sect. A2, includes redistribution of state variables due to movement of the water or ice in which they reside. The majority of these calculations (e.g., advection and diffusion of water and ice) take place outside the biological module and follow the default ROMS behavior for biological tracer variables. The one exception is the exchange of NO₃ and IceNO3, NH₄ and IceNH4, and PhL and IcePhL due to the formation or loss of ice in a grid cell.

The second process type we term source-minus-sink processes (Sect. A3); these processes take place within a single depth layer and involve transfer of biomass from one state variable to another. For notation, each source-minus-sink flux process is represented in this document as a function of the source and sink state variables, respectively. For example, Abc(X, Y) is the flux rate of material from group X to group Y via the Abc process.

The final category (Sect. A4) is vertical movement, where the concentration of state variables is redistributed within the water column due to the sinking or rising movement of the state variables within the water (note that this is separate from vertical advection of tracers due to movement of the water itself).

The three types of processes are calculated sequentially in the code, such that changes due to ice loss and formation are calculated first, followed by the total rate of change due to source-minus-sink processes, and finally redistribution due to vertical movement.

Several of the equations in this section rely on state or diagnostic variables that come from the physical model or from the ROMS grid geometry. See Table A2 for a description of these variables and their notation in this document. Additional parameters derived from biological input parameters are listed in Table A3.

To express the concentration of a biological state variable X, we use the [X] notation, with units corresponding to those in Table A1. All intermediate fluxes are expressed in terms of mg C for simplicity. Conversions between units assume constant stoichiometry for all living and detrital groups.

Because many of the variable names used in these equations involve multi-letter and mixed-case notations, we have chosen to use dot notation for all instances of multiplication in this document. Please note that this indicates simple element-by-element multiplication in the BESTNPZ code, not a dot product.

A2 Ice formation and loss

Although the primary ROMS sea ice model tracks ice presence in terms of fraction grid cell coverage, the BESTNPZ biological module uses a simpler scheme where ice presence is treated as a binary condition. When the ice thickness (h_ice) in a grid cell is greater than the prescribed thickness of the skeletal ice layer ($h_{s_{\mathrm{ice}}}$) and the grid cell has at least 50 % ice cover (as indicated by the a_ice variable), we assume the entirety of that grid cell now supports the ice-related biological processes in a thin layer of skeletal ice covering the entire grid cell and located just above the free surface of that grid cell.

If sea ice appears in a grid cell between the previous time step and the current time step, the large phytoplankton, nitrate, and ammonium in the top layer are redistributed such that the top water column layer and the skeletal ice layer receive equal concentrations of each tracer by volume.

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }[\mathrm{IcePhL}{]}_{\mathrm{ice}}={\displaystyle \frac{[\mathrm{PhL}{]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\text{(A2)}& {\displaystyle }& {\displaystyle }[\mathrm{PhL}{]}_N={\displaystyle \frac{[\mathrm{PhL}{]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\text{(A3)}& {\displaystyle }& {\displaystyle }[\mathrm{IceNO}\mathrm{3}{]}_{\mathrm{ice}}={\displaystyle \frac{{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\text{(A4)}& {\displaystyle }& {\displaystyle }{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_N={\displaystyle \frac{{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\text{(A5)}& {\displaystyle }& {\displaystyle }[\mathrm{IceNH}\mathrm{4}{]}_{\mathrm{ice}}={\displaystyle \frac{{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\text{(A6)}& {\displaystyle }& {\displaystyle }{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}_N={\displaystyle \frac{{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}_N\cdot h_N}{h_{s_{\mathrm{ice}}}+h_N}}\end{array}$$

Likewise, when ice disappears between the previous step and the current one, all material in the skeletal ice layer is moved to the top layer of the water.

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }[\mathrm{PhL}{]}_N={\displaystyle \frac{[\mathrm{PhL}{]}_N\cdot h_N+[{\mathrm{IcePhL}}_{\mathrm{ice}}]\cdot h_{s_{\mathrm{ice}}}}{h_N}}\text{(A8)}& {\displaystyle }& {\displaystyle }{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_N={\displaystyle \frac{{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}_N\cdot h_N+\left[\mathrm{IceNO}{\mathrm{3}}_{\mathrm{ice}}\right]\cdot h_{s_{\mathrm{ice}}}}{h_N}}\text{(A9)}& {\displaystyle }& {\displaystyle }{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}_N={\displaystyle \frac{{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}_N\cdot h_N+\left[\mathrm{IceNH}{\mathrm{4}}_{\mathrm{ice}}\right]\cdot h_{s_{\mathrm{ice}}}}{h_N}}\text{(A10)}& {\displaystyle }& {\displaystyle }[\mathrm{IcePhL}{]}_{\mathrm{ice}}=\mathrm{0}\text{(A11)}& {\displaystyle }& {\displaystyle }[\mathrm{IceNO}\mathrm{3}{]}_{\mathrm{ice}}=\mathrm{0}\text{(A12)}& {\displaystyle }& {\displaystyle }[\mathrm{IceNH}\mathrm{4}{]}_{\mathrm{ice}}=\mathrm{0}\end{array}$$

A3 Source-minus-sink processes

Unless otherwise specified, all processes detailed in this section are specific to a single layer. We have used the subscript k when defining each layer-specific flux rate; k can be either the index of a water column layer or “sice” to indicate the skeletal ice layer. This notation distinguishes between rates that are specific to a single layer and use volumetric units (mg C m⁻³ d⁻¹) and those that exchange material between layers and are expressed in total flux across a boundary (mg C m⁻² d⁻¹). To avoid clutter, we have chosen not to apply these subscripts to any remaining layer-dependent variables (such as state variable concentrations), but these are to be assumed for all pelagic and ice layer components.

A3.1 Light attenuation in water

The model assumes that radiation is attenuated with depth as follows:

$$\begin{array}{}\text{(A13)}& I_z=f_{\mathrm{PAR}}\cdot I_{\mathrm{0}}\cdot \exp \left(-K_{\mathrm{PAR}}\cdot (\mathit{\zeta }-z)\right),\end{array}$$

where f_PAR is the fraction of surface light that is photosynthetically available, I₀ is the surface irradiance, and K_PAR is the light attenuation coefficient for photosynthetically active radiation (i.e., 400–700 nm). Incoming radiation is supplied by the physical model and converted to photon flux.

The attenuation coefficient is itself the sum of attenuation from clear water, chlorophyll, and other sediment and organic material:

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{K}_{\mathrm{PAR}}& =K_{\mathrm{w}}+K_{\mathrm{A}}\cdot {\left({\displaystyle \frac{\left[\mathrm{PhL}\right]}{{\mathrm{ccr}}_{\mathrm{L}}}}+{\displaystyle \frac{\left[\mathrm{PhS}\right]}{{\mathrm{ccr}}_{\mathrm{S}}}}\right)}^{K_{\mathrm{B}}}+K_{\mathrm{C}}\\ & +K_{{\mathrm{D}}_{\mathrm{1}}}\cdot h^{K_{{\mathrm{D}}_{\mathrm{2}}}}.\end{array}\end{array}$$

The first two terms in Eq. (A14) are derived from Morel et al.'s (1988) analysis of light attenuation in Case I waters. The final terms (i.e., the K_D portion plus K_C constant) add additional attenuation based on the depth of the water column; this approximates the assumption that sediment and organic material is higher near the coastline than in open water. The power law formula was chosen based on a fit to satellite-derived inherent optical properties across the Bering Sea domain (see Appendix B for further details).

The K_PAR parameter is also used to calculate shortwave radiation decay in the physical model. By default, ROMS uses the equations of Paulson and Simpson (1977), which consider the differing attenuation length scales for blue-green light versus shorter- and longer-length wavelengths that are primarily attenuated in the upper 5–10 m of the water column. When coupling to the BESTNPZ model, we modify this equation to substitute our custom PAR attenuation length scale for the blue-green portion of the spectrum:

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}{f}_{\mathrm{swdk}}& =\left(\mathrm{1}-a_{\mathrm{frac}}\right)\cdot \exp \left(-K_{\mathrm{PAR}}\cdot (\mathit{\zeta }-z)\right)\\ & +a_{\mathrm{frac}}\cdot \exp \left({\displaystyle \frac{-(\mathit{\zeta }-z)}{a_{{\mathit{\mu }}_{\mathrm{1}}}}}\right).\end{array}\end{array}$$

The values of a_frac and $a_{{\mathit{\mu }}_{\mathrm{1}}}$ correspond to R and ζ₂ in Paulson and Simpson (1977), with our K_PAR replacing Paulson et al.'s (1977) ζ₁; we use the parameter values for Case I waters.

A3.2 Gross primary production

Primary production for both small and large phytoplankton is governed by the same set of equations. The maximum photosynthetic growth rate per unit chlorophyll (mg C (mg Chl a)⁻¹ d⁻¹) is a function of temperature and defined in terms of each group's doubling rate D_i and doubling rate exponent D_p (Frost, 1987). The maximum uptake rate is calculated in both carbon-specific and chlorophyll-specific units:

$$\begin{array}{}\text{(A16)}& {\displaystyle }& {\displaystyle }{P}_{\max }=\left({\mathrm{2}}^{\left(D_{\mathrm{i}}\times {\mathrm{10}}^{\left(D_{\mathrm{p}}\cdot T\right)}\right)}-\mathrm{1}\right)\text{(A17)}& {\displaystyle }& {\displaystyle }{P}_{\max }^{*}=P_{\max }\cdot \mathrm{ccr}.\end{array}$$

This rate is moderated by light and nutrient limitation. Light limitation uses a hyperbolic tangent function, after Jassby and Platt (1976):

$$\begin{array}{}\text{(A18)}& {\mathrm{Lim}}_I=\tanh \left({\displaystyle \frac{\mathit{\alpha }\cdot I_z}{P_{\max }^{*}}}\right).\end{array}$$

Nutrient limitation is based on the availability of nitrate, iron, and ammonium. Nitrate and ammonium limitation terms follow Frost and Franzen (1992), with nitrate uptake inhibited by ammonium when the latter is high relative to its half-saturation parameter:

$$\begin{array}{}\text{(A19)}& {\displaystyle }& {\displaystyle }{\mathrm{Lim}}_{{\mathrm{NO}}_{\mathrm{3}}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\left(k_{\mathrm{1}}+\left[{\mathrm{NO}}_{\mathrm{3}}\right]\right)\left(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{k_{\mathrm{2}}}\right)}}\text{(A20)}& {\displaystyle }& {\displaystyle }{\mathrm{Lim}}_{{\mathrm{NH}}_{\mathrm{4}}}={\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{k_{\mathrm{2}}+\left[{\mathrm{NH}}_{\mathrm{4}}\right]}}.\end{array}$$

Iron limitation follows a similar Michaelis–Menton form but with an additional term to enforce saturation at a critical threshold value:

$$\begin{array}{}\text{(A21)}& {\mathrm{Lim}}_{\mathrm{Fe}}=\min \left(\mathrm{1.0},\mathit{ϵ}+{\displaystyle \frac{\left[\mathrm{Fe}\right]}{k_{\mathrm{Fe}}+\left[\mathrm{Fe}\right]}}\cdot {\displaystyle \frac{k_{\mathrm{Fe}}+F{e}_{\mathrm{Crit}}}{F{e}_{\mathrm{Crit}}}}\right).\end{array}$$

Nitrate uptake is controlled by the minimum limitation factor between light, nitrate, and iron, while ammonium uptake is limited by light or ammonium:

$$\begin{array}{}\text{(A22)}& {\displaystyle }& {\displaystyle }\mathrm{Gpp}{\left({\mathrm{NO}}_{\mathrm{3}},X\right)}_k=P_{\max }\cdot \left[X\right]\cdot \min \left({\mathrm{Lim}}_{{\mathrm{NO}}_{\mathrm{3}}},{\mathrm{Lim}}_{\mathrm{Fe}},{\mathrm{Lim}}_I\right)\text{(A23)}& {\displaystyle }& {\displaystyle }\mathrm{Gpp}{\left({\mathrm{NH}}_{\mathrm{4}},X\right)}_k=P_{\max }\cdot \left[X\right]\cdot \min \left({\mathrm{Lim}}_{{\mathrm{NH}}_{\mathrm{4}}},{\mathrm{Lim}}_I\right).\end{array}$$

Primary production also occurs in the ice layer when ice is present. In the ice layer, production is a function of light, nutrient limitation, brine salinity, and temperature, following Jin et al. (2006). Light limitation uses the following photosynthesis–irradiance curve; unlike the pelagic production, this one includes strong photoinhibition at higher light levels:

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}{\mathrm{Lim}}_{I_{\mathrm{ice}}}& =\left(\mathrm{1}-\exp \left(-{\mathit{\alpha }}_{I_{\mathrm{b}}}\cdot {\displaystyle \frac{I_{\mathrm{0}}\cdot f_{\mathrm{PAR}}}{c_I}}\right)\right)\\ & \cdot \exp \left(-{\mathit{\beta }}_I\cdot {\displaystyle \frac{I_{\mathrm{0}}\cdot f_{\mathrm{PAR}}}{c_I}}\right),\end{array}\end{array}$$

where $\frac{I_{\mathrm{0}}\cdot f_{\mathrm{PAR}}}{c_I}$ is the photosynthetically active radiation converted to W m⁻².

As in the water column, nitrate limitation uses Michaelis–Menten uptake dynamics with ammonium inhibition (with f_r denoting the f ratio between new (nitrate) and regenerated (ammonium) production):

$$\begin{array}{}\text{(A25)}& {\displaystyle }\begin{array}{rl}{\mathrm{Lim}}_{N_{\mathrm{ice}}}& ={\displaystyle \frac{\left[\mathrm{IceNO}\mathrm{3}\right]}{k_{\mathrm{1}}+\left[\mathrm{IceNO}\mathrm{3}\right]}}\cdot \exp (-\mathit{\psi }\cdot [\mathrm{IceNH}\mathrm{4}\left]\right)\\ & +{\displaystyle \frac{\left[\mathrm{IceNH}\mathrm{4}\right]}{k_{\mathrm{2}}+\left[\mathrm{IceNH}\mathrm{4}\right]}}\end{array}\text{(A26)}& {\displaystyle }{f}_{\mathrm{r}}={\displaystyle \frac{\frac{\left[\mathrm{IceNO}\mathrm{3}\right]}{k_{\mathrm{1}}+\left[\mathrm{IceNO}\mathrm{3}\right]}\cdot \exp (-\mathit{\psi }\cdot [\mathrm{IceNH}\mathrm{4}\left]\right)}{{\mathrm{Lim}}_{N_{\mathrm{ice}}}}}.\end{array}$$

Brine salinity (S_b) is not tracked explicitly by the ice model, so instead it is estimated based on a piecewise polynomial fit to ice temperature (T_i, tracked by the ice model), following Arrigo et al. (1993):

$$\begin{array}{}\text{(A27)}& S_{\mathrm{b}}=c_{\mathrm{0}}+c_{\mathrm{1}}\cdot T_{\mathrm{i}}+c_{\mathrm{2}}\cdot T_{\mathrm{i}}^{\mathrm{2}}+c_{\mathrm{3}}\cdot T_{\mathrm{i}}^{\mathrm{3}}.\end{array}$$

	c0	c1	c2	c3
$T_{\mathrm{i}}\ge -\mathrm{22.9}$	−3.9921	−22.7	−1.0015	−0.019956
$-\mathrm{44.0}<T_{\mathrm{i}}<-\mathrm{22.9}$	206.24	−1.8907	−0.060868	−0.0010247
$T_{\mathrm{i}}\le -\mathrm{44.0}$	−4442.1	−277.86	−5.501	−0.03669

The effect of salinity on ice algae growth rate is also a polynomial fit (Arrigo and Sullivan, 1992):

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}{\mathit{\xi }}_{s_{\mathrm{b}}}=\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{2}}& +\mathrm{3.012}\times {\mathrm{10}}^{-\mathrm{2}}\cdot S_{\mathrm{b}}\\ & +\mathrm{1.0342}\times {\mathrm{10}}^{-\mathrm{3}}\cdot S_{\mathrm{b}}^{\mathrm{2}}\\ & -\mathrm{4.6033}\times {\mathrm{10}}^{-\mathrm{5}}\cdot S_{\mathrm{b}}^{\mathrm{3}}\\ & +\mathrm{4.926}\times {\mathrm{10}}^{-\mathrm{7}}\cdot S_{\mathrm{b}}^{\mathrm{4}}\\ & -\mathrm{1.659}\times {\mathrm{10}}^{-\mathrm{9}}\cdot S_{\mathrm{b}}^{\mathrm{5}}.\end{array}\end{array}$$

When running in climatological ice mode (CLIM_ICE_1D), where no explicit ice temperature is modeled, ${\mathit{\xi }}_{s_{\mathrm{b}}}=\mathrm{1.0}$.

The final primary production calculation is then

$$\begin{array}{}\text{(A29)}& {\displaystyle }\begin{array}{rl}\mathrm{Gpp}(\mathrm{IceNO}\mathrm{3},\mathrm{IcePhL}{)}_{\mathrm{ice}}=& {\mathit{\mu }}_{\mathrm{0}}\cdot \exp (\mathrm{0.0633}\cdot T_{k=N})\cdot {\mathit{\xi }}_{s_{\mathrm{b}}}\\ & \cdot \min \left({\mathrm{Lim}}_{I_{\mathrm{ice}}},{\mathrm{Lim}}_{N_{\mathrm{ice}}}\right)\\ & \cdot \left[\mathrm{IcePhL}\right]\cdot f_{\mathrm{r}}\end{array}\text{(A30)}& {\displaystyle }\begin{array}{rl}\mathrm{Gpp}(\mathrm{IceNH}\mathrm{4},\mathrm{IcePhL}{)}_{\mathrm{ice}}=& {\mathit{\mu }}_{\mathrm{0}}\cdot \exp (\mathrm{0.0633}\cdot T_{k=N})\cdot {\mathit{\xi }}_{s_{\mathrm{b}}}\\ & \cdot \min \left({\mathrm{Lim}}_{I_{\mathrm{ice}}},{\mathrm{Lim}}_{N_{\mathrm{ice}}}\right)\\ & \cdot \left[\mathrm{IcePhL}\right]\cdot \left(\mathrm{1}-f_{\mathrm{r}}\right).\end{array}\end{array}$$

In this case, T_k=N is the temperature of the top water layer, used to approximate the temperature of the ice itself.

A3.3 Grazing and predation

Pelagic grazing and predation fluxes are a function of a grazer's or predator's maximum ingestion rate (e_Y), its total prey availability, prey-specific feeding preferences (fp_XY; see Table A7), and the water temperature, using the multiple resource Holling type 3 functional response of Ryabchenko et al. (1997):

$$\begin{array}{}\text{(A31)}& \begin{array}{rl}\mathrm{Gra}(X,Y{)}_k& =Q_Y^{\left(\frac{T-Q_{TY}}{\mathrm{10}}\right)}\cdot e_Y\cdot \left[Y\right]\\ & \cdot {\displaystyle \frac{{\mathrm{fp}}_{XY}\cdot [X{]}^{\mathrm{2}}}{f_Y+\sum _Z\left({\mathrm{fp}}_{ZY}\cdot [Z{]}^{\mathrm{2}}\right)}},\end{array}\end{array}$$

where Y refers to the predator group, X is a specific prey group, and Z refers to the set of all prey groups of that predator. Note that some of the pelagic groups can graze on ice algae; when preyed upon, the ice algae concentration is adjusted as though it were located in the surface layer of the water:

$$\begin{array}{}\text{(A32)}& [\mathrm{IcePhL}{]}_k=\left\{\begin{array}{ll}[\mathrm{IcePhL}{]}_{\mathrm{ice}}\cdot \frac{h_{s_{\mathrm{ice}}}}{h_k},& k=N\\ \mathrm{0},& \mathrm{otherwise}.\end{array}\right.\end{array}$$

The maximum ingestion rates, e_Y, are constant for all groups except large-bodied copepods (NCaS and NCaO). These groups can be parameterized to perform seasonal diapause, and during periods of downward migration, their ingestion rates are dropped to e_Y=0 d⁻¹. See Sect. A4 for a description of the diapause time-of-year calculation.

Benthic processes in BESTNPZ are based on a greatly simplified version of the European Regional Seas Ecosystem model (ERSEM) (Ebenhöh et al., 1995). Benthic infauna graze on pelagic detritus and phytoplankton located within a certain distance of the bottom (currently hardcoded to dw=1.0 m). The feeding fluxes are defined as follows:

$$\begin{array}{}\text{(A33)}& {\displaystyle }& {\displaystyle }[X{]}_{\mathrm{ben}}=\underset{-h}{\overset{-h+dw}{\int }}[X]\mathrm{d}z\text{(A34)}& {\displaystyle }& {\displaystyle }{F}_X={\displaystyle \frac{{\left({\mathrm{fp}}_{X_B}\cdot [X{]}_{\mathrm{ben}}\right)}^{\mathrm{2}}}{{\mathrm{fp}}_{X_B}\cdot [X{]}_{\mathrm{ben}}+L_P}}\text{(A35)}& {\displaystyle }& {\displaystyle }\mathrm{Gra}(X,\mathrm{Ben})=Q_B^{\left(\frac{T-Q_{T_B}}{\mathrm{10}}\right)}\cdot e_{\mathrm{Ben}}\cdot \left[\mathrm{Ben}\right]\cdot {\displaystyle \frac{F_X}{\sum _Z{F}_Z+f_{P_B}}}.\end{array}$$

As in the pelagic grazing equation, X refers to a single pelagic prey group (PhS, PhL, Det, or DetF), and Z refers to the full set of these four groups. A weight factor, w_k,X, is calculated to distribute these losses proportionately across the water column (see Sect. A3.9):

$$\begin{array}{}\text{(A36)}& w_{k,X}={\displaystyle \frac{\underset{z_{\mathrm{lo},k}}{\overset{\min \left(z_{\mathrm{hi},k},-h+dw\right)}{\int }}\left[X\right]\mathrm{d}z}{[X{]}_{\mathrm{ben}}}},\end{array}$$

where z_lo,k and z_hi,k are the lower and upper depth limits of layer k.

Benthic infauna also graze on benthic detritus, following the same equation but with different parameters for prey threshold and half-saturation values:

$$\begin{array}{}\text{(A37)}& {\displaystyle }& {\displaystyle }{F}_X={\displaystyle \frac{{\left({\mathrm{fp}}_{X_B}\cdot \left[X\right]\right)}^{\mathrm{2}}}{{\mathrm{fp}}_{X_B}\cdot \left[X\right]+L_D}}\text{(A38)}& {\displaystyle }& {\displaystyle }\mathrm{Gra}(X,\mathrm{Ben})=Q_B^{\left(\frac{T-Q_{T_B}}{\mathrm{10}}\right)}\cdot e_{\mathrm{Ben}}\cdot \left[\mathrm{Ben}\right]\cdot {\displaystyle \frac{F_X}{F_X+f_{D_B}}}.\end{array}$$

Here, X refers to a single prey group, DetBen.

Note that the water column grazing fluxes are in units of mg C m⁻³ d⁻¹, while the benthic feeding fluxes are in mg C m⁻² d⁻¹.

A3.4 Egestion and excretion

Egestion fluxes associated with grazing and predation in the water column are a simple fraction of total prey eaten:

$$\begin{array}{}\text{(A39)}& \mathrm{Ege}(Y,\mathrm{Det}[\mathrm{F}]{)}_k=\left(\mathrm{1}-{\mathit{\gamma }}_Y\right)\cdot \sum _Z\mathrm{Gra}(Z,Y{)}_k.\end{array}$$

Egestion fluxes from the microzooplankton group (MZL) go to the slow-sinking detrital pool (Det); all other egestion fluxes go to the fast-sinking detrital pool (DetF).

Infauna egestion and excretion is a bit more complex; it is proportional to the prey eaten, with differing rates for detrital verus phytoplankton prey. The flux is split evenly, with half going to benthic detritus (DetBen) and half to NH₄.

$$\begin{array}{}\text{(A40)}& {\displaystyle }\begin{array}{rl}\mathrm{Exc}(\mathrm{Ben},\mathrm{DetBen})& =\mathrm{0.5}\cdot \left({\mathrm{ex}}_D\cdot \sum _{X=\det }\mathrm{Gra}(X,\mathrm{Ben})\right.\\ & \left.+{\mathrm{ex}}_P\cdot \sum _{X=\mathrm{phyto}}\mathrm{Gra}(X,\mathrm{Ben})\right)\end{array}\text{(A41)}& {\displaystyle }\begin{array}{rl}\mathrm{Exc}\left(\mathrm{Ben},{\mathrm{NH}}_{\mathrm{4}}\right)& =\mathrm{0.5}\cdot \left({\mathrm{ex}}_D\cdot \sum _{X=\det }\mathrm{Gra}(X,\mathrm{Ben})\right.\\ & \left.+{\mathrm{ex}}_P\cdot \sum _{X=\mathrm{phyto}}\mathrm{Gra}(X,\mathrm{Ben})\right)\end{array}\end{array}$$

As with benthic grazing, these benthic excretion fluxes are in units of mg C m⁻² d⁻¹. The flux to NH₄ is assumed to return to the bottom water column layer and is converted to a volumetric flux based on the thickness of that layer (see Sect. A3.9).

A3.5 Respiration

All pelagic producers and consumers except jellyfish respire following the temperature-dependent formulation of Arhonditsis and Brett (2005). The phytoplankton and microzooplankton groups maintain a constant basal metabolic (bm) rate:

$$\begin{array}{}\text{(A42)}& \mathrm{Res}{\left(X,{\mathrm{NH}}_{\mathrm{4}}\right)}_k=\exp \left(k_{t_{\mathrm{b}}}\cdot \left(T-T_{\mathrm{ref}}\right)\right)\cdot \mathrm{bm}\cdot \left[X\right].\end{array}$$

The larger zooplankton groups substitute a basal metabolic rate that includes a starvation response when prey is scarce:

$$\begin{array}{}\text{(A43)}& \mathrm{Res}{\left(X,{\mathrm{NH}}_{\mathrm{4}}\right)}_k=\exp \left(k_{t_{\mathrm{b}}}\cdot \left(T-T_{\mathrm{ref}}\right)\right)\cdot B_{\mathrm{met}}\cdot \left[X\right],\end{array}$$

where

$$\begin{array}{}\text{(A44)}& B_{\mathrm{met}}=\left\{\begin{array}{ll}\mathrm{bm}\cdot \left(\frac{\sum _Z\left({\mathrm{fp}}_{ZY}\cdot [Z{]}^{\mathrm{2}}\right)}{\mathrm{0.01}}\right)& \sum _Z\left({\mathrm{fp}}_{ZY}\cdot [Z{]}^{\mathrm{2}}\right)<\mathrm{0.01}\\ \mathrm{bm}& \mathrm{otherwise}.\end{array}\right.\end{array}$$

(The summation relates to the total available prey; see Sect. A3.3 for details.)

The large copepod groups (NCaS and NCaO) also include a diapause adjustment, such that their bm rate is reduced to 10 % of the bm parameter value during periods of downward migration. See Sect. A4 for details of the time-of-year calculation for diapause.

Jellyfish respiration also follows a temperature-dependent formula, after Uye and Shimauchi (2005):

$$\begin{array}{}\text{(A45)}& \mathrm{Res}{\left(\mathrm{Jel},{\mathrm{NH}}_{\mathrm{4}}\right)}_k=Q_{\mathrm{r}}^{\left(\frac{T-Q_{T_{\mathrm{r}}}}{\mathrm{10}}\right)}\cdot \mathrm{bm}\cdot \left[\mathrm{Jel}\right].\end{array}$$

Infaunal respiration includes terms for both basal metabolism and active metabolism proportional to grazing:

$$\begin{array}{}\text{(A46)}& {\displaystyle }\begin{array}{rl}\mathrm{Res}\left(\mathrm{Ben},{\mathrm{NH}}_{\mathrm{4}}\right)& =Q_B^{\left(\frac{T-Q_{T_B}}{\mathrm{10}}\right)}\cdot \mathrm{bm}\cdot \left[\mathrm{Ben}\right]\\ & +\mathrm{am}\cdot \left(\left(\mathrm{1}-{\mathrm{ex}}_D\right)\cdot \sum _{X=\det }\mathrm{Gra}(X,\mathrm{Ben})\right.\\ & +\left.\left(\mathrm{1}-{\mathrm{ex}}_P\right)\cdot \sum _{X=\mathrm{phyto}}\mathrm{Gra}(X,\mathrm{Ben})\right).\end{array}\end{array}$$

Finally, ice algae respiration uses a metabolic rate linearly proportional to its maximum growth rate, after Jin et al. (2006):

$$\begin{array}{}\text{(A47)}& \begin{array}{rl}\mathrm{Res}(\mathrm{IcePhl},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}& =r\cdot {\mathit{\mu }}_{\mathrm{0}}\cdot \exp \left(\mathrm{0.0633}\cdot T_{k=N}\right)\\ & \cdot \left[\mathrm{IcePhL}\right].\end{array}\end{array}$$

As with ice algae production, the surface water temperature is used as a proxy for ice temperature when calculating the temperature component of this rate.

A3.6 Mortality and senescence

Non-predatory mortality losses for phytoplankton groups are formulated as a linear closure term:

$$\begin{array}{}\text{(A48)}& \mathrm{Mor}(X,\mathrm{Det}{)}_k=m_L\cdot [X].\end{array}$$

Microzooplankton losses have the option of following either a linear closure as above (MZLM0LIN flag defined) or a quadratic closure:

$$\begin{array}{}\text{(A49)}& \mathrm{Mor}(X,\mathrm{Det}{)}_k=m_Q\cdot [X{]}^{\mathrm{2}}.\end{array}$$

Note that when switching between the linear and quadratic formulations, the relevant input parameter for the MZL group switches between mMZL and mpredMZL.

All larger zooplankton groups use a temperature-mediated quadratic closure term:

$$\begin{array}{}\text{(A50)}& \mathrm{Mor}(X,\mathrm{DetF}{)}_k=Q_Y^{\left(\frac{T-Q_{T_Y}}{\mathrm{10}}\right)}\cdot m_Q\cdot [X{]}^{\mathrm{2}}.\end{array}$$

Non-predatory mortality fluxes from the phytoplankton and microzooplankton groups go to the slow-sinking detritus, while all other non-predatory mortality losses go to the fast-sinking detritus.

The benthic infauna group includes both a linear and quadratic mortality function: the former to represent senescence and the latter as a predation closure term.

$$\begin{array}{}\text{(A51)}& \begin{array}{rl}\mathrm{Mor}(\mathrm{Ben},\mathrm{DetBen})& =Q_B^{\left(\frac{T-Q_{T_B}}{\mathrm{10}}\right)}\cdot \left(m_L\cdot \left[\mathrm{Ben}\right]\right.\\ \\ & \left.+m_Q\cdot [\mathrm{Ben}{]}^{\mathrm{2}}\right)\end{array}\end{array}$$

Finally, ice algae use a linear mortality rate with a temperature dependence, following Jin et al. (2006):

$$\begin{array}{}\text{(A52)}& \mathrm{Mor}(\mathrm{IcePhl},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}=\exp \left(\mathrm{rg}\cdot T_{k=N}\right)\cdot m_{L_{\mathrm{0}}}\cdot [\mathrm{IcePhL}].\end{array}$$

A3.7 Remineralization and nitrification

Detrital remineralization is proportional to the temperature and the nitrogen content of detritus, after Kawamiya et al. (2000):

$$\begin{array}{}\text{(A53)}& \mathrm{Rem}{\left(X,{\mathrm{NH}}_{\mathrm{4}}\right)}_k=\left(P_{v_{\mathrm{0}}}\cdot \exp \left(P_{v_T}\cdot T\right)\cdot \left[X\right]\cdot \mathit{\xi }\right)/\mathit{\xi }.\end{array}$$

The conversion from nitrogen content back to carbon content is done to maintain unit consistency with the other between-group fluxes, using the assumption that all living and detrital groups maintain identical C : N : Fe stoichiometry. Both fast- and slow-sinking detritus use the same parameters for this process.

Nitrification rate in the water column is also influenced by temperature (Arhonditsis and Brett, 2005):

$$\begin{array}{}\text{(A54)}& \begin{array}{rl}\mathrm{Nit}{\left({\mathrm{NH}}_{\mathrm{4}},{\mathrm{NO}}_{\mathrm{3}}\right)}_k& =\left(n_{\mathrm{0}}\cdot \exp \left(-k_{\mathrm{tntr}}\cdot {\left(T-T_{\mathrm{opt}}\right)}^{\mathrm{2}}\right)\right.\\ & \left.\cdot \left[{\mathrm{NH}}_{\mathrm{4}}\right]\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{k_{\mathrm{Nit}}+\left[{\mathrm{NH}}_{\mathrm{4}}\right]}}\right)/\mathit{\xi }.\end{array}\end{array}$$

Nitrification in the ice is a simple linear function of ammonium concentration (Jin et al., 2006):

$$\begin{array}{}\text{(A55)}& \mathrm{Nit}(\mathrm{IceNH}\mathrm{4},\mathrm{IceNO}\mathrm{3}{)}_{\mathrm{ice}}=\left(N_{\mathrm{nit}}\cdot \left[\mathrm{IceNH}\mathrm{4}\right]\right)/\mathit{\xi }.\end{array}$$

As with remineralization, the final nitrification rate values are converted to carbon units simply for bookkeeping purposes; they will be converted back to nitrogen units when used in the final rate-of-change equations (see Sect. A3.9.)

A3.8 Ice interface convective exchange

As an ice sheet grows, dense brine is released from the skeletal ice layer at its base and replaced with seawater; this results in a convective exchange of water and nutrients between the ice and surface water. The BESTNPZ model follows Jin et al. (2006) and sets the rate of this exchange using a polynomial function of ice growth rate:

$$\begin{array}{}\text{(A56)}& \begin{array}{rl}& T_{\mathrm{wi}}=\\ & \left\{\begin{array}{ll}\mathrm{720}\cdot \mathrm{86}\mathrm{400}\cdot \left(\mathrm{4.9}\times {\mathrm{10}}^{-\mathrm{6}}\cdot \left(-\frac{\mathrm{d}H}{\mathrm{d}t}\right)-\mathrm{1.39}\times {\mathrm{10}}^{-\mathrm{5}}\cdot {\left(-\frac{\mathrm{d}H}{\mathrm{d}t}\right)}^{\mathrm{2}}\right)& \\ \frac{\mathrm{d}H}{\mathrm{d}t}\le \mathrm{0}\\ \mathrm{72}\cdot \mathrm{86}\mathrm{400}\cdot \left(\mathrm{9.667}\times {\mathrm{10}}^{-\mathrm{11}}+\mathrm{4.49}\times {\mathrm{10}}^{-\mathrm{6}}\cdot \left(\frac{\mathrm{d}H}{\mathrm{d}t}\right)-\mathrm{1.39}\times {\mathrm{10}}^{-\mathrm{5}}\cdot {\left(\frac{\mathrm{d}H}{\mathrm{d}t}\right)}^{\mathrm{2}}\right)& \\ \frac{\mathrm{d}H}{\mathrm{d}t}>\mathrm{0},\end{array}\right.\end{array}\end{array}$$

where $\frac{\mathrm{d}H}{\mathrm{d}t}$ is the rate of change of ice thickness (m s⁻¹) between the current time step and the previous one. The resulting exchange rate, T_wi, is expressed in m d⁻¹.

The exchange in nutrients then becomes a function of the difference in concentrations in the surface layer of water versus the ice layer. Phytoplankton can be washed out of the skeletal ice layer but not in, so the exchange of ice algae and large phytoplankton assumes a concentration of 0 in the surface water:

$$\begin{array}{}\text{(A57)}& {\displaystyle }& {\displaystyle }\mathrm{Twi}\left(\mathrm{IceNO}\mathrm{3},{\mathrm{NO}}_{\mathrm{3}}\right)=T_{\mathrm{wi}}\cdot \left(\left[\mathrm{IceNO}\mathrm{3}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}\right]\right)/\mathit{\xi }\text{(A58)}& {\displaystyle }& {\displaystyle }\mathrm{Twi}\left(\mathrm{IceNH}\mathrm{4},{\mathrm{NH}}_{\mathrm{4}}\right)=T_{\mathrm{wi}}\cdot \left(\left[\mathrm{IceNH}\mathrm{4}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}\right]\right)/\mathit{\xi }\text{(A59)}& {\displaystyle }& {\displaystyle }\mathrm{Twi}(\mathrm{IcePhL},\mathrm{PhL})=T_{\mathrm{wi}}\cdot \left[\mathrm{IcePhL}\right].\end{array}$$

This equation results in a rate of exchange of material across the boundary (mg C m⁻² d⁻¹) that can have either a positive value (net movement from ice to water) or a negative value (net movement from water to ice). As in previous sections, the nutrient transport is converted to carbon units here purely for bookkeeping purposes.

A3.9 Total rate of change

The total rate of change for each state variable due to source-minus-sink processes is calculated as a sum of the rates detailed in the previous sections (Sect. A3.2–A3.8). Recall that in the previous sections, all flux rates taking place within the pelagic water column layers, or within the ice layer, were expressed in volumetric units (mg C m⁻³ d⁻¹), while all processes in the benthos or across the water–ice or water–benthos boundaries were expressed in per-area units (mg C m⁻² d⁻¹). Terms wrapped in square brackets apply only to the top layer (k=N), while terms wrapped in curly brackets apply only to the bottom layer (k=1).

$$\begin{array}{}\text{(A60)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{{\mathrm{NO}}_{\mathrm{3}}}_k& =\left(\mathrm{Nit}{\left({\mathrm{NH}}_{\mathrm{4}},{\mathrm{NO}}_{\mathrm{3}}\right)}_k-\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL})\end{array}}\mathrm{Gpp}{\left({\mathrm{NO}}_{\mathrm{3}},X\right)}_k\right.\\ & \left.+\left[{\displaystyle \frac{\mathrm{Twi}\left(\mathrm{IceNO}\mathrm{3},{\mathrm{NO}}_{\mathrm{3}}\right)}{h_k}}\right]\right)\cdot \mathit{\xi }\end{array}\end{array}$$

$$\begin{array}{}\text{(A61)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{{\mathrm{NH}}_{\mathrm{4}}}_k& =\left(\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\\ \mathrm{Cop},\mathrm{NCaS},\\ \mathrm{NCaO},\mathrm{EupS},\\ \mathrm{EupO},\mathrm{Jel})\end{array}}\mathrm{Res}{\left(X,{\mathrm{NH}}_{\mathrm{4}}\right)}_k\right.\\ & +\sum _{\begin{array}{c}X\in (\mathrm{Det},\\ \mathrm{DetF})\end{array}}\mathrm{Rem}{\left(X,{\mathrm{NH}}_{\mathrm{4}}\right)}_k\\ & -\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL})\end{array}}\mathrm{Gpp}{\left({\mathrm{NH}}_{\mathrm{4}},X\right)}_k\\ & -\mathrm{Nit}{\left({\mathrm{NH}}_{\mathrm{4}},{\mathrm{NO}}_{\mathrm{3}}\right)}_k+\left\{{\displaystyle \frac{\mathrm{Exc}{\left(\mathrm{Ben},{\mathrm{NH}}_{\mathrm{4}}\right)}_k}{h_k}}\right\}\\ & +\left\{{\displaystyle \frac{\mathrm{Res}{\left(\mathrm{Ben},{\mathrm{NH}}_{\mathrm{4}}\right)}_k}{h_k}}\right\}\\ & \left.+\left[{\displaystyle \frac{\mathrm{Twi}{\left(\mathrm{IceNH}\mathrm{4},{\mathrm{NH}}_{\mathrm{4}}\right)}_k}{h_k}}\right]\right)\cdot \mathit{\xi }\end{array}\end{array}$$

$$\begin{array}{}\text{(A62)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{PhS}}_k& =\sum _{\begin{array}{c}X\in ({\mathrm{NO}}_{\mathrm{3}},\\ {\mathrm{NH}}_{\mathrm{4}})\end{array}}\mathrm{Gpp}(X,\mathrm{PhS}{)}_k\\ & -\sum _{\begin{array}{c}X\in (\mathrm{MZL},\\ \mathrm{Cop},\mathrm{NCaS},\\ \mathrm{NCaO},\mathrm{EupS},\\ \mathrm{EupO})\end{array}}\mathrm{Gra}(\mathrm{PhS},X{)}_k-\mathrm{Mor}(\mathrm{PhS},\mathrm{Det}{)}_k\\ & -\mathrm{Res}{\left(\mathrm{PhS},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\\ & -{\displaystyle \frac{\mathrm{Gra}(\mathrm{PhS},\mathrm{Ben})\cdot w_{k,\mathrm{PhS}}}{h_k}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A63)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{PhL}}_k& =\sum _{\begin{array}{c}X\in ({\mathrm{NO}}_{\mathrm{3}},\\ {\mathrm{NH}}_{\mathrm{4}})\end{array}}\mathrm{Gpp}(\mathrm{X},\mathrm{PhL}{)}_k\\ & -\sum _{\begin{array}{c}X\in (\mathrm{MZL},\\ \mathrm{Cop},\mathrm{NCaS},\\ \mathrm{NCaO},\mathrm{EupS},\\ \mathrm{EupO})\end{array}}\mathrm{Gra}(\mathrm{PhL},X{)}_k-\mathrm{Mor}(\mathrm{PhL},\mathrm{Det}{)}_k\\ & -\mathrm{Res}{\left(\mathrm{PhL},{\mathrm{NH}}_{\mathrm{4}}\right)}_k-{\displaystyle \frac{\mathrm{Gra}(\mathrm{PhL},\mathrm{Ben})\cdot w_{k,\mathrm{PhL}}}{h_k}}\\ & +\left[{\displaystyle \frac{\mathrm{Twi}(\mathrm{IcePhL},\mathrm{PhL})}{h_k}}\right]\end{array}\end{array}$$

$$\begin{array}{}\text{(A64)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{MZL}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL})\end{array}}\mathrm{Gra}(X,\mathrm{MZL}{)}_k-\mathrm{Ege}(\mathrm{MZL},\mathrm{Det}{)}_k\\ & -\sum _{\begin{array}{c}X\in (\mathrm{Cop},\\ \mathrm{NCaS},\mathrm{NCaO},\\ \mathrm{EupS},\mathrm{EupO})\end{array}}\mathrm{Gra}(\mathrm{MZL},X{)}_k\\ & -\mathrm{Mor}(\mathrm{MZL},\mathrm{Det}{)}_k-\mathrm{Res}{\left(\mathrm{MZL},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A65)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{Cop}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\left[\mathrm{IcePhL}\right])\end{array}}\mathrm{Gra}(X,\mathrm{Cop}{)}_k\\ & -\mathrm{Ege}(\mathrm{Cop},\mathrm{DetF}{)}_k-\sum _{\begin{array}{c}X\in (\mathrm{EupS},\\ \mathrm{EupO},\mathrm{Jel})\end{array}}\mathrm{Gra}(\mathrm{Cop},X{)}_k\\ & -\mathrm{Mor}(\mathrm{Cop},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{Cop},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A66)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{NCaS}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\\ \left[\mathrm{IcePhL}\right])\end{array}}\mathrm{Gra}(X,\mathrm{NCaS}{)}_k\\ & -\mathrm{Ege}(\mathrm{NCaS},\mathrm{DetF}{)}_k-\mathrm{Gra}(\mathrm{NCaS},\mathrm{Jel}{)}_k\\ & -\mathrm{Mor}(\mathrm{NCaS},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{NCaS},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A67)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{EupS}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\\ \mathrm{Cop},\left[\mathrm{IcePhL}\right])\end{array}}\mathrm{Gra}(X,\mathrm{EupS}{)}_k\\ & -\mathrm{Ege}(\mathrm{EupS},\mathrm{DetF}{)}_k-\mathrm{Gra}(\mathrm{EupS},\mathrm{Jel}{)}_k\\ & -\mathrm{Mor}(\mathrm{EupS},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{EupS},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A68)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{NCaO}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\\ \left[\mathrm{IcePhL}\right])\end{array}}\mathrm{Gra}(X,\mathrm{NCaO}{)}_k\\ & -\mathrm{Ege}(\mathrm{NCaO},\mathrm{DetF}{)}_k-\mathrm{Gra}(\mathrm{NCaO},\mathrm{Jel}{)}_k\\ & -\mathrm{Mor}(\mathrm{NCaO},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{NCaO},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A69)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{EupO}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL},\\ \mathrm{Cop},\left[\mathrm{IcePhL}\right])\end{array}}\mathrm{Gra}(X,\mathrm{EupO}{)}_k\\ & -\mathrm{Ege}(\mathrm{EupO},\mathrm{DetF}{)}_k-\mathrm{Gra}(\mathrm{EupO},\mathrm{Jel}{)}_k\\ & -\mathrm{Mor}(\mathrm{EupO},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{EupO},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A70)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{Det}}_k& =\mathrm{Ege}(\mathrm{MZL},\mathrm{Det}{)}_k+\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL},\mathrm{MZL})\end{array}}\mathrm{Mor}(X,\mathrm{Det}{)}_k\\ & -\sum _{X\in (\mathrm{EupS},\mathrm{EupO})}\mathrm{Gra}(\mathrm{Det},X{)}_k\\ & -\mathrm{Rem}{\left(\mathrm{Det},{\mathrm{NH}}_{\mathrm{4}}\right)}_k-{\displaystyle \frac{\mathrm{Gra}(\mathrm{Det},\mathrm{Ben})\cdot w_{k,\mathrm{Det}}}{h_k}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A71)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{DetF}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{Cop},\\ \mathrm{NCaS},\mathrm{NCaO},\\ \mathrm{EupS},\mathrm{EupO},\\ \mathrm{Jel})\end{array}}\mathrm{Ege}(X,\mathrm{DetF}{)}_k\\ & +\sum _{\begin{array}{c}X\in (\mathrm{Cop},\\ \mathrm{NCaS},\mathrm{NCaO},\\ \mathrm{EupS},\mathrm{EupO},\\ \mathrm{Jel})\end{array}}\mathrm{Mor}(X,\mathrm{DetF}{)}_k\\ & -\sum _{\begin{array}{c}X\in (\mathrm{EupS},\\ \mathrm{EupO})\end{array}}\mathrm{Gra}(\mathrm{DetF},X{)}_k-\mathrm{Rem}{\left(\mathrm{DetF},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\\ & -{\displaystyle \frac{\mathrm{Gra}(\mathrm{DetF},\mathrm{Ben})\cdot w_{k,\mathrm{DetF}}}{h_k}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A72)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{Jel}}_k& =\sum _{\begin{array}{c}X\in (\mathrm{Cop},\\ \mathrm{NCaS},\mathrm{NCaO},\\ \mathrm{EupS},\mathrm{EupO})\end{array}}\mathrm{Gra}(X,\mathrm{Jel}{)}_k-\mathrm{Ege}(\mathrm{Jel},\mathrm{DetF}{)}_k\\ & -\mathrm{Mor}(\mathrm{Jel},\mathrm{DetF}{)}_k-\mathrm{Res}{\left(\mathrm{Jel},{\mathrm{NH}}_{\mathrm{4}}\right)}_k\end{array}\end{array}$$

$$\begin{array}{}\text{(A73)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{Fe}}_k=\left(\sum _{\begin{array}{c}X\in (\mathrm{PhS},\\ \mathrm{PhL})\end{array}}\mathrm{Gpp}{\left({\mathrm{NO}}_{\mathrm{3}},X\right)}_k\right)\cdot \mathrm{Fe}C\end{array}$$

$$\begin{array}{}\text{(A74)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{Ben}=& \sum _{\begin{array}{c}X\in (\mathrm{Det},\\ \mathrm{DetF},\mathrm{PhS},\\ \mathrm{PhL},\mathrm{BenDet})\end{array}}\mathrm{Gra}(X,\mathrm{Ben})\\ & -\sum _{\begin{array}{c}X\in ({\mathrm{NH}}_{\mathrm{4}},\\ \mathrm{BenDet})\end{array}}\mathrm{Exc}(\mathrm{Ben},X)-\mathrm{Mor}(\mathrm{Ben},\mathrm{BenDet})\\ & -\mathrm{Res}\left(\mathrm{Ben},{\mathrm{NH}}_{\mathrm{4}}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(A75)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{BenDet}& =\mathrm{Exc}(\mathrm{Ben},\mathrm{BenDet})+\mathrm{Mor}(\mathrm{Ben},\mathrm{BenDet})\\ & -\mathrm{Gra}(\mathrm{BenDet},\mathrm{Ben})-\mathrm{Rem}\left(\mathrm{BenDet},{\mathrm{NH}}_{\mathrm{4}}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(A76)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{IcePhL}& =\sum _{\begin{array}{c}X\in (\mathrm{IceNO}\mathrm{3},\\ \mathrm{IceNH}\mathrm{4})\end{array}}\mathrm{Gpp}(X,\mathrm{IcePhL}{)}_{\mathrm{ice}}\\ & -\left(\sum _{\begin{array}{c}X\in (\mathrm{Cop},\\ \mathrm{NCaS},\mathrm{NCaO},\\ \mathrm{EupS},\mathrm{EupO})\end{array}}\mathrm{Gra}(\mathrm{IcePhL},X{)}_k\right)\cdot {\displaystyle \frac{h_k}{h_{s_{\mathrm{ice}}}}}\\ & -\mathrm{Mor}(\mathrm{IcePhL},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}\\ & -\mathrm{Res}(\mathrm{IcePhL},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}\\ & -{\displaystyle \frac{\mathrm{Twi}(\mathrm{IcePhL},\mathrm{PhL})}{h_{s_{\mathrm{ice}}}}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A77)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{IceNO}\mathrm{3}& =\left(\mathrm{Nit}(\mathrm{IceNH}\mathrm{4},\mathrm{IceNO}\mathrm{3}{)}_{s_{\mathrm{ice}}}\right.\\ & -\mathrm{Gpp}(\mathrm{IceNO}\mathrm{3},\mathrm{IcePhL}{)}_{\mathrm{ice}}\\ & \left.-{\displaystyle \frac{\mathrm{Twi}\left(\mathrm{IceNO}\mathrm{3},{\mathrm{NO}}_{\mathrm{3}}\right)}{h_{s_{\mathrm{ice}}}}}\right)\cdot \mathit{\xi }\end{array}\end{array}$$

$$\begin{array}{}\text{(A78)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{IceNH}\mathrm{4}& =\left(\mathrm{Res}(\mathrm{IcePhL},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}\right.\\ & +\mathrm{Mor}(\mathrm{IcePhL},\mathrm{IceNH}\mathrm{4}{)}_{\mathrm{ice}}\\ & -\mathrm{Gpp}(\mathrm{IceNH}\mathrm{4},\mathrm{IcePhL}{)}_{\mathrm{ice}}\\ & -\mathrm{Nit}(\mathrm{IceNH}\mathrm{4},\mathrm{IceNO}\mathrm{3}{)}_{\mathrm{ice}}\\ & \left.-{\displaystyle \frac{\mathrm{Twi}\left(\mathrm{IceNH}\mathrm{4},{\mathrm{NH}}_{\mathrm{4}}\right)}{h_{s_{\mathrm{ice}}}}}\right)\cdot \mathit{\xi }.\end{array}\end{array}$$

A4 Vertical movement and exchanges

All vertical movement in the BESTNPZ model is calculated using a piecewise parabolic method and weighted non-oscillatory scheme, following the sediment settling code from a ROMS sediment model (Warner et al., 2008). This scheme allows for fast sinking speeds that may cause material to cross multiple layers, without being constrained by the Courant–Friedrichs–Lewy (CFL) criterion.

We have modified this scheme slightly to allow a zero-flux boundary condition to be imposed at a specified depth. We also allow the use of a vertical velocity rather sinking speed; a negative velocity implies sinking, while a positive velocity indicates rising.

A4.1 Sinking of phytoplankton and detritus

Phytoplankton and detrital groups (PhS, PhL, Det, and DetF) are subject to vertical settling, with a constant sinking speed for each state variable.

When material crosses the water–benthic boundary, it is assumed that 20 % of the material becomes biologically unavailable (Walsh et al., 1981; Walsh and McRoy, 1986), and 1 % is lost to denitrification (personal communication with David Shull via Gibson and Spitz, 2011). The remaining 79 % of the flux across the boundary is transferred to the benthic detritus (DetBen). Note that the “denitrification” flux is not tracked explicitly but simply subtracted from the flux reaching benthic detritus; the biomass associated with both burial and denitrification is lost from the system.

A4.2 Copepod diapause

Copepod diapause is included in the BESTNPZ model by imposing seasonal movement through the water column on both large-bodied copepod groups (NCaS and NCaO), coupled with modifications to their feeding and respiration rates during periods of downward movement. The timing of copepod diapause is specified via input parameters for sinking start (S_start), sinking end (S_end), rising start (R_start), and rising end (R_end) day for each group, all specified as days of the year. These four parameters combine to define periods of downward-directed movement ($S_{\mathrm{start}}\le t_{\mathrm{doy}}\le S_{\mathrm{end}}$), upward-directed movement ($R_{\mathrm{start}}\le t_{\mathrm{doy}}\le R_{\mathrm{end}}$), and no directed movement ($R_{\mathrm{end}}<t_{\mathrm{doy}}<S_{\mathrm{start}}$ and $S_{\mathrm{end}}<t_{\mathrm{doy}}<R_{\mathrm{start}}$) (during all times, both groups are still subject to passive advection and diffusion).

Earlier versions allowed specification of S_start, S_end, R_start, and R_end for the offshore group (NCaO) only, with the onshore group automatically lagged by 30 d behind the offshore group. In the current version of the code, the 30 d lag is assumed when all four parameters are set to 0 for the onshore group (this maintains backward compatibility with older input files that do not include values for the onshore group; BESTNPZ sets missing input parameters to zero by default). Diapause can be turned off for either group by setting all four timing parameters to the same non-zero value.

Onshore copepods migrate to 200 m depth during their diapause period. During downward movement, a zero-flux condition is set at 200 m or at the bottom boundary, whichever is shallower. When migrating upward, a zero-flux condition is applied to the top of the surface layer.

Offshore copepods migrate to 400 m depth for their diapause period. A zero-flux condition is set at 400 m. In shallower waters, any biomass that crosses the bottom boundary is transferred to the benthic detritus (DetBen) group. During upward movement, a zero-flux condition is applied to the top of the surface layer.

A4.3 Euphausiid diel vertical migration

Diel vertical migration is currently implemented for on-shelf euphausiids through the (still experimental) EUPDIEL compilation flag. When defined, sinking and rising velocities are applied such that EupS moves at a hardcoded speed of 100 m d⁻¹ toward a target depth, defined as the shallowest depth layer where photosynthetically active radiation (PAR ⋅ I₀) is less than 0.5 E m⁻² d⁻¹. This option was turned off during the simulations detailed in this study.

A5 Analytical relaxation of state variables

Iron concentrations throughout the water column are initialized using a vertical profile that prescribes values above 50 m and below 300 m, with a linear interpolation between these two depths. The values of the shallow- and deep-water limits are a function of the bottom depth in a grid cell, with higher values in shallow water and lower values in deep water. The primary source of iron in this region is the sediment, and therefore this gradient between shallow and deep water is intended to capture the iron differences between on-shelf and off-shelf regions.

Iron-related biogeochemical processes are not included in this model. Instead, the iron state variable is continuously nudged towards these prescribed vertical profiles on an annual timescale. The nudging process is implemented using the generic ROMS framework for climatological nudging, with the TNUDG input parameter controlling the strength of the relaxation calculations.

Figure A1 Schematic of the BESTNPZ model. Circular nodes represent state variables (gold indicates nutrient, green indicates producer, blue indicates consumer, brown indicates detritus). Edges (lines) represent fluxes between state variables and curve clockwise from source node to sink node. Edge colors indicate process type: green indicates primary production, blue indicates grazing and predation, brown indicates egestion, gold indicates respiration, red indicates remineralization, pink indicates nitrification, orange indicates non-predatory mortality, tan indicates excretion, purple indicates convective exchange, gray indicates sinking to seafloor, and navy indicates freezing/melting of ice.

Table A1 Biological state variables in the BESTNPZ model.

Thimijan and Heins (1983)

Table A2 Notation and description for variables derived from the physical model.

Table A3 Notation and description for input parameters applicable to all biomass rate-of-change processes.

Table A4 Notation and description for input parameters related to light attenuation.

Table A5 Notation and description for input parameters related to gross primary production.

Table A6 Notation and description for input parameters related to grazing and predation.

Table A7 Feeding preferences matrix, indicating feeding preferences for each predator (columns) on each prey (rows).

Table A8 Notation and description for input parameters related to egestion and excretion.

Table A9 Notation and description for input parameters related to respiration. (See Table A6 for pelagic prey preferences and infauna Q₁₀ parameters, Table A8 for infauna excretion fraction parameters, and Table A5 for ice algae growth rate parameter.)

Table A10 Notation and description for input parameters related to non-predatory mortality. (See Table A6 for Q₁₀ parameters.)

Table A11 Notation and description for input parameters related to remineralization and nitrification.

Table A12 Notation and description for input parameters related to vertical movement.

Table A13 Notation and description for input parameters related to state variable nudging.

Appendix B: Light attenuation and limitation in BESTNPZ

Light attenuation in the Bering10K-BESTNPZ model has undergone a variety of modifications between the Hermann et al. (2013, 2016) publications, and between the Hermann et al. (2016) publication and this one. In this section, we clarify the motivations behind these changes and their effects on both biological and physical state variables in the fully coupled model.

Light attenuation formulations are used in a typical, physics-only implementation of ROMS in order to determine the attenuation of shortwave radiation in the water column and distribute surface heat fluxes appropriately. When a biological module is added, any light attenuation calculations necessary for photosynthesis (typically focused only on the photosynthetically active wavelengths) are coded separately, independent of the shortwave attenuation code in the physical core of ROMS, and with no direct feedback of the biology on the physics.

The Hermann et al. (2013) version of BESTNPZ followed this default setup; the attenuation coefficient of photosynthetically active radiation for the biological model followed Gibson and Spitz (2011):

$$\begin{array}{}\text{(B1)}& {\displaystyle }& {\displaystyle }{I}_{\mathrm{PAR}}\left(z\right)=I_{\mathrm{PAR}}\left(\mathrm{0}\right)\exp (-Kz)\text{(B2)}& {\displaystyle }& {\displaystyle }K=k_{\mathrm{ext}}+k_{\mathrm{A}}{C}^{\mathrm{0.0428}},\end{array}$$

where C is the concentration of chlorophyll a (mg Chl a m⁻³), k_ext is the clear-water attenuation coefficient, and the k_A term is attenuation due to light absorption by chlorophyll; the k_ext and k_A values, as well as the C-term exponent, are derived from an empirical fit of open-ocean chlorophyll concentrations versus measured attenuation (Morel, 1988).

The physical model relied on the ROMS default formula, which applies separate attenuation length scales to blue-green wavelengths and other wavelengths following Paulson and Simpson (1977):

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}I\left(z\right)& =I\left(\mathrm{0}\right)\left(\left(\mathrm{1}-a_{\mathrm{frac}}\right)\exp \left(-z/a_{{\mathit{\mu }}_{\mathrm{1}}}\right)\right.\\ & \left.+a_{\mathrm{frac}}\exp \left(-z/a_{{\mathit{\mu }}_{\mathrm{2}}}\right)\right).\end{array}\end{array}$$

The values for a_frac (unitless), $a_{{\mathit{\mu }}_{\mathrm{1}}}$ (m), and $a_{{\mathit{\mu }}_{\mathrm{2}}}$ (m) are set via an internal lookup table based on one of several water classification types. These water type classes account for the varying chlorophyll, sediment, and CDOM concentrations in different types of environments. Early versions of the Bering10K domain alternated between Class I (a_frac=0.58, $a_{{\mathit{\mu }}_{\mathrm{1}}}=\mathrm{0.35}$, $a_{{\mathit{\mu }}_{\mathrm{2}}}=\mathrm{23.0}$) and Class III (a_frac=0.78, $a_{{\mathit{\mu }}_{\mathrm{1}}}=\mathrm{1.4}$, $a_{{\mathit{\mu }}_{\mathrm{2}}}=\mathrm{7.9}$) parameters; the Hermann et al. (2013) study used the Class III option. Note that in this setup, the amount of chlorophyll in the water column assumed by the parameters of Eq. (B3) is independent from the chlorophyll levels modeled by the coupled biological model.

A more harmonious approach to light attenuation was adopted in the Hermann et al. (2016) version of the model. As one of several adjustments aimed at addressing a warm bias in heating of the water column, the attenuation equation in the physical model replaced the $a_{{\mathit{\mu }}_{\mathrm{2}}}$ attenuation length scale parameter with the PAR-related length scale (K⁻¹) from the biological model; this provided direct feedback of phytoplankton on heat absorption in the water column and better encompassed the varying light conditions expected between the low-productivity basin and high-productivity shelf region. A rough approximation of attenuation due to sediment was also added as a function of bottom depth (h) in both the biological and physical calculations for light attenuation:

$$\begin{array}{}\text{(B4)}& K=k_{\mathrm{ext}}+k_{\mathrm{A}}{C}^{\mathrm{0.0428}}+\mathrm{2.0}\exp \left(\mathrm{0.05}h\right).\end{array}$$

Hermann et al. (2016) also increased the value of the k_A parameter (from 0.0518 to 0.121) and adjusted the values of a_frac and $a_{{\mathit{\mu }}_{\mathrm{1}}}$ to Class I values (now hardcoded within the source code).

The new biological feedback resulted in only a very small change in water column heat content, and the model warm bias was later addressed through other adjustments to surface boundary conditions. However, the addition of the sediment attenuation term strongly affected the biological calculations. Even in the absence of grazing and under nutrient-replete, peak-irradiance conditions, the Hermann et al. (2016) set of parameters raises the compensation depths for both large and small phytoplankton to very shallow depths; in the shallowest waters of our domain, primary production is only possible in the upper 4–5 m of the water column (Fig. B1). While high sediment attenuation is possible in a few specific locations in the Bering10K domain, near the mouths of the Yukon, Anadyr, and Kuskokwim rivers during peak streamflow, there is no evidence to support such strong sediment shading throughout the entire inner domain.

During early validation of the circa Hermann et al. (2016) model for this study (prior to any adjustments for bug fixes in the code), it was also noted that light was the key limiting factor in the deep basin region. The sediment attenuation term is negligible in this deep-water region, and the compensation depth is reasonable in this location. However, the coarse vertical resolution of the model domain resulted in a surface layer nearly 40 m thick. The discretization used by this model assumes that all processes are approximately linear within a given layer, and that the layer midpoint can be therefore be used to approximate characteristics of the layer. In the case of light levels in these deep basin surface layers, this assumption did not hold, and the resulting light levels at the midpoints of the surface layers were too low to support any significant growth of phytoplankton.

For this study, all parameters related to light were revisited. A new sediment attenuation term was derived empirically from satellite inherent optical property measurements from the region (Fig. B2). This remains a rough estimate that does not consider the seasonal variability of sediment and detrital matter that contribute to this term but alleviates the previously excessive limitation in shallow water. Attenuation-related parameters for clear-water attenuation and chlorophyll-based attenuation were updated to values consistent with literature (k_ext=0.034 following Morel et al., 2007, and k_A=0.0518 and k_B=0.428 following Morel, 1988). It was noted that the value used for photosynthetically active radiation fraction (PAR_frac=0.5) in the Gibson and Spitz (2011) and Hermann et al. (2013, 2016) publications was higher than most field-based and analytical estimates. We adjusted this parameter to a value of 0.42 based on examination of satellite-derived PAR versus our model's surface boundary condition shortwave radiation values. Finally, the vertical resolution in the physical model was increased from 10 layers to 30 layers to allow for better resolution of mixed layer dynamics and light attenuation in the basin.

The updated light attenuation equations and parameters, coupled with the higher vertical resolution in the newer simulations, better reflect the true mechanisms controlling light levels in the Bering Sea; the bathymetry-following artifacts in phytoplankton spatial patterns are no longer present (Fig. B3). However, the lower light limitation in certain parts of the domain has exposed some previously overlooked deficiencies in micro- and macronutrient limitation in the BESTNPZ model. Observations suggest that the deep basin should be a low-production region, a pattern that is present but for the wrong reasons in the Hermann et al. (2016) version of the model and absent in the updated, more mechanistically sound version of the model. For some applications, it may be more useful to have the correct gradient between on- and off-shelf populations, even if it is for the incorrect reason, than to have mechanistically sound light limitation; however, users should be fully aware of this discrepancy before analyzing any biological output from either the earlier Hermann et al. (2016) simulations or this newer set of simulations using the updated code.

Figure B1 Compensation depth (depth where photosynthesis balances phytoplankton loss terms) and euphotic depth (depth of 1 % PAR) in three versions of the BESTNPZ model, under low (light-colored) and high (dark-colored) temperature conditions, assuming a surface irradiance of 200 W m⁻² and nutrient-replete waters. The “no Chl or grazing” scenario (blue) assumes losses due to respiration and non-predatory mortality only, with no attenuation due to phytoplankton itself. The “low Chl and grazing” scenario (green) assumes attenuation due to 0.5 mmol Chl a m⁻³ throughout the water column and adds a constant grazing loss rate of 0.05 d⁻¹. Gray lines indicate boundaries of the ROMS depth levels in the 10-layer (dark) and 30-layer (light) versions of the model. (Note that values are very nearly identical for small and large phytoplankton due to the proportional scaling of rate parameters between the groups.)

Figure B2 Attenuation due to sediment derives from a power law fit of satellite estimations of absorption due to gelbstoff and detrital material (entire mission composite from VIIRS, 2012–2018) versus bottom depth. Subpanels show the same data in linear and logarithmic space.

Figure B3 Mean depth-integrated chlorophyll in the Hermann et al. (2016) version of BESTNPZ versus this study. The contour line indicates the 2500 m isobath.

Appendix C: Coupling BESTNPZ to ROMS: compilation options and output variables

The following appendix provides instructions for setting up input for, compiling, and retrieving output from a simulation of ROMS coupled to the BESTNPZ biological model. At present, the BESTNPZ code has only been coupled to the Bering10K domain, and a few quirks of its internal coding prevent it from being coupled to other physical domains in its current form. Therefore, this user manual should be considered to be specific to the Bering10K-BESTNPZ ROMS model setup.

The primary input parameters specific to the BESTNPZ are described in Appendix A. The remaining setup options relate to the compilation flags and the indices of input and output state variables, described in the following sections.

C1 Compilation flags

Within the ROMS source code, C preprocessing flags are used to selectively compile the code. Several C preprocessing flags appear throughout the BESTNPZ code; those that currently exist in the master branch of the source code are detailed in the following table.

A few of these compilation options resulted from the gradual addition of new features (such as the addition of benthic and ice state variables) that later became, for all intents and purposes, permanent additions to the BESTNPZ model. These are noted in the table as “recommended to always define”. While we have attempted to keep the source code flexible to all compilation options, we rarely test for compilation stability with these undefined, and all recent model validation work assumes these model features are present.

C2 Output variables and input parameter indices

The BESTNPZ module allows for a large number of variables to be added to the various ROMS output files (history, average, and station files).

Regardless of compilation options, the biological tracer variables are saved and available as output variables. Water column tracers are specified via the H/Sout(idTvar) input parameters. The exact tracer variables available, and their positions within the idTvar array, depend on whether the ICE_BIO, BENTHIC, IRON_LIMIT, and JELLY compilation flags are defined; see Table C2 for details. Benthic tracers can be specified for output through the H/Sout(idBvar) input; iBen= = 1 and iDetBen = 2 within this array. The ice variables are turned on and off with the variable-specific input parameters H/Sout(idIcePhL), H/Sout(idIceNO3), and H/Sout(idIceNH4).

If the STATIONARY flag is defined, many more intermediate diagnostic variables are saved internally and available for output. These are controlled by the idTSvar index array to Hout and Sout. Note that when running with a subdivided time step (BioIter > 1), these diagnostic variables will reflect the values calculated during the final subdivision only.

Older versions of the code included a second compilation flag, STATIONARY2 (with the corresponding index array idTS2var), to define two-dimensional stationary diagnostics. More recent versions of the code no longer include 2-D diagnostic variables. However, the code structure is still in place for this if it becomes necessary in the future.

See Table C3 for a comprehensive description of all output variables associated with the BESTNPZ model.

Table C1 Compilation flags related to the BEST_NPZ model.

Table C2 Indices in the idTvar array. Note that Hout(idTvar) appears separately in the ocean.in file and BPARNAM file, with the latter including only biological active tracers. The Sout(idTvar) input appears only once, in the STANAME file, and biological tracer indices begin at NAT + 1; NAT is here assumed to be 2, for temperature and salinity. In all cases, note the skip between iPhL and iMZL due to the now-deprecated small microzooplankton group.

Table C3 Output variables available in the BEST_NPZ and FEAST models. The variable type indicates whether each variable is two- or three-dimensional, and whether they are located in the center of each grid cell (ρ), on the corners (ψ), or along an edge (u or v); currently, all biological variables are ρ variables.