Seasonal patterns in phytoplankton biomass across the northern and deep Gulf of Mexico: a numerical model study

Biogeochemical models that simulate realistic lower-trophic-level dynamics, including the representation of main phytoplankton and zooplankton functional groups, are valuable tools for improving our understanding of natural and anthropogenic disturbances in marine ecosystems. Previous three-dimensional biogeochemical modeling studies in the northern and deep Gulf of Mexico (GoM) have used only one phytoplankton and one zooplankton type. To advance our modeling capability of the GoM ecosystem and to investigate the dominant spatial and seasonal patterns of phytoplankton biomass, we configured a 13-component biogeochemical model that explicitly represents nanophytoplankton, diatoms, micro-, and mesozooplankton. Our model outputs compare reasonably well with observed patterns in chlorophyll, primary production, and nutrients over the Louisiana–Texas shelf and deep GoM region. Our model suggests silica limitation of diatom growth in the deep GoM during winter and near the Mississippi delta during spring. Model nanophytoplankton growth is weakly nutrient limited in the Mississippi delta year-round and strongly nutrient limited in the deep GoM during summer. Our examination of primary production and net phytoplankton growth from the model indicates that the biomass losses, mainly due to zooplankton grazing, play an important role in modulating the simulated seasonal biomass patterns of nanophytoplankton and diatoms. Our analysis further shows that the dominant physical process influencing the local rate of change of model phytoplankton is horizontal advection in the northern shelf and vertical mixing in the deep GoM. This study highlights the need for an integrated analysis of biologically and physically driven biomass fluxes to better understand phytoplankton biomass phenologies in the GoM.

1 Introduction

The Gulf of Mexico (GoM) is characterized by large spatial differences in plankton productivity and biomass, ranging from the oligotrophic Loop Current to the highly productive northern shelf. Productivity in this last region is strongly influenced by river runoff. The Mississippi–Atchafalaya (MS-A) river system is the largest river input with a mean river discharge of 21 524 m³ s⁻¹ (Aulenbach et al., 2007), contributing more than 80 % of the entire dissolved inorganic nitrogen (DIN) load into the northern GoM (Xue et al., 2013). The large plankton production and vertical stratification driven by the MS-A river system discharge promote the development of a hypoxic bottom layer a few meters thick off Louisiana and Texas during summer (Obenour et al., 2013). This hypoxic layer can negatively impact metabolism and growth of fish and invertebrates (Rosas et al., 1998; Craig and Crowder, 2005), and disturb species distribution and composition (Craig, 2012). The influence of river runoff on plankton production substantially decreases offshore (Green and Gould, 2008). In the oligotrophic deep GoM, the spatiotemporal patterns in phytoplankton biomass are mainly associated with seasonal changes in thermal stratification and mesoscale ocean dynamics (e.g., Muller-Karger et al., 2015).

Multiple ocean–biogeochemical modeling studies have been conducted in the northern GoM to understand the drivers of phytoplankton biomass variability, carbon export, nutrient cycling, and bottom hypoxia variability. Green et al. (2008) configured a zero-dimensional Lagrangian model of the Mississippi (MS) river plume, which included two types of phytoplankton (small and large sized), two types of zooplankton (micro- and mesozooplankton), bacteria, detritus, ammonium, and nitrate. This study derived distinct production patterns for small- and large-sized phytoplankton production, concluding that primary production was mainly limited by physical dilution of nitrate, light attenuation, and the sinking of diatoms (large phytoplankton). More complex modeling efforts for the region include a series of three-dimensional (3-D), fully coupled ocean–biogeochemical models, based on Fennel's biogeochemical model (Fennel et al., 2006). The original Fennel's model formulation included ammonium, nitrate, phytoplankton, chlorophyll, zooplankton (representing mesozooplankton), and two detritus types as state variables. Fennel et al. (2011) examined the underlying factors determining seasonal patterns in phytoplankton biomass in the Louisiana–Texas shelf and concluded that phytoplankton production was not nitrogen limited near the MS delta. They also showed that zooplankton grazing played an important role in defining phytoplankton biomass changes and speculated that physical transport of phytoplankton could impact biomass seasonality. Xue et al. (2013) configured Fennel's model for the entire GoM, describing main spatiotemporal patterns in plankton biomass and DIN in the coastal and oceanic domains. However, since they did not investigate underlying drivers (production, biomass losses) of phytoplankton biomass as was done in Fennel et al. (2011), less is known about the factors modulating the seasonality of phytoplankton in the deep GoM.

Significant differences in plankton production and carbon export can be expected between food webs dominated by small-sized (nanophytoplankton, microzooplankton) and large-sized (diatoms, mesozooplankton) plankton components. Sedimentation rates are enhanced (decreased) in diatom (small phytoplankton)-based food webs, and therefore changes in phytoplankton composition could influence bottom remineralization processes (Dortch and Whiteledge, 1992; Dagg et al., 2003; Green et al., 2008; Zhao and Quigg, 2014). In addition, changes in phytoplankton composition may modulate trophodynamics, which can impact the reproductive success of upper trophic levels and therefore modulate marine population abundance (Rykaczewski and Checkley, 2008). In the GoM, 3-D regional ocean–biogeochemical models that include more than one plankton functional group have been implemented only for the western Florida shelf (Walsh et al., 2003). New modeling efforts are required to examine spatiotemporal patterns of main phytoplankton functional groups across the northern and deep GoM. A key modeling aspect is the characterization of diatoms and nanophytoplankton growth. It is well known that (1) nanophytoplankton uptake nutrients more efficiently than diatoms; (2) diatoms can achieve greater growth rates than nanophytoplankton in nutrient-rich environments; and (3) diatoms require silicate as an additional nutrient for frustule formation (Litchman et al., 2006; Falkowski and Oliver, 2007). These differences should be considered when simulating phytoplankton responses to changes in nutrient availability.

