2.1.1 NEMURO model description

The biogeochemical model for this study is based on NEMURO (Kishi et al., 2007) but has been modified and parameterized to more accurately reflect the ecology of the GoM (herein called NEMURO-GoM). NEMURO is a concentration-based, lower-trophic-level ecosystem model originally developed and parameterized for the North Pacific. Like most marine functional-group biogeochemical models, it is structured around simplified representations of the lower food web originating from earlier nutrient–phytoplankton–zooplankton models (Fasham et al., 1990; Franks, 2002; Riley, 1946). Complexity is added through additional state variables and transfer functions with the specific goal of resolving dynamics within the nutrient, phytoplankton and zooplankton pools. In total, NEMURO has 11 state variables: six nonliving state variables – nitrate (NO₃), ammonium (NH₄), dissolved organic nitrogen (DON), particulate organic nitrogen (PON), silicic acid (Si(OH)₄) and particulate silica (opal); two phytoplankton state variables – small (SP) and large phytoplankton (LP); and three zooplankton state variables – small (SZ), large (LZ) and predatory zooplankton (PZ).

Each biological state variable in NEMURO is an aggregated representation of taxonomically diverse plankton groups that function similarly in the ecosystem. The phytoplankton community is modeled as two functional types of obligate autotrophs: small phytoplankton (SP, predominantly cyanobacteria and picoeukaryotes in the GoM) and large phytoplankton (LP, diatoms). Small zooplankton (SZ) represent heterotrophic protists, and metazoan zooplankton are divided into suspension-feeding mesozooplankton (LZ) and predatory zooplankton (PZ). Here, we assume that LZ and PZ are nonmigratory. Heterotrophic bacteria are implicitly represented by temperature-dependent decomposition rates, which represent nitrification and remineralization processes. NEMURO uses nitrogen as a model currency since it is the major limiting macronutrient in much of the ocean. Silica is also included as a potentially colimiting nutrient for diatoms (i.e., LP).

By default, sinking in NEMURO is restricted to PON and opal, and benthic processes are not included. Here, because of the large shelf area in the GoM, we implemented a simple diagenesis of PON∕opal to NO₃∕SiO₄ and removal of PON∕opal through sedimentation, where 1 % of the flux sinking out of bottom cell was removed and 10 % converted back into NO₃∕SiO₄. However, we found that this had no significant impact on the simulated surface Chl or mesozooplankton biomass on the shelf. The inclusion of a more complex sediment diagenesis model (including denitrification) would have added further realism (Fennel et al., 2011). However, our main focus was to evaluate zooplankton dynamics in the oligotrophic region where higher trophic levels that depend on mesozooplankton secondary production may experience food limitation and where benthic processes are negligible.

NEMURO was chosen for the present study because it distinguishes SZ, LZ and PZ, permitting a detailed analysis of dynamics for multiple functional types in the GoM zooplankton community. During initial GoM simulations, default NEMURO parameterizations, configured for the North Pacific (Kishi et al., 2007), substantially overestimated surface nitrate, surface Chl and mesozooplankton biomass relative to observations. We attribute these differences to (1) substantially higher temperatures in the GoM compared with the North Pacific, which significantly increase decomposition and growth rates in the model, resulting in higher nutrient recycling and elevated near-surface stocks of phytoplankton and zooplankton; and (2) distinct differences in taxonomic composition of phytoplankton and zooplankton communities in the GoM and North Pacific, with significant differences in key parameter values associated with growth and grazing.

For more details on the specific processes represented and the interactions between state variables in NEMURO, we direct readers to Kishi et al. (2007). All model equations are provided in the Supplement to this paper. Biogeochemical model forcing, initial and open boundary conditions are also outlined in Supplement Sect. S1. Briefly, daily-average shortwave radiation fields obtained from Climate Forecast System Reanalysis (CFSR) were used to force light limitation of phytoplankton. Once a final parameter set was determined (see Sect. 2.1.3), initial and open boundary conditions for all state variables were prescribed from a spun up idealized one-dimensional version of NEMURO-GoM. After initializing, the three-dimensional model was spun up over four years before conducting the full 20-year experiment. River nutrient input from the Mississippi River was prescribed using nitrate samples collected by United States Geological Survey (USGS) and due to a lack of observations for other rivers was prescribed for all 37 rivers represented in the model.

2.1.2 Modifications to default NEMURO model