The present study explores underlying factors determining spatial and seasonal patterns in phytoplankton biomass across the coastal and ocean domains in the GoM, using an ocean–biogeochemical model that explicitly simulates small- and large-sized plankton groups. After validating the model results with available observations, we examine the main seasonal patterns of phytoplankton biomass. Our main goals are (1) to describe the spatiotemporal patterns in growth limitation for diatoms and nanophytoplankton and (2) to evaluate the coupled role of biological (phytoplankton production and biological losses) and physical (advection and turbulent diffusion of biomass) processes as drivers of phytoplankton seasonality. This study complements Fennel et al. (2011) on phytoplankton variability in the northern GoM, by adding complexity to the modeled lower-trophic-level dynamics, extending the description of phytoplankton growth limitation patterns to the deep GoM, and quantifying the role of advection and diffusion.

2 Data and model

2.1 Data

Monthly mean composite fields of Sea-Viewing Wide Field-of-View Sensor (SeaWiFS, 1998–2011) and Moderate Resolution Imaging Spectroradiometer (MODIS, 2003–2014) chlorophyll a were retrieved from the Institute for Marine and Remote Sensing, University of South Florida (http://imars.usf.edu, last access: 17 October 2017). These data were processed using the NASA OC4 and OC3 band ratio algorithms (O'Reilly et al., 2000). All products followed the latest implementation of the atmospheric correction based on Ding and Gordon (1995). In situ observations of chlorophyll and nutrients for the Louisiana–Texas shelf were obtained from the Coastal Waters Consortium (CWC; Rabalais, 2015; Smith, 2015; Parson et al, 2015). Chlorophyll observations in the deep GoM were derived from Autonomous Profiling Explorer (APEX) measurements collected during the Lagrangian Approach to Study the Gulf of Mexico Deep Circulation project (Hamilton and Leidos, 2017). Nutrient observations in the deep GoM were obtained from water samples collected in the Gulf of Mexico and East Coast Carbon Cruises (GOMECC; Wanninkhof et al., 2014). Observed primary production rates are derived from measurements collected by Lehrter et al. (2009) in the delta and Texas shelf, and Biggs (1992) and Sanchez (1992) in the deep GoM.

2.2 Model description

We use a 13-component biogeochemical model (hereinafter referred to as GoMBio) that simulates nitrogen (N) and silica (Si) cycling. The model includes nitrate (NO₃), ammonium (NH₄), nanophytoplankton (small phytoplankton, PS), diatom (large phytoplankton, PL), chlorophyll of nanophytoplankton and diatom (ChlS and ChlL), microzooplankton (small zooplankton, ZS), mesozooplankton (large zooplankton, ZL), small and large detritus (DS and DL), opal, labile dissolved organic nitrogen (DON), and silicate (SiOH₄). Small detritus is particulate nitrogen linked to ZS egestion and small plankton (PS + ZS) mortality, while large detritus is particulate nitrogen associated with ZL egestion and large plankton (PL + ZL) mortality. Opal is non-living particulate Si linked to diatom mortality and zooplankton egestion. The state variables NO₃, NH₄, PS, PL, ZS, ZL, DS, DL, and DON are simulated in terms of millimoles of nitrogen per cubic meter (mmol N m⁻³), silicate and opal in terms of millimoles of silicate per cubic meter (mmol Si m⁻³), and ChlS and ChlL in terms of milligrams of chlorophyll per cubic meter (mg chlorophyll m⁻³). The model does not include phosphate as a limiting nutrient for phytoplankton growth. Although previous studies have indicated the existence of phosphate limitation near the MS-A deltas during May–July (Sylvan et al. 2006, 2007; Laurent et al., 2012; Laurent and Fennel, 2014; Fennel and Laurent, 2017), we focus here on the role of N and Si, as observational studies suggest that N and Si can modulate phytoplankton production and composition across the northern GoM (Dortch and Whitledge, 1992; Nelson and Dortch, 1996; Lohrenz et al., 1997; 2008; Rabalais et al., 2002; Zhao and Quigg, 2014).

GoMBio describes the following processes: (1) phytoplankton growth as a function of temperature, light, NO₃, and NH₄, including NH₄ inhibition of NO₃ uptake; (2) silicate limitation of PL growth; (3) photo-acclimation; (4) phytoplankton exudation; (5) ZS grazing on PS and PL; (6) ZL grazing on PS and PL, and predation on ZS; (7) zooplankton egestion and zooplankton excretion; (8) phytoplankton and zooplankton mortality; (9) nitrification; (10) detritus, DON, and opal remineralization; 11) detritus, diatoms, and opal sinking; and (12) sediment coupled nitrification–denitrification (instantaneous remineralization). Processes 1, 3, 9, and 12 follow Fennel et al. (2006, 2011) formulations, while processes 2, 4–8, and 10–11 follow Kishi et al. (2007) formulations. Descriptions of the model equations and parameters are included in the Supplement. Model parameter values are presented in Table 1.

Table 1 Model parameter values.