To improve realism for application to the GoM, five structural changes were made to the original NEMURO model. First, we removed the SP to LZ grazing pathway. The original SP state variable for the North Pacific represents nanophytoplankton (e.g., coccolithophores), which can be important prey of copepods and other mesozooplankton. In the GoM, however, cyanobacteria and picoeukaryotes (too small for direct feeding by most mesozooplankton) comprise much of the phytoplankton biomass and hence are represented as SP in our model. In addition to adding ecological realism, this change in direct trophic connection between SP and LZ allowed the model to produce a more realistic LP-dominated phytoplankton community on the shelf (see Discussion).

Next, quadratic mortality was replaced with linear mortality for all biological state variables with the exception of predatory zooplankton (PZ). In biogeochemical models, quadratic mortality is often used for numerical stability and/or to represent implicit loss terms to an unmodeled parasite or predator that covaries in abundance with its prey (e.g., viral lysis of phytoplankton or predation by unmodeled higher predators) (Anderson et al., 2015). However, grazing mortality is explicitly modeled in NEMURO, and viral mortality is generally not a substantial loss term for bulk phytoplankton (Staniewski and Short, 2018). Quadratic mortality was retained for PZ to account for predation pressure of unmodeled predators (e.g., planktivorous fish). During the model tuning process, we found that removal of quadratic mortality from the four other plankton functional groups was an important parameterization change that allowed the model to simulate more realistic mesozooplankton biomass in the oligotrophic GoM (see Discussion).

The default ammonium inhibition term and light limitation functional form in NEMURO were replaced in NEMURO-GoM with more widely adopted parameterizations. The exponential ammonium inhibition term in the nitrate limitation function was replaced with the term described by Parker (1993), as has been done in previous PBM studies (Fennel et al., 2006) due to the nonmonotonic behavior of the default NEMURO ammonium inhibition term. At high NO₃ concentrations, the default term is known to generate unrealistic phytoplankton nutrient uptake patterns in which total nutrient uptake (i.e., uptake of NO₃ plus uptake of NH₄) can actually decrease despite increases in NH₄ (and constant NO₃).

Light limitation in NEMURO is based on an optimal light parameterization that implicitly includes photoinhibition. This formulation was replaced with the Platt et al. (1980) functional form that allows one to explicitly control the amount of photoinhibition, which can be important in the GoM where surface irradiances are high. Additionally, the Platt functional form is commonly used, and thus parameter values are easier to find for comparison (e.g., initial slope of the PI curve, α). This formulation is also implemented in newer versions of NEMURO, such as the code used in the Regional Ocean Modeling System (ROMS) NEMURO biogeochemical package.

Finally, to account for photoacclimation and more accurately simulate deep chlorophyll maximum (DCM) dynamics, we replaced the constant C:Chl parameter with a variable C:Chl model where ratios for SP and LP were allowed to vary based on the formulation described by Li et al. (2010), which considers both light and nutrient limitation (see Supplement). The Li et al. (2010) equations build on a previously constructed dynamic regulatory model of phytoplankton physiology which describes C:Chl variability under balanced growth and nutrient-saturated conditions at constant temperature (see Geider et al., 1998). Herein, “default” NEMURO includes the modified ammonium inhibition, light formulation and variable C:Chl model.

2.1.3 NEMURO-GoM model tuning procedure

In total, NEMURO includes 71 parameters, 23 of which were modified in the present study. For initial model tuning, we used an idealized one-dimensional model designed to mimic the oligotrophic GoM. To guide our tuning procedure, we relied on a semiquantitative approach where the one-dimensional model solution was evaluated based on five ecosystem benchmarks. Target values for benchmarks and other ecosystem attributes were determined from observations or a theoretical basis. Ecosystem benchmarks included surface Chl, mesozooplankton biomass, DCM depth, DCM magnitude, and SP:LP ratio. Surface Chl and mesozooplankton biomass were chosen as benchmarks to evaluate the realism of plankton biomass in the model. The DCM depth and magnitude were chosen to evaluate the vertical structure of the simulated ecosystem, and the SP:LP ratio was used to gauge the realism of the plankton community composition (i.e., high SP:LP is expected in the oligotrophic GoM). The model was also tuned by considering the relative magnitudes of loss terms for phytoplankton (grazing, mortality, respiration and excretion), total protistan zooplankton grazing relative to mesozooplankton grazing, and surface and deep nitrate concentrations. We outline each parameter change, justification and the resulting impact on the ecosystem benchmarks simulated by the idealized one-dimensional model in Sect. S3. Where possible, we modified parameters in groups so that relative changes were consistent throughout the model (e.g., doubling all zooplankton mortality terms). After tuning in the one-dimensional model, parameter sets were implemented into the full three-dimensional model where additional tuning was performed. Once a final parameter set was determined we conducted a parameter sensitivity analysis over 18 individual experiments to identify impacts of parameter changes from default NEMURO values (Sect. S5).

2.2.1 Description of the offline numerical environment

To run large numbers of three-dimensional simulations efficiently for basin-scale tuning, NEMURO-GoM was run offline using the MITgcm offline tracer advection package. MITgcm was selected as it contains convenient packages for running offline simulations (McKinley et al., 2004). That is, the dynamical equations of motion are not computed during the NEMURO-GoM integration, but rather the physical prognostic variables (i.e., temperature, salinity and three-dimensional velocity fields) are prescribed from daily-averaged flow fields saved from a previous hydrodynamic model integration. This allows the recycled use of flow fields, leaving only the tracer equations to be computed. In the offline MITgcm package, the prognostic variables provide input to an advection scheme and mixing routine that conservatively handles offline advection and diffusion of the biogeochemical tracer fields. MITgcm has many options for linear and nonlinear advection schemes. Here we use a third-order direct space–time flux-limiting scheme. Sub-grid-scale mixing of the biogeochemical fields is handled offline through the nonlocal K-profile parameterization (KPP) package based on mixing schemes developed by Large et al. (1994). For more information about the MITgcm packages, we direct readers to the MITgcm manual (http://mitgcm.org/, last access: 29 June 2020).

There are two main advantages to running PBMs in an offline environment: (1) the momentum equations are not integrated during the model run; and (2) the physical time step is no longer bound by the dynamical Courant–Friedrichs–Lewy (CFL) numerical stability criterion, which together significantly reduces the computational cost. Instead, the stability of the tracer advection scheme and timescales needed to resolve biological/physical processes of interest set the limits on the time steps and prescription frequencies of flow fields. When the physical time step is shorter than the flow field prescription frequency, a simple linear interpolation of the flow fields is performed between time steps. Offline simulations of tracer advection have been found to closely resemble online runs (that is, computed together with the integration of the hydrodynamic model's prognostic equations) when the three-dimensional flow fields are prescribed at a frequency that is at or below the inertial period ($T=\mathrm{2}\mathit{\pi }/f$, T_GoM > 24 h) for a region (Hill et al., 2005).

In the present study, the offline time step (30 min) is an order of magnitude greater than the hydrodynamic model's (HYCOM-GoM, described in Sect. 2.2.2) baroclinic time step (120 s). For reference, HYCOM-GoM required ∼76 d to run to completion on 64 parallel cores. These time requirements would increase considerably with the 11 additional biogeochemical tracers used in NEMURO. In contrast, NEMURO-GoM ran significantly faster, taking a total of ∼50 h on 80 parallel cores. Offline models offer a valuable tool for integrating PBMs particularly as spatial resolution and complexity in these models continue to increase (e.g., DARWIN, Follows et al., 2007; GENOME, Coles et al., 2017). While computationally advantageous, however, offline simulations have inherently greater input and output (I/O) demands that can become bottlenecks in some applications. Issues with conservation can also arise as three-dimensional advection schemes are only approximately positive definite.

2.2.2 Description of the offline dynamical fields

The NEMURO-GoM model is forced by daily-averaged three-dimensional velocity, temperature and salinity fields from a preexisting 20-year (1993–2012) HYCOM (Chassignet et al., 2003) regional GoM hindcast (H-GoM). H-GoM is based on version 2.2.99B of the HYCOM code, originally provided by the Naval Oceanographic Office (NAVOCEANO) Major Shared Resource Center. H-GoM was run at 1/25^∘ (∼4 km) horizontal resolution with 36 vertical hybrid coordinate layers and assimilated historic, in situ and satellite observations. The domain encompasses the entire GoM and extends south of the Mexico—Cuba Yucatán Channel to 18^∘ N and as far east as 77^∘ W (Fig. 1). Further details on H-GoM (experiment ID: GOMu0.04/expt_50.1), including model forcing and the main model configuration file (i.e., blkdat.input_501), can be found at https://www.hycom.org (last access: 29 June 2020).

The H-GoM flow fields were mapped from the HYCOM native hybrid vertical coordinate to z levels used by the MITgcm. NEMURO-GoM was configured for 29 vertical z levels (10 m intervals from 0 to 150 m; 25 m intervals from 150 to 300 m; 50 m intervals from 300 to 500; and 1000, 2000 and ∼4000 m). Mapping was performed by computing total zonal and meridional transports across the lateral boundaries of each MITgcm grid cell (e.g., 0–10 m bin; which may include multiple HYCOM layers) and then dividing by the area of the respective cell face. This vertical mapping approach is consistent as both HYCOM and MITgcm use an Arakawa C-grid orientation for model variables. The H-GoM bathymetry was adjusted such that no partial cells existed in the domain to avoid thin cells. The continuity equation was subsequently used to calculate vertical velocities. The use of transports in this approach ensures conservation and approximately identical profiles of vertical velocity to those in H-GoM fields. For mapping of temperature and salinity fields (used in the KPP mixing routine and for scaling biological temperature-dependent rates), a simple linear interpolation was performed.

2.3.1 Surface chlorophyll observations

A benchmark for surface Chl was determined using the Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) product from the Ocean Biology Processing Group (OBPG) of the National Aeronautics and Space Administration (NASA). The product used here is the mapped, level-3, daily, 9 km resolution images from 4 September 1997 to 10 December 2010 processed according to the algorithm of Hu et al. (2012). To compute model–data point-to-point comparisons, we take the corresponding daily-averaged simulated surface Chl field and interpolate to the SeaWiFS grid before applying the daily cloud coverage mask corresponding to the matching SeaWiFS image. In total 4291 daily images consisting of 22 244 513 nonzero cell values (herein referred to SeaWiFS measurements) were used to validate NEMURO-GoM. Approximately 500–1200 daily model–data point-to-point comparisons were made for each SeaWiFS grid cell (Fig. 1).

Figure 1(a–e) Spatial and temporal coverage of observational datasets used for model validation. Total number of nonzero values available in daily SeaWiFS images from September 1997 to 10 December 2010 along with nitrate profiles from the World Ocean Database (black dots) and the location of in situ measurements collected during May 2017 (circles) and May 2018 (triangles) process study cruises (a). Total annual sampling of the SEAMAP surveys from 1983 to 2017 with samples overlapping with the PBM simulation period denoted in red (b). SEAMAP samples organized by month (d). Total sample density within each $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ box (c) along with the number of years with at least one sample (e). The 1000 m isobaths and coastline are denoted by black continuous lines.

2.3.2 Mesozooplankton biomass observations

To evaluate model mesozooplankton biomass estimates, we used plankton tows collected during Southeast Area Monitoring and Assessment Program (SEAMAP) surveys in the northern and central GoM. In total, 11 781 plankton tows were collected from 1983 to 2017, with two main annual surveys in the spring (offshore) and fall (shelf) (Fig. 1). On average, SEAMAP collected approximately 300 samples per year with a specific sampling array offshore and more general sampling coverage on the shelf. In total, 6835 samples were used for direct point-to-point model–data comparisons. Samples were collected using standard gear consisting of a 61 cm diameter bongo frame fitted with two 333 µm mesh nets. The nets were fished in a double-oblique tow pattern from the surface down to 200 or 5 m off the bottom and back to the surface. Simultaneous samples were also collected using a 202 µm mesh net during 82 tows. Of these samples, roughly half were collected in the oligotrophic GoM. The average ratio between biomass measured in the 333 and 202 µm bongo tows (0.5093±0.12) was used to convert 333 µm samples so that direct comparisons could be made with model mesozooplankton (LZ+PZ) biomass fields. Zooplankton biomass was originally quantified as displacement volumes (DVs). Carbon mass (CM) equivalents were subsequently calculated as log₁₀(CM) = (log₁₀(DV) +1.434)/0.820 (Wiebe, 1988; Moriarty and O'Brien, 2013).

For comparison to in situ biomass measurements, size ranges were assumed for zooplankton state variables (see Sect. 2.3.4). In NEMURO-GoM, SZ represent heterotrophic protists (e.g., ciliates) and early stages of mesozooplankton, which are typically < 200 µm in the ocean. Hence, SZ are assumed to range in size from 2 to 200 µm. LZ are considered to represent small mesozooplankton ranging in size from 0.2 to 1.0 mm, while PZ represent large mesozooplankton ranging in size from 1.0 to 5.0 mm. For comparison to the SEAMAP climatology the model LZ and PZ fields were first converted to units of carbon assuming Redfield C:N ratio and summed over the entire simulation before being depth averaged to the bottom or 200 m. For point-to-point model–data comparisons, total mesozooplankton biomass fields were interpolated to SEAMAP sample locations/times and then depth averaged to the corresponding sample tow depth.

2.3.3 Observed vertical profiles of chlorophyll and nitrate

Depth profiles of Chl were also collected during SEAMAP surveys using a Sea-Bird WETStar fluorometer attached to a CTD. Calibration of the fluorimeter was infrequent, and thus profiles were used to determine the depth of the fluorescence maxima for comparisons to DCM depths in the model. In total, 2435 profiles were collected from 2003 to 2012, with 1052 profiles overlying bottom depths > 1000 m. Profiles were available for earlier SEAMAP surveys; however, no standard QA/QC protocol for fluorometer data was in place prior to 2003.

To evaluate DCM magnitudes in the model, we used 145 fluorescence profiles collected during May 2017 and 2018 process study cruises (see Sect. 2.3.4). The fluorometer was attached to a CTD and calibrated using 126 in situ Chl samples. Chl concentrations were determined from filtered samples collected at depths ranging from 5 to 115 m using high-performance liquid chromatography (HPLC). Since the cruise sampling does not overlap with our NEMURO-GoM simulation period, model–data comparisons were made for all 20 years of the model run using sample locations and time of the year. This was also done with other field measurements from the process cruises (see Sect. 2.3.4).

For model–data comparisons of nitrate, we utilized profiles from the World Ocean Database (WOD). In total, 96 profiles were available during our simulation period and located in the oligotrophic GoM (> 1000 m isobath). Profiles were collected during all months except March and December, with the majority of samples collected during May, July and August (Fig. 1a).

2.3.4 Biomass and rate measurements from process study cruises

Although in situ rate measurements are made much less frequently than biological standing stock measurements, they offer very powerful constraints for validating the internal dynamics of a biogeochemical model (Franks, 2009). Consequently, we made phytoplankton and zooplankton rate measurements on two cruises in the open-ocean GoM in May 2017 and 2018 and used these measurements to validate the model (Fig. 1a). On the process study cruises, we utilized a quasi-Lagrangian sampling scheme to investigate plankton dynamics in the oligotrophic GoM. Two drifting arrays (one sediment trap array and one in situ incubation array) were deployed to serve as a moving frame of reference during ∼4 d studies (cycles) characterizing the water parcel (Landry et al., 2009; Stukel et al., 2015). During these cycles, we measured daily profiles of Chl, photosynthetically active radiation, phytoplankton growth rates and productivity, protistan grazing rates, and size-fractionated mesozooplankton biomass and grazing rates.

Size-fractionated mesozooplankton biomass and grazing rates were determined from daily day–night paired oblique ring-net tows (1 m diameter, 202 µm mesh). In total, 40 oblique bongo net tows (16 in 2017 and 24 in 2018) sampled the oligotrophic GoM mesozooplankton community from near surface to a depth ranging from 100 to 135 m. Upon recovery, the sample was anesthetized using carbonated water, split using a Folsom splitter, filtered through a series of nested sieves (5, 2, 1, 0.5 and 0.2 mm), filtered onto preweighed 200 µm Nitex filters, rinsed with isotonic ammonium formate to remove sea salt and flash frozen in liquid nitrogen. In the lab, defrosted samples were weighed for total wet weight and subsampled in duplicate (wet weight removed) for gut fluorescence analyses. The remaining wet sample was dried and subsequently reweighed and combusted for CHN analyses to determine total dry weight and C and N biomasses. Gut fluorescence subsamples were homogenized using a sonicating tip, extracted in acetone, and measured for Chl and phaeopigments using the acidification method. The phaeopigment concentrations in the zooplankton guts were the basis for calculated grazing rates using gut turnover times based on temperature relationships for mixed zooplankton assemblages. For additional details, see Décima et al. (2011, 2016).

Protistan grazing rates were measured using the two-point, mini-dilution variant of the microzooplankton grazing dilution method (Landry et al., 1984, 2008; Landry and Hassett, 1982). Briefly, one 2.8 L polycarbonate bottle was gently filled with whole seawater taken from six depths (from the surface to the depth of the mixed layer). A second 2.8 L bottle was then filled with 33 % whole seawater and 67 % 0.2 µm filtered seawater. Both bottles were then placed in mesh bags and incubated in situ at natural depths for 24 h. These experiments were conducted on each day of the ∼4 d cycle. After 24 h, the bottles were retrieved and filtered onto glass fiber filters, and Chl concentrations were determined using the acidification method (Strickland and Parsons, 1972). Net growth rates (k=ln(Chl_final/Chl_init)) in each bottle were determined relative to initial Chl samples. Phytoplankton specific mortality rates resulting from the grazing pressure of protists were calculated as m= (k_d–k₀)/(1–0.33), where k_d is the growth rate in the dilute bottle and k₀ is the growth rate in the control bottle. Phytoplankton specific growth rates were calculated as $\mathit{\mu }=k_{\mathrm{0}}+m$. For additional details, see Landry et al. (2016) and Selph et al. (2016). Phytoplankton net primary production was quantified at the same depths by ${\mathrm{H}}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{3}}^{-}$ uptake experiments. Triplicate 2.8 L polycarbonate bottles and a fourth “dark” bottle were spiked with ${\mathrm{H}}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{3}}^{-}$ and incubated in situ for 24 h at the same sampling depths as for the dilution experiments. Samples were then filtered and the ¹³C:¹²C ratios of particulate matter determined by isotope ratio mass spectrometry.