The model domain encompasses the entire GoM and is based on the Regional Ocean Model System (ROMS) (Shchepetkin and McWilliams, 2005). The model's horizontal resolution is about 8 km and has 37 sigma-coordinate (bathymetry-following) vertical levels. Boundary conditions are Flather (Flather, 1976) and Chapman (Chapman, 1985) for the barotropic velocity and free surface, respectively, and a combination of radiation and nudging for the baroclinic velocity and tracers (Marchesiello et al., 2001). Tidal constituents were not included in the model. The open-boundary nudging timescale is 4 days for the incoming signal and 90 days for the outgoing signal. A third-order upstream scheme and a fourth-order Akima scheme are used for horizontal and vertical momentum advection, respectively. Multidimensional positive definitive advection transport algorithm (MPDATA) is used for horizontal and vertical tracer advection (Smolarkiewicz and Margolini, 1998). Horizontal viscosity and diffusivity are set to 1 m² s⁻¹, increasing gradually to 4 m² s⁻¹ in a 100 km wide sponge layer at the open boundaries to reduce signal reflection problems. The Mellor and Yamada 2.5-level closure scheme is used for vertical turbulence (Galperin et al., 1988). Initial and open-boundary conditions are derived from a 25 km resolution Modular Ocean Model basin-scale model for the Atlantic Ocean (Liu et al., 2015), which includes the Tracers of Ocean Phytoplankton with Allometric Zooplankton (TOPAZ) as biogeochemical model (Dunne et al., 2010). Since TOPAZ does not include zooplankton as a state variable, we assumed zooplankton correspond to 20 % of the total phytoplankton biomass, assigning 30 % to mesozooplankton and 70 % to microzooplankton, assuming that the microzooplankton is the dominant zooplankton component near the model open boundaries for the Gulf of Mexico. Sensitivity simulations indicated that changes to those allocations do not affect greatly the derived plankton biomass patterns.