4.2.1 Surface chlorophyll discrepancies

Within our model–data comparisons of surface Chl we find that NEMURO-GoM reproduces important patterns in both the oligotrophic and shelf region, the latter of which, apart from the northern shelf, has not been well resolved by previous PBMs (e.g., Gomez et al., 2018; Xue et al., 2013). The absence of a shelf Chl signature may, in some cases, be overly attributed to bias in satellite measurement due to high concentrations of colored dissolved organic matter (CDOM). While a clear shelf signature is well resolved in NEMURO-GoM, the model–data mismatch is greater on the shelf compared to oligotrophic regions. This is an expected result considering that the model incorporates climatological river forcing while actual variability is much more complex. Furthermore, the absence of CDOM in the model likely contributes to the overestimation of phytoplankton biomass on the shelf.

In future studies, the inclusion of daily nutrient data like that produced for the Mississippi River by USGS starting in 2011 is needed for PBMs to better resolve variability on the shelf. Including benthic processes, such as denitrification (Fennel et al., 2006), may also reduce model–data mismatch in shelf regions. Implementing more realistic light attenuation (e.g., wavelength-specific light attenuation or inclusion of CDOM) could further improve estimates of phytoplankton biomass on the shelf as primary production can be sensitive to different light attenuation formulations (Anderson et al., 2015). In our model, it was difficult to simulate deep DCMs in the oligotrophic GoM while also simulating DCMs on the shelf that were shallow enough to maintain high nitrate. This may reflect the need for more realistic light attenuation in the model.

Quantifying uncertainty in C:Chl ratios is also an important task moving forward as future PBMs will likely continue to depend heavily on satellite Chl for the bulk of model validation. To correctly parameterize C:Chl ratios, more in situ samples are needed that resolve changes in phytoplankton light-harvesting pigments along gradients from coastal to oligotrophic regions and from the surface to the DCM. Without these observations, it is difficult to gauge mismatches between model and satellite ocean color products (e.g., discrepancies on the Campeche Bank in our model) or in situ profiles of Chl. In addition to validation, these measurements are needed to avoid erroneous model tuning. For instance, a model that exhibits significant mismatch with respect to surface Chl may in fact accurately estimate carbon-based phytoplankton biomass while using unrealistic C:Chl ratios. One could arrive at incorrect conclusions about regional ecosystem dynamics as a result of modifying model parameters or structure in an effort to better fit Chl observations. Given the importance of C:Chl ratios in PBMs, future studies should quantify uncertainty in modeled Chl through sensitivity experiments focused on C:Chl model parameters and formulation, with explicit comparison to direct field measurements of phytoplankton C:Chl.

In our model, the most noticeable surface Chl model–data mismatch occurs on the southern GoM shelf (Campeche Bank), where the model consistently overestimates surface Chl. This bias was also notable in the GoM PBM implemented by Damien et al. (2018), particularly in winter. We believe this discrepancy is driven by a combination of errors involving overestimation of shelf mixing by the hydrodynamic model, entrainment of high-Chl water (given the overestimated DCM magnitude in the model) or errors in the open boundary conditions which result in an overestimation of upwelled nutrients/biomass near the Yucatán Peninsula that are transported westward by shelf currents. We found that the Campeche Bank model–data mismatch was reduced when open boundary conditions included nitracline depths of greater than 100 m. This may reflect realistic in situ conditions considering that Caribbean water entering the GoM is highly oligotrophic. During our process cruises, nitrate was often undetectable above 100 m in samples collected near the Loop Current (Angela Knapp, personal communication, 2018).