The model is forced with monthly surface water flux; daily shortwave and longwave radiation; and 6 h resolution air temperature, sea level pressure, humidity, and winds from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis product (0.75^∘ resolution; Dee et al., 2011). Surface net heat flux and wind stress are estimated using bulk parameterization. River runoff from 54 river sources (35 in the US) is explicitly represented. Daily water discharges from US rivers were retrieved from the US Geological Survey (USGS) river gauges (https://waterdata.usgs.gov, last access: 21 August 2017). Climatologies from Mexican river discharges were derived from He et al. (2011), Munoz-Salinas and Castillo (2015), and Martinez-Lopez and Zavala-Hidalgo (2009). Monthly observations of dissolved inorganic nutrients (nitrate, ammonia, silicate) and organic nitrogen in the MS-A rivers were retrieved from the USGS (http://toxics.usgs.gov, last access: 21 August 2017; Aulenbach et al., 2007). Following Yu et al. (2015), the MS-A particulate organic nitrogen (PON) was determined as the difference between unfiltered and filtered total Kjendahl nitrogen (TKN), while the dissolved organic nitrogen (DON) was estimated as the difference between filtered TKN and ammonia. Only 10 % of the estimated DON was incorporated into the model as labile DON, considering that most of the observed MS-A DON corresponds to refractory material (Green et al., 2006). Riverine PON was assigned to the small detritus pool. For river sources other than the MS-A, dissolved inorganic nutrients and organic nitrogen concentrations are prescribed as climatological averages (USGS; Dunn, 1996; He et al., 2011; Livingstone, 2015). Because submarine groundwater discharge (SGD) is a significant source of nitrogen off the west Florida shelf (Hu et al., 2006), we included SGD NH₄ fluxes based on rates reported by Swarzenski et al. (2007). We assumed that SGD NH₄ fluxes occurred in regions shallower than 30 m, decreasing exponentially from 0.694 mmol m⁻² day⁻¹ at 10 m (minimum model depth) to 0.069 mmol m⁻² day⁻¹ at 30 m. Surface photosynthetically active radiation (PAR) is assumed to be 43 % of the surface shortwave radiation. Light attenuation includes a salinity-dependent coefficient (K_salt) as in Fennel et al. (2011).

Figure 1 Model domain and bathymetry. Polygons A, B, and C depict the MS delta, Texas shelf, and deep-ocean region, respectively, selected to describe plankton patterns. Gray contours show the 20, 30, 50, and 200 m isobaths.

A 40-year model spin-up was completed before starting the historical simulation. To spin up the model, we used the basin-model boundary conditions and the ERA-Interim surface fluxes of randomly selected years from 1979 to 2014, following Lee et al. (2011). After spin-up, the model was run continuously from January 1979 until December of 2014, with monthly averaged fields saved.

3 Results

The ocean–biogeochemical model reproduces reasonably well main patterns of temperature, salinity, sea level anomaly, and eddy kinetic and biogeochemical variables. A model–data comparison of selected physical variables is presented in the Supplement. In the following section we perform a validation for chlorophyll, diatom-to-total-chlorophyll ratio, primary production, and nutrients.

Figure 2 Spatial patterns of model and satellite chlorophyll. Comparison between surface chlorophyll concentration (mg m⁻³) derived from model outputs (a, c) and SeaWiFS (b, d) during summer and winter. Gray line depicts the 200 m isobath.

Figure 3 Monthly chlorophyll time series derived from model (GoMBio, red line), SeaWiFS (green line), and MODIS (blue line). Gray area depicts the model chlorophyll range. Correlation coefficient between model and satellite time series is indicated at each panel. Black dots with vertical bars depict the monthly mean and interquartile range of in situ chlorophyll from the Coastal Waters Consortium dataset (a, b) and APEX profiling floats (Lagrangian Approach to Study the Gulf of Mexico Deep Circulation project, c). Ticks on the ordinate mark 1 January of each year.

3.1 Biogeochemical model–data comparison

Modeled surface chlorophyll agreed qualitatively well in the spatiotemporal patterns with the satellite chlorophyll (Fig. 2). The main differences between model and satellite chlorophyll are in the coastal region. Those differences can be explained (in part) by satellite chlorophyll overestimation, due to the high concentration of dissolved colored organic matter and sediments associated with river runoff (Hu et al., 2000; Del Castillo et al., 2001; Gilbes et al., 2002; D'Sa and Miller, 2003). The greatest chlorophyll concentration values are within the MS river delta, and the lowest values within the region influenced by the Loop Current. Significant seasonal differences are evident in the oceanic region, with minimum chlorophyll during summer (June–August), when thermal vertical stratification is the strongest, and maximum chlorophyll during winter (December–February) and early spring (March), concomitant with the greatest surface cooling and wind-driven mixing (Muller-Karger et al., 1991; 2015). To compare temporal patters from model outputs and satellite observations, we derived monthly time series of chlorophyll in three regions: MS delta; Texas shelf and western part of the Louisiana shelf (for simplicity hereinafter referred to as Texas shelf); and the deep-ocean area, encompassing 25–27.5^∘ N, 85.5–92^∘ W (see regions in Fig. 1). The MS delta and the Texas shelf are two productive regions strongly influenced by the MS-A river runoff, whereas the deep-ocean box is an oligotrophic region often influenced by the Loop Current. The simulated chlorophyll time series are strongly correlated with the satellite chlorophyll time series, reproducing main seasonal and interannual patterns (Fig. 3). However, the model tends to underestimate the long-term mean of satellite chlorophyll in the MS delta and Texas shelf (the ratio of satellite chlorophyll to model chlorophyll ranges from 1.98 to 2.80; see Table 2). An underestimation of model chlorophyll is also evident when we contrast the modeled time series with in situ observations from the CWC during spring–summer (black dots in Fig. 3a), although the CWC and simulated chlorophyll tend to agree well during fall and winter, suggesting that satellite sensors could be overestimating surface chlorophyll in these two seasons. This is not surprising for shelf waters influenced by river runoff, as previous studies in the northern GoM have reported that the satellite chlorophyll overestimates in situ chlorophyll by a factor of 2 to 4 (e.g., Nababan et al., 2011). In the oceanic region, the simulated chlorophyll overestimates the long-term mean of SeaWiFS and MODIS chlorophyll by 12 and 22 %, respectively, while in situ chlorophyll estimation based on APEX profiling floats (black dots in Fig. 3c) closely match the model-derived patterns.

Table 2 Long-term mean and standard deviation of model and satellite chlorophyll.

Figure 4 EOF analysis of chlorophyll anomalies. (a, b) First EOF mode of surface model chlorophyll (a) and SeaWiFS chlorophyll (b). (c) Principal component associated with the first EOF mode of model, SeaWiFS, and MODIS chlorophyll. Correlation coefficient between model and satellite PC1 series is indicated in (c).

We evaluated the model's ability to reproduce interannual patterns of chlorophyll by performing empirical orthogonal decomposition of chlorophyll anomaly time series (anomaly refers to monthly outputs/observations with the monthly climatological mean subtracted). The first empirical orthogonal function (EOF1) of model chlorophyll is consistent with the EOF1 from SeaWiFS (Fig. 4a, b) and MODIS (not shown). EOF1 is eminently a coastal pattern, with the greatest values located near the MS-A deltas. The main differences between model and satellite EOF1 are located in the northwestern Florida region, where model chlorophyll is much lower than SeaWiFS chlorophyll, probably linked to a misrepresentation of the interannual variability in riverine nutrient load. The interannual variability of the first principal component (PC1) time series (which represents the temporal variability of EOF1) of model chlorophyll is well correlated to the PC1 time series of SeaWiFS (r=0.66) and MODIS (0.59).

Figure 5 Chlorophyll ratios. (a) Climatological mean of the diatom-to-total-phytoplankton ratio of chlorophyll; (b) comparison between observations and model-derived diatom-to-total-phytoplankton ratio of chlorophyll in the Louisiana shelf (coastal stations A and B depicted in panel a) during April (green) and August (red) of 2010–2012. Vertical and horizontal bars depict ±1 SD. Mean and standard deviation of observed chlorophyll ratio are derived from values reported by Zhao and Quigg (2014).

The model's skill in reproducing the patterns in phytoplankton composition is evaluated through the diatom-to-total-chlorophyll ratios reported by Zhao and Quigg (2014) for two coastal stations off Louisiana (stations A and B in Fig. 5a). The model tends to overestimate the diatom ratio at station A (29.04^∘ N–89.56^∘ W) and underestimate it at station B (28.59^∘ N–92.00^∘ W), but the differences are reasonably small considering the large variability in the observed diatom ratios (Fig. 5b). This variability can be associated with strong mesoscale variability across the MS delta (e.g., Marta-Almeida et al., 2013), which is not reflected in the monthly outputs of our 8 km resolution model. In terms of temporal variability, the model is able to reproduce the observed decline in the diatom ratio during summer.

Model-derived estimations of vertically integrated primary production were compared with observed rates, assuming a carbon-to-nitrogen ratio of 6.625 to express model production in grams of carbon per square meter per day. The temporal variability of the simulated production rates agrees reasonably well with the observed seasonal pattern, though a model underestimation is evident during late summer (Fig. 6a). The interquartile range of model production is 0.87–1.5, 0.32–0.47, and 0.13–0.23 g C m⁻² d⁻¹ for the MS delta, Texas shelf, and deep-ocean region, respectively, which are within the range of production estimated by Lehrter et al. (2009) in the northern shelf and Biggs (1992) and Sanchez (1992) in the deep GoM (Fig. 6b).

Figure 6 (a) Modeled and observed time series of primary production for the Mississippi delta and Texas shelf; gray and light blue shades depict the model production ranges for the Mississippi delta and Texas shelf, respectively; (b) box plots of primary production in the Mississippi delta, Texas shelf, and deep Gulf region derived from observations and model (GoMBio) outputs during spring–summer. Red lines, bottom and top edges of the boxes, and whiskers represent the median, interquartile interval, and non-outlier range, respectively.

Simulated and observed time series of nitrate and silicate at station C6 (28.86^∘ N, 90.46^∘ W; Louisiana shelf) are shown in Fig. 7a, b (observations only available for May–October). The seasonal change in surface nitrate is reproduced well by the model, which displays values >10 mmol m⁻³ during spring and <2 mmol m⁻³ in summer. On the other hand, the simulated surface silicate concentration show a poor agreement with the observed values, which in part could be explained by the relatively weak silicate seasonality and strong mesoscale variability over the Louisiana shelf. To further examine the ability of the model to reproduce coastal patterns in nitrate and silicate concentration on the Louisiana–Texas shelf, we evaluated the relationship between surface salinity and surface nutrient during spring–summer (Fig. 7c, d). Both nitrate and silicate show conservative mixing linked to the Mississippi and Atchafalaya river discharge. The model reproduces well the observed salinity–silicate relationship, while the similarity between the modeled and observed salinity–nitrate relationship is less clear. It is likely that additional observations are required to objectively visualize the observed pattern. However, our simulated salinity–nitrate relationship is consistent with observations by Sylvan et al. (2006) and modeling results by Fennel et al. (2011)(see their Fig. 4).

Figure 7 (a, b) Model–data comparison of nitrate and silicate time series at station C6 (28.86^∘ N, 90.46^∘ W); (c, d) relationship between surface salinity and both nitrate and silicate concentration over the Louisiana–Texas shelf. Modeled values are shown as two-dimensional histogram (color scale), and observations as blue marks.

Nitrate and silicate measurements collected in most oceanic stations of the Mississippi and Tampa lines from GOMECC cruises 1 and 2 were used to evaluate the model's ability to simulate nitrate and silicate patterns in the deep GoM. The modeled nutrient profiles (red lines) reproduce well the depleted nitrate and silicate levels in the upper 30 m, as well as the strong vertical gradient linked to the nutricline over 30–300 m depth (blue dots) (Fig. 8a, d). Some model overestimation of nitrate and silicate is seen at depth >300 m, but that bias most likely has a limited impact on the nutrient concentration in the upper 100 m layer. A better model–observation agreement is observed at the station on the Tampa line.

Figure 8 Comparison between profiles of nitrate and silicate derived from model outputs (red lines) and GOMECC data (blue dots). The model's climatological mean and range for July are also shown as black line and yellow area, respectively. (a, b, e, f) show the profiles associated with the most oceanic station from the Mississippi line and (c, d, g, h) for the Tampa line for GOMECC cruises 1 (July 2007) and 2 (July 2012).

3.2 Phytoplankton biomass patterns

The model–data comparison shown in the previous section, along with the physical model validation presented in the Supplement, indicates that the model is able to reproduce dominant ocean–biogeochemical processes and consequently could be used to explore the underlying factors modulating spatiotemporal changes in diatom and nanophytoplankton biomass. In this section we describe the main seasonal patterns in phytoplankton biomass in the three selected regions shown in Fig. 1. Subsequently we examine the driving factors modulating the phytoplankton biomass seasonality.

The model-derived patterns in plankton biomass have important regional differences in terms of seasonality. To illustrate this, we estimated monthly climatologies of phytoplankton concentration from the surface to 30 m depth (or bottom depth if <30 m) within the MS delta, Texas shelf, and deep-ocean regions (Fig. 9a, c; regions depicted in Fig. 1). Total phytoplankton is the greatest during March–April in the MS delta and Texas shelf, and February–March in the oceanic region, and smallest during August in the three regions. The timing and amplitude of the seasonal maxima differ significantly between phytoplankton components. In the MS delta, diatoms peak in February while nanophytoplankton peak in April–May. In the Texas shelf, the spring phytoplankton maximum is mainly driven by nanophytoplankton. Diatoms do not have a marked spring peak like in the MS delta, displaying two maxima in February (the greatest) and June. In the oceanic region, both nanophytoplankton and diatom peak in February–March, with nanophytoplankton clearly dominating upon diatom (>80 %).

Figure 9 Climatological seasonal cycle of phytoplankton biomass in the 30 m upper layer from the Mississippi delta, Texas shelf, and deep ocean (regions depicted in Fig. 1, gray polygons).

3.2.1 Limitation factors and growth

To investigate the drivers of phytoplankton growth variability, we derive climatological patterns for the nutrient limitation factors (L_P; Eqs. A1.5 and A2.7), the light limitation factors (f_P; Eqs. A1.6 and A2.8), the temperature-dependent growth rates (V_p; Eqs. A1.3 and A2.3), and the specific growth rate (SGR, which is the product of L_P, f_P, and V_p). It is important to note that the nutrient and light limitation factors ranges from 0 to 1, with 0 indicating non-growth and 1 indicating no limitation. This implies that growth limitation is inversely related to the limitation factors. Seasonal changes in L_PS and L_PL for the MS delta, Texas shelf, and deep ocean are depicted in Fig. 10a and b. In the MS delta and Texas shelf, the model nutrient limitation factors are the greatest (i.e., the weakest limitation) during February–April and the smallest (i.e., the strongest limitation) during September–November, reflecting the seasonality in river discharge along the northern shelf (the maximum river discharge in Louisiana and Texas is during April and March, respectively, and the minimum in August–September). A secondary peak in the nutrient limitation factors is observed during July in the Texas shelf, which can be related to wind-driven upwelling and a secondary peak in river discharge during summer. In the deep-ocean region, the nutrient limitation factors are maxima during January–March and minima during July–October, a pattern associated with the seasonal cycle in thermal stratification and mixing (enhanced mixed in winter, enhanced stratification in summer). Significant differences exist between the magnitude of L_PS and L_PL. The L_PS/L_PL ratio is ∼ 1.5 in the MS delta, ∼ 2 in the Texas shelf, and ∼ 3 in the deep GoM. Unlike nanophytoplankton, diatoms can be considerably nutrient limited in the MS delta region. The monthly climatologies of the ratio of silica limitation to nitrogen limitation (SLF : NLF) are used to evaluate whether diatoms are nitrogen limited (SLF : NLF >1) or silica limited (SLF : NLF <1) (Figs. 10c and S10 in the Supplement). Overall SLF : NLF is predominantly >1 in the three regions, implying that diatoms are mainly nitrogen limited. However, SLF : NLF shows values near or smaller than 1 during December–April in the deep Gulf and during February–April in the MS delta, indicating that both nitrogen and silica can limit model diatom growth.

Figure 10 Growth limitation and specific growth rates for nanophytoplankton (PS) and diatoms (PL): (a, b) nutrient limitation factors; (c) ratio of silica limitation to nitrogen limitation (SLF : NLF; for diatoms only); (d, e) light limitation factors; (f, g) temperature-dependent growth; (h, i) specific growth rates. Factors were averaged in the upper 30 m layer from the Mississippi delta, Texas shelf, and deep-ocean regions (depicted in Fig. 1, gray polygons).

Besides nutrients, light and temperature influence model phytoplankton growth. The strongest light limitation is in the MS delta, and the weakest is in the deep GoM (Fig. 10d, e), but the regional differences in light limitation are much smaller than those for nutrient limitation. Seasonally, light limitation is weakest during April in the coastal regions and during May in the deep ocean. Conversely, light limitation is the strongest during August and December in the coastal regions, and during December in the deep ocean. In the coastal regions, the decline in light limitation during June–August can be linked to increased light attenuation, driven by the offshore spread of low-salinity and phytoplankton-rich waters by wind-driven upwelling. The temperature-dependent growth rate (V_p) displays the largest amplitude in the coastal regions, with a maximum in August and minimum in January–February (Fig. 10f, g). The ratio between the maximum and minimum V_p is ∼ 2.3 in the coastal regions and ∼ 1.4 in the deep ocean.

Figure 11 Phytoplankton production, net phytoplankton growth (production minus biological losses), and grazing estimated for the upper 30 m layer of the Mississippi delta (a, d, g), Texas shelf (b, e, h), and deep ocean (c, f, i) (regions depicted in Fig. 1, gray polygons). Grazing terms are microzooplankton upon nanophytoplankton (PS2ZS) and diatoms (PL2ZS), and mesozooplankton upon nanophytoplankton (PS2ZL) and diatoms (PL2ZL).

The interplay among nutrient, light, and temperature conditions determines the model phytoplankton SGR. The seasonal pattern in the SGR shows differences between coastal and oceanic domains (Fig. 10h, i). In the coastal regions, the inverse relationship between V_p and both light and nutrient limitation factors during March–August determines the greatest SGR in June–July, while the small V_p and light limitation factors during December–February determine the minimum SGR in December–January. In the deep-ocean region, the SGR seasonality is mainly driven by nutrient and light limitation. The maximum SGR is in February (concomitant with the maximum nutrient limitation factor), while the minimum SGR is in June–August, the latter driven by the strong nutrient limitation during summer.

3.2.2 Biomass sources and losses

Now we explore how the patterns in phytoplankton production and losses influence the patterns in phytoplankton biomass. We showed that the model SGR is the maximum during June–July in the coastal regions and during February in the deep ocean (Fig. 10h, i). We may expect that the seasonal changes in production reflect the changes in SGR, since production is the product between SGR and phytoplankton biomass. The link between SGR and production is evident in the deep ocean, as SGR and production have maxima in February and minima during July–September (Fig. 11c). However, in the MS delta and Texas shelf, the simulated production peaks occur 2–3 months earlier than the SGR peaks (Fig. 11a, b). This necessarily implies that biomass losses due to biological (grazing, mortality, exudation) and physical (advection/diffusion) processes play an important role in modulating production seasonality during spring–summer.

Figure 12 (a–c) Phytoplankton biomass budget: the advection + mixing term represents the sum of advection and turbulent diffusion of phytoplankton biomass, the net phytoplankton growth is production minus biological losses, and the local rate of change is the balance between net phytoplankton growth and the advection + mixing term. Right y axis (red) is for the local rate of change, and left y axis (blue) is for the net phytoplankton growth and the advection + mixing term. (d–f) Components of the advection + mixing term: Hadv, Vadv, and Vmix correspond to horizontal advection, vertical advection, and vertical mixing, respectively. Horizontal mixing can be neglected in the budget analysis, as it is 2 orders of magnitude smaller than other physical term components. Patterns are averages within the upper 30 m ocean layer from the Mississippi delta, Texas shelf, and deep ocean (regions depicted in Fig. 1).

To evaluate how biologically driven processes influence the seasonal patterns in model phytoplankton biomass, we calculated the net phytoplankton growth, which is the balance between production and biological losses (Fig. 11d, f). The net phytoplankton growth displays distinct patterns for each phytoplankton component and region. The maximum net growth for diatoms is in January–February in the MS delta, December–January in the Texas shelf, and February in the deep ocean, while the maximum net growth for nanophytoplankton is in April in the MS delta, February in the Texas shelf, and January in the deep ocean. The net growth for diatoms and nanophytoplankton begins to decline before the production maximum. Moreover, in the Texas shelf, the net growth is negative during the production maximum. In the three regions, the net growth for total phytoplankton (diatoms plus nanophytoplankton) is positive in November–February, has a marked decline in spring, and is negative in May–August. The seasonality of the net phytoplankton growth contrasts with the pattern in the SGR in the MS delta and Texas shelf, as SGR is minimum in December–January and maximum in June. All these features suggest that the seasonal changes in model phytoplankton biomass are strongly modulated by biological losses. Zooplankton grazing is the dominant biological loss term (Fig. 11g, i), markedly prevailing upon mortality and exudation (not shown). Microzooplankton exert the strongest grazing pressure on nanophytoplankton biomass, and mesozooplankton on diatoms, with the grazing patterns closely following the patterns in production. The seasonal patterns for microzooplankton (mesozooplankton) grazing upon nanophytoplankton (diatoms) closely follow the patterns in nanophytoplankton (diatom) production. Peaks in micro- and mesozooplankton grazing are concomitant or lag by 1 month the peak in nanophytoplankton and diatom production.

The seasonal patterns in net phytoplankton growth do not completely explain the seasonal changes in model phytoplankton biomass. To fully elucidate the local phytoplankton biomass change, the role of physically driven fluxes of phytoplankton biomass needs to be examined. To this effect, we estimate the advection + mixing term, which represents the sum of advection and turbulent diffusion of phytoplankton biomass, and compare it with the net phytoplankton growth (Fig. 12a–c). The balance between these two terms determines the local rate of change of phytoplankton biomass. The net phytoplankton growth is generally inversely related to the advection + mixing term, implying that the biologically driven changes tend to be offset by the physically driven changes. Furthermore, the net phytoplankton growth is generally larger than the advection + mixing term, and consequently the sign of the local rate of change is mainly determined by the biological component. The few exceptions are the positive growth during September in the MS delta and during April and September in the Texas shelf, and the negative growth in March–April in the deep-ocean region. In the last case, the physically driven fluxes influence not only the amplitude of the monthly biomass change but also the timing of the seasonal maxima. In the MS delta, the greatest magnitude for the advection + mixing term is during January–April, representing biomass losses mostly linked to horizontal advection (Fig. 12d). The advection can be related to the downstream export of phytoplankton-rich water associated with the MS river plume. A substantial fraction of phytoplankton biomass from the MS-A delta is transported to the Texas shelf, which explains the positive advection + mixing term during March–June (Fig. 12b, e). In the deep ocean, the greatest magnitude for the advection + mixing term is in December–February, representing biomass losses mainly driven by turbulent vertical diffusion (Fig. 12c, f). The close similitude between the magnitude of the advection + mixing term and the net phytoplankton growth determines a much smaller local rate of change in the deep ocean than in the coastal regions (about 1 order of magnitude).

4 Discussion

We configured an ocean–biogeochemical model for the GoM that explicitly represents two types of phytoplankton and zooplankton, and nitrogen and silica as limiting nutrients for phytoplankton growth. Our model reproduces reasonably well the main physical and biochemical patterns, although an underestimation of the mean surface chlorophyll is evident in the northern shelf, especially at bottom depth <20 m. A comparison with in situ chlorophyll observations suggests that part of the model–satellite chlorophyll disagreement could be linked to chlorophyll overestimation by satellite sensors during fall–winter. Realistic representations of phytoplankton variability in regions with strong physical and biochemical gradients, like those in the northern GoM, are challenging. Previous modeling efforts on the Louisiana–Texas shelf based on Fennel's model reproduced better the mean satellite chlorophyll condition than our model (e.g., Fennel et al., 2011; Laurent et al., 2012). However, Fennel's model tends to overestimate satellite chlorophyll by a factor >3 in the deep-ocean region during winter, which could be linked to misrepresentation of microzooplankton grazing (see Sect. 4 in the Supplement). We acknowledge that additional components and processes could be included in our model – such as phosphorus cycling, iron limitation, and nitrogen fixation – to represent more realistic biogeochemical dynamics. We also recognize that more observational studies will be required to constrain better our model parameters, as well as the biogeochemical fluxes between land and ocean. Nevertheless, we believe that the current model configuration can capture well enough the seasonal dynamics of diatoms and nanophytoplankton biomass in the GoM. It is known that variations in phytoplankton composition can have important repercussion for the ecosystem, including changes in upper-trophic-level dynamics, carbon export (carbon export is enhanced in diatom-dominated food webs), and bottom hypoxia (Dagg et al., 2003; Green et al., 2008). Therefore, modeling efforts exploring variability in phytoplankton components, such as this study, are needed to advance our understanding of ecosystem variability in the GoM.

We examined the main model phytoplankton biomass patterns and explored the underlying factors explaining biomass variability, following a similar approach to that used by Fennel et al. (2011). We used a constant-depth layer (0–30 m), whereas Fennel et al. (2011) calculated seasonal patterns in a seasonally variable mixed-depth layer (∼ 10 m in summer to ∼ 40 m in winter). We chose a constant-depth layer because it makes the biomass budget analysis more straightforward. It is also worthwhile to mention that an important fraction of primary production can be distributed below the mixed layer in spring–summer (Yu et al., 2015). Our growth limitation analysis compared distinct regions in terms of phytoplankton production and river runoff influence, including the oligotrophic deep GoM, a region that has received less attention in previous modeling studies. We found that nutrient limitation displayed the largest spatial differences compared to other limiting factors (light and temperature). Although the model indicated that the main limiting nutrient for model diatom is nitrogen, silicate also can limit model diatom production in the deep GoM during winter and in the MS delta during spring. The latter agrees with observations of severe silica depletion during spring in the Louisiana shelf (Dortch and Whitledge, 1992; Nelson and Dortch, 1996). Although observational studies suggested the occurrence of silica limitation in the MS delta decades ago with a potential link to anthropogenically driven declines in the MS river Si : N ratio (Turner and Rabalais, 1991), this is the first modeling study to evaluate the role of silica as a driver of diatom growth in the region. The implication for silica and nitrogen limitation in the Louisiana–Texas shelf is that changes in the MS-A river nutrient load can modulate changes in diatom production, influencing phytoplankton composition.

The simulated SGR patterns showed important difference between coastal and oceanic domains. Nutrients, light, and temperature are important in modulating the seasonal SGR changes on the northern shelf, while nutrients and light are the dominant factors driving the SGR seasonality in the deep GoM. The monthly averages for the SGR in small and large phytoplankton range within 0.28–0.85 and 0.18–0.57 day⁻¹ in the coastal regions, with the maximum (minimum) values in June–July (December–January). These SGR values are within the observational range reported by Fahnestiel et al. (1995) and similar to model estimations by Fennel et al. (2011). In the oceanic region, the SGR range for nanophytoplankton is 0.17–0.40 day⁻¹, with the maximum (minimum) values in February–March (June–September). Consistent with Fennel et al. (2011), we found that zooplankton grazing plays a leading role in modulating phytoplankton biomass seasonality. This is especially evident in the coastal regions, where the net phytoplankton growth is negative (biomass decrease) in summer and positive (biomass increase) in winter, i.e., opposite to the pattern in the SGR.

Our study examined the coupled role of biologically (production and biological losses) and physically (advection and vertical mixing) driven biomass fluxes. Previous studies suggested the importance of advection and diffusion as drivers of biomass changes in the GoM (e.g., Dagg et al., 2003; Green et al., 2008; Fennel et al., 2011). However, a quantification of these dynamics in biogeochemical model has not been done in the region. We found that the seasonal patterns in model phytoplankton biomass are largely determined by small imbalances between biologically and physically driven fluxes, the latter mainly associated with horizontal advection in the Louisiana–Texas shelf and turbulent vertical diffusion in the deep GoM. Consequently, we cannot obtain a proper understanding of biomass seasonality when the physically driven biomass fluxes are excluded from the analysis. Disentangling the processes influencing phytoplankton seasonality is a complex task, as the mechanisms acting as physical loss terms can also influence the balance between production and biological losses. That is the case for turbulent vertical diffusion, which modulates the vertical distribution of nutrients (impacting phytoplankton production) and zooplankton (impacting zooplankton grazing) (Behrenfeld, 2010).

Finally, future projections of environmental scenarios suggest substantial increases in both river runoff and thermal stratification in the northern GoM due to anthropogenic climate change (Tao et al., 2014; Liu et al., 2015). Therefore, how such environmental disturbances acting at multiple timescales can alter the subtle imbalances between primary production and biological losses (or between biological and physical driven biomass fluxes) is a topic that deserves further attention.

5 Summary and conclusions

A coupled ocean–biogeochemical model was configured for the GoM to examine underlying mechanisms determining spatial and seasonal variability in diatoms and nanophytoplankton biomass. We investigated the factors modulating the specific growth rate (SGR) and explored the seasonal changes in biologically and physically driven biomass fluxes. We found that model diatom growth was ∼ 40 % and >80 % nutrient limited in the Louisiana shelf and deep GoM, respectively, whereas model nanophytoplankton growth was ∼ 10 % and 40–85 % nutrient limited. Our model indicates that diatom growth is mainly limited by nitrogen. However, silica limitation can occur in the deep GoM during winter, and in the MS delta during spring. The interplay among nutrient, light, and temperature determined the SGR seasonal timing (max/min) in the Louisiana–Texas shelf, while nutrient and light determined the simulated SGR seasonal timing in the deep GoM. Primary production in the model was not only driven by changes in SGR but also influenced by biomass losses linked to zooplankton grazing. Moreover, the net phytoplankton growth (i.e., the balance between primary production and biological losses) revealed top-down control of phytoplankton biomass. The physically driven biomass fluxes, mainly associated with horizontal advection in the Louisiana–Texas shelf and turbulent vertical diffusion in the deep GoM, played a key role in modulating amplitude and phase in the seasonal phytoplankton biomass cycle. These results stress the importance of an integrated analysis of biologically and physically driven biomass fluxes to better characterize phytoplankton biomass phenologies.