Although modifying the boundary conditions may be justified, deepening the nitracline at the boundaries made it increasingly difficult to sustain realistic surface phytoplankton biomass in the oligotrophic GoM. This may point to the importance of nitrogen-fixing cyanobacteria, which provide an alternative source of new nitrogen (other than upwelling and mixing) that could be supporting phytoplankton at the surface given the strong stratification and deep nitraclines in the GoM. In the process of model tuning, we noticed that increasing the DON pool by increasing the PON to DON decomposition rate was necessary to maintain both relatively deep nitraclines and realistic surface Chl by providing a slow leeching of ammonium near the surface through bacterial communities. The need for this slow production of ammonium in surface layers may compensate for nitrogen fixation, which is not included in NEMURO (Holl et al., 2007; Mulholland et al., 2006). In future studies, including diazotrophs as a separate phytoplankton functional type would be essential for evaluating the importance of nitrogen fixation in the GoM.

Despite the model–data mismatch on the Campeche Bank, this discrepancy appears to have little impact on the rest of the GoM. However, the model overestimates surface Chl in the southwestern GoM, which can likely be attributed to entrainment of high-Chl water originating from the Campeche Bank. Locally, the ecological impact is likely more significant. Higher phytoplankton biomass would be expected to support higher mesozooplankton grazing rates and secondary production. Indeed, some of the highest model rates of secondary production occur on the Campeche Bank. Hence, the surface Chl model–data mismatch may lead to an overestimation of secondary production for this region.

4.2.2 Deep chlorophyll maximum discrepancies

Since most PBMs focus on validating against satellite-derived surface chlorophyll, the dynamics of the DCM is often insufficiently investigated. Consequently, many models predict DCM depths that are far too shallow. Identifying this issue in the literature proved to be difficult because most studies do not provide profiles of simulated Chl (an exception is the recent GoM PBM by Damien et al., 2018). We note that DCM depths in the DIAZO model (Stukel et al., 2014) were often quite shallow or completely nonexistent in the portion of the domain that included the oligotrophic GoM region. Underestimates of DCM depth in the unmodified COBALT biogeochemical model have also been identified (Moeller et al., 2019). In our investigation of the GoM PBM implemented by Gomez et al. (2018), we found that DCMs in the oligotrophic region were commonly shallow and weak. In the default NEMURO simulation, DCM depths in the oligotrophic region were typically at a depth of 25 m, which is much shallower than observed (SEAMAP: 80±25 m, Process cruises: 107±21 m). While this issue may seem insignificant, particularly if a study is focused on mixed-layer dynamics, accurate placement of the DCM can have profound impacts on PBM behaviors, because the DCM is typically colocated with the nitracline. Unrealistically shallow DCMs and nitraclines permit unrealistically high nitrate fluxes into the surface layer following mixing events; thus, validating the DCM in PBMs is critical.

For these reasons, we devoted substantial effort to tuning phytoplankton dynamics in the DCM. Modifications to α (the slope of the photosynthesis–irradiance curve) and attenuation coefficients allowed the model to estimate deeper more realistic DCM depths. Inclusion of a variable C:Chl module was also implemented to better resolve the DCM. However, an additional issue was present in the default NEMURO simulations, the NEMURO-GoM, and every simulation that we attempted. In all simulations that formed DCMs, the location of the DCM was always colocated with a maximum in phytoplankton specific growth rate, even though field measurements indicate that phytoplankton growth rates and NPP are either relatively constant with depth or decline in the DCM. This is not surprising, given the low photon flux at the base of the euphotic zone and the energetic demands required to upregulate cellular density of light-harvesting pigments. Additionally, our field measurements show that the DCM was not associated with a biomass maximum (biomass was fairly constant with depth), suggesting that DCM formation in the GoM is physiologically driven.

We believe this DCM dynamical issue was responsible, in part, for the underestimation of specific phytoplankton growth and microzooplankton grazing rates by the model despite estimating higher NPP (Fig. 4d). An underestimation of C:Chl ratios or overestimation of phytoplankton biomass may also contribute to the model–data mismatch at the DCM. The ecological impact of this mismatch is likely to be greatest in the model during the summer months when the mixed-layer depth is shallow and the DCM reaches its maximum depth. Elevated growth rates or biomass at the DCM during this period could significantly decrease nitrate flux into the mixed layer, causing the model to underestimate near-surface primary and secondary production.

Future PBM studies need to focus more effort on resolving ecological dynamics responsible for the formation of the DCM. Errors originating from the hydrodynamic model or use of temporally averaged velocity fields used in offline models may also contribute to model–data mismatch at the DCM. Vertical mixing is particularly important to PBMs but is often poorly validated in hydrodynamic models. Greater coordination between physical modelers and biologists to constrain vertical fluxes should be considered an important avenue for improving PBM simulations moving forward.

4.2.3 Mesozooplankton grazing discrepancies

Novel to this study, model estimates of mesozooplankton biomass are shown to agree closely with observations on the shelf and in the oligotrophic GoM. To our knowledge, this study includes the first quasi-regional zooplankton biomass model validation in a PBM. Our model also provides the first model–data comparisons of size-specific zooplankton biomass and grazing rates for the GoM. Such comparisons provide valuable insights into the potential biases of traditional functional-group biogeochemical models pertaining to zooplankton dynamics (Everett et al., 2017). While NEMURO-GoM shows broad agreement with zooplankton observations, some model–data mismatch occurs, particularly for mesozooplankton grazing rates.

We identify three factors that may explain the model–data mismatch for mesozooplankton grazing rates. The first and most obvious factor is the temporal sampling discrepancy as measurements were collected outside our model simulation period. Model–data mismatch may also arise from inaccuracies in the field measurements. During our process cruises, the zooplankton gut pigment measurements were based solely on phaeopigment content due to phytodetrital aggregates and Trichodesmium colonies found in our zooplankton net tows, which can lead to substantial contamination. Thus, true mesozooplankton grazing rates were likely underestimated because undegraded Chl can be abundant in the foreguts of mesozooplankton. Furthermore, the gut pigment approach assumes that any group of mesozooplankton has a constant gut throughput time (as a function of temperature), which is an oversimplification.

Uncertainties in model grazing formulations could also contribute to model–data mismatch (Gentleman et al., 2003a; Sailley et al., 2015). Future in situ grazing measurements are needed to enable an objective selection of grazing formulations and parameter values. In particular, field studies that shed light on prey selectivity would be useful for parameterizing PBMs with multiple mesozooplankton functional groups, such as NEMURO-GoM. Such studies are challenging, however, because the difficulty of making in situ grazing measurements on mesozooplankton, combined with the inherent uncertainty of these measurements, can make it challenging to differentiate between, for instance, Ivlev and Holling's disk grazing formulations (e.g., Fig. 4 in Morrow et al., 2018). Nevertheless, differences in parameterizing grazing can lead to substantially different model behavior (Anderson et al., 2010; Sailley et al., 2015; Wainwright et al., 2007). In NEMURO-GoM, secondary production and dietary preferences of the mesozooplankton community are both strongly influenced by model grazing formulation. While we carefully chose parameterizations that gave reasonable fits to extensive field datasets of zooplankton biomass and grazing rates, this does not preclude the possibility that other functional forms would have more accurately simulated zooplankton dynamics. Hence, future PBMs should investigate how different grazing formulations impact zooplankton dynamics in the region. We especially recommend collaborations between experimentalists (potentially using new techniques such as DNA metabarcoding of gut contents) and modelers to develop synthetic approaches with the potential to quantitatively assess the realism of different grazing formulations.

Clear model–data mismatch is also evident in the proportion of grazing mediated by PZ and LZ. This may be due to the fact that PZ is by default explicitly defined and parameterized as a higher-trophic-level mesozooplankton that can feed on LZ. In reality, while there is a correlation between size and trophic level in the ocean, many predatory zooplankton are < 1 mm, and many suspension-feeding zooplankton are > 1 mm; hence, the overlap of taxonomic groups with different functional roles and sizes makes it difficult to directly compare model categories to field data. For example, shelf suspension-feeding zooplankton are likely larger than their counterparts in the oligotrophic GoM although their functional role in the ecosystem does not change between environments.

The ecological impact of the model's potential overestimation of LZ grazing rates is most likely to manifest through an increase in the ratio of secondary production to mesozooplankton biomass. Since both LZ and PZ biomasses are accurately modeled by NEMURO-GoM, the overestimation of grazing rates suggests that LZ turnover times may be too high, thus leading to higher estimates of secondary production. However, this interpretation may oversimplify the complex interactions within pelagic protistan communities. In the oligotrophic region where our model overestimates LZ grazing rates, the model indicates that heterotrophic protists comprise approximately half of the LZ diet. Thus, overestimates of grazing on LP do not necessarily lead to overestimates in total consumption if < 1 mm zooplankton derive substantial nutrition from nonphototrophic sources in the field. Furthermore, the model's construction (i.e., LZ and PZ are functional groups, while the field data are size classes) suggests that part of the model–data mismatch in Fig. 5d may result from the presence of some suspension feeders (i.e., LZ) in the > 1 mm zooplankton size class and some carnivorous zooplankton (i.e., PZ) in the < 1 mm zooplankton size class. In this case, the model may simply attribute too high of a LP:SZ prey ratio to LZ. If this is the issue, the model's estimate of LZ secondary production may be accurate but its trophic level too low (or, conversely, the trophic level of PZ too high). Direct assessments of zooplankton trophic position (e.g., by compound-specific isotopic analysis of amino acids, Chikaraishi et al., 2009; Décima et al., 2017) may help resolve these issues.