Hydrodynamic and biochemical impacts on the development of hypoxia in the Louisiana–Texas shelf – Part 1: roles of nutrient limitation and plankton community

A three-dimensional coupled hydrodynamic–biogeochemical model with multiple nutrient and plankton functional groups was developed and adapted to the Gulf of Mexico to investigate the role of nutrients and the complexity of plankton community in dissolved oxygen (DO) dynamics. A 15-year hindcast was achieved covering the period of 2006–2020. Extensive model validation against in situ data demonstrates that the model was capable of reproducing vertical distributions of DO and spatial distributions of bottom DO concentration, as well as their interannual variations. The study demonstrates that bottom DO dynamics and hypoxia evolution are significantly influenced by both physical processes and local biochemistry, with sedimentary oxygen consumption and vertical diffusion identified as key contributors. Summer hydrodynamics play a critical role in nutrient distribution and limitation: a notable expansion of Si limitation was simulated when coastal currents shifted eastward or northward. This effect, especially pronounced on the western part of the Louisiana–Texas shelf, underscores the importance of nutrient limitation in shaping DO dynamics. The model identifies a bi-peak primary production pattern in spring and early summer, aligned with satellite chlorophyll a variations, attributed to the complexity of the plankton community and interactions among different plankton groups. Our findings emphasize the necessity of integrating sophisticated plankton community dynamics into biogeochemical models to understand primary production variability and its impact on bottom hypoxia.

1 Introduction

The Louisiana–Texas (LaTex) shelf in the northern Gulf of Mexico (nGoM) has one of the most notorious recurring hypoxia zones in the world (bottom dissolved oxygen (DO) < 2 mg L⁻¹, Rabalais et al., 2002; Rabalais et al., 2007a; Justić and Wang, 2014). Historical observations show that hypoxia usually emerges in mid-May and persists through mid-September (Rabalais et al., 1999, 2002). The hypoxic zone can cover as much as 23 000 km² and has a volume of up to 140 km³ (Rabalais and Turner, 2019; Rabalais and Baustian, 2020). Although nitrogen (N) is the ultimate limiting nutrient, phosphorus (P) load reduction also leads to a significant reduction of the hypoxia area (Fennel and Laurent, 2018). Transient P limitation on the shelf (Laurent et al., 2012; Sylvan et al., 2007) was deemed to be associated with the delayed onset and reduction of the hypoxia area. Sensitivity experiments with hypoxia area reduction to different nutrient reduction strategies suggested that to meet the hypoxic area reduction goal (< 5000 km² in a 5-year running average) set by the Mississippi River/Gulf of Mexico Watershed Nutrient Task Force (2008), a dual-nutrient strategy with a reduction of 48 % of total N and inorganic P would be the most effective way (Fennel and Laurent, 2018).

Coastal eutrophication in the LaTex shelf leads to a high rate of microbial respiration and depletion of DO (Rabalais et al., 2007b). Incubation studies in the LaTex shelf have suggested that sediment oxygen consumption (SOC) accounted for 20 ± 4 % (Murrell and Lehrter, 2011) to 25 ± 5.3 % (McCarthy et al., 2013) of below-pycnocline respiration, nearly 7-fold greater than the corresponding percentage in waters overlying sediments (3.7 ± 0.8 %, about 20 cm above sediments in McCarthy et al., 2013). The numerical study by Fennel et al. (2013) calculated the corresponding SOC fraction, which reached 60 % when applying the water respiration rates of Murrell and Lehrter (2011) and sediment respiration rates of Rowe et al. (2002). Another numerical study (Yu et al., 2015) also pointed out that on the LaTex shelf, oxygen consumption at the bottom water layer was more associated with SOC than water column respiration. According to in situ data and statistical analysis, SOC can be estimated using the bottom temperature and DO concentration (e.g., Hetland and DiMarco, 2008). Nevertheless, many numerical studies have treated SOC as only associated with the abundance of organic matter in the sediment (e.g., Justić and Wang, 2014; Fennel et al., 2006, 2011). An instantaneous remineralization parameterization by Fennel et al. (2006, 2011) estimated SOC as a function of sediment detritus and phytoplankton. Using this scheme, Große et al. (2019) found that the simulated SOC was supported by Mississippi N supply (51 ± 9 %), Atchafalaya N supply (33 ± 9 %), and open-boundary N supply (16 ± 2 %). However, the instantaneous remineralization parameterization tends to overestimate SOC at the peak of phytoplankton blooms while underestimating SOC after the blooms. In a realistic environment, there should be a lag between the blooms and the peak SOC (Fennel et al., 2013). Developments of coupled sediment–water models have emphasized the importance of biogeochemical processes in sediments for SOC dynamics and evolution of bottom hypoxia in the shelf (Moriarty et al., 2018; Laurent et al., 2016). However, coupled sediment–water models are computationally more expensive than a simplified parameterization of SOC. Especially for long-term simulations and time-sensitive forecasts, it is crucial to balance the model's efficiency with its complexity.

In addition to SOC and excess nutrient supply from the rivers, water column stratification also plays an important role in regulating the variability of bottom DO concentration in the LaTex shelf. Strong stratification prohibits DO ventilation and thus reduces DO supply to the bottom water layer (Hetland and DiMarco, 2008; Bianchi et al., 2010; Fennel et al., 2011, 2013, 2016; Justić and Wang, 2014; Wang and Justić, 2009; Feng et al., 2014; Yu et al., 2015; Laurent et al., 2018). On the shelf, the Mississippi plume and the Atchafalaya plume introduce buoyancy, leading to a stable water column and weak DO ventilation processes (Mattern et al., 2013; Fennel and Testa, 2019). Due to the different distances from major river mouths, the influence of freshwater-induced buoyancy varies along the shelf. Moreover, the transport and deposition processes of organic matter are affected by the coastal alongshore current systems, resulting in an SOC gradient across the shelf. For instance, Hetland and DiMarco (2008) pointed out that in the west of Terrebonne Bay, where stratification is usually weak, bottom hypoxia is mainly controlled by bottom respiration.

The phytoplankton blooms on the LaTex shelf mainly result from cyanobacteria and diatoms (Wawrik and Paul, 2004; Schaeffer et al., 2012; Chakraborty et al., 2017). In the Mississippi River plume, diatoms were found to be the most diverse algal class, accounting for over 42 % of all unique genotypes observed (Wawrik and Paul, 2004). Cruise data in the nGoM indicated that diatoms accounted for ∼ 50 % to ∼ 65 % (inner-shelf) and ∼ 33 % to ∼ 64 % (mid-shelf) of chlorophyll a in winter and spring and ∼ 30 % to ∼ 46 % (inner-shelf) during summer and fall, respectively (Chakraborty and Lohrenz, 2015). A field survey documented the fact that the biovolume contribution of diatoms to the total phytoplankton could have been as high as 80 % and 70 % during the upwelling seasons in 2013 and 2014, respectively (Anglès et al., 2019). While a lot of existing studies indicate that N and P were more limited than silicon (Si) on the shelf (e.g., for cruises in 2004 in Quigg et al., 2011; for cruises in 2012 in Zhao and Quigg, 2014; for cruises in 1984, 1994, 2005, 2010, and 2011 in Turner and Rabalais, 2013), Si limitation has also been reported in both plume and shelf water. A bioassay study on sampled collected in spring and summer 2004 showed signs of co-limitation of N, P, and Si at multiple sites (Quigg et al., 2011). Based on cruise studies in the plume of the Mississippi River in 1992 and 1993, strong Si limitation in spring was found due to the increasing N : Si ratio in the Mississippi River water (Nelson and Dortch, 1996). Cruise measurements in 1987 and 1988 also suggested the likelihood of Si limitation, which sometimes overwhelmed the N limitation (Dortch and Whitledge, 1992).

Numerical studies for hypoxia in the LaTex shelf were developed mostly incorporating nutrient flows of N and P only (e.g., Fennel et al., 2006, 2011, 2013; Laurent et al., 2012; Laurent and Fennel, 2014; Fennel and Laurent, 2018; Justić et al., 2003, 2007; Justić and Wang, 2014; Große et al., 2019; Moriarty et al., 2018). In addition, many existing models utilized an oversimplified lower-trophic-level model (one phytoplankton + one zooplankton function group or only one phytoplankton group). The recycling of nutrients in water columns and the associated biogeochemical processes, which may be important to hypoxia evolution (e.g., in the Chesapeake Bay by Testa and Kemp, 2012), could be oversimplified. Moreover, we noticed that there was a bi-peak primary production pattern observed by satellite and modeled by Gomez et al. (2018) (see comparisons of modeled and satellite-derived chlorophyll a concentration in that work). Their biogeochemical model incorporated a more complex community (two phytoplankton + three zooplankton function groups) than other oversimplified models where the bi-peak pattern was hardly captured (e.g., Fennel et al., 2011). The temporal variation of shelf primary production can further induce corresponding changes in DO concentration and in the bottom hypoxia. In this study, we aimed to investigate the possible Si limitation and to assess the impacts of the complexity of the plankton community on DO dynamics and bottom hypoxia development. We adapted and modified a coupled physical–biogeochemical model covering the entire Gulf of Mexico (GoM) by introducing the oxygen and P cycles to the North Pacific Ecosystem Model for Understanding Regional Oceanography (NEMURO, Kishi et al., 2007). The model has two phytoplankton and three zooplankton functional groups for a more comprehensive representation of the plankton community. We also modified the instantaneous remineralization parameterization by adding a conceptual sedimentary organic pool (represented by a sedimentary particulate organic N pool, PON_sed; Fig. 1) to allow the accumulation of organic matter in the sediment. The influence of the community is represented in the biogeochemical processes in water columns and sediments and will eventually be reflected in the bottom DO variability.

2 Methods

2.1 Coupled hydrodynamic–biogeochemical model

We adapted the three-dimensional, free-surface, topography-following community model, the Regional Ocean Model System (ROMS, version 3.7), on the platform of Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) modeling system (Warner et al., 2010) to the GoM (Gulf–COAWST). ROMS solves finite-difference approximations of Reynolds-averaged Navier–Stokes equations by applying hydrostatic and Boussinesq approximations with a split explicit time-stepping algorithm (Haidvogel et al., 2000; Shchepetkin and McWilliams, 2005, 2009). The biogeochemical model applied is primarily based on NEMURO developed by Kishi et al. (2007). NEMURO is a concentration-based, lower-trophic-level ecosystem model developed and parameterized for the North Pacific. The original NEMURO has 11 concentration-based state variables, including nitrate (NO₃), ammonium (NH₄), small and large phytoplankton biomass (PS and PL), microzooplankton biomass, mesozooplankton biomass, predatory zooplankton biomass (ZS, ZL, and ZP), particulate and dissolved organic N (PON and DON), particulate silica (Opal), and silicic acid (Si(OH)₄). NEMURO is known for its capability to distinguish ZS, ZL, and ZP and to provide a detailed analysis of the dynamics of different functional groups. It was widely used in studies of plankton biomass on regional scales (Fiechter and Moore 2009; Gomez et al., 2018; Shropshire et al., 2020). The embedded Si cycle permits the inclusion of a diatom group (i.e., PL), one of the dominant phytoplankton groups in the LaTex shelf.

2.2 Model modification

In a recent effort, Shropshire et al. (2020) adapted and modified NEMURO to the GoM with five structural changes. (1) The grazing pathway of ZL on PS was removed since, in the GoM, the PS group is predominated by cyanobacteria and picoeukaryotes, which are too small for direct feeding by most mesozooplankton (i.e., ZL). (2) A linear function of mortality was applied for PS, PL, ZS, and ZL, while quadratic mortality was used for ZP, accounting for predation pressure of unmodeled predators, like planktivorous fish. (3) The ammonium inhibition term in the nitrate limitation function was no longer considered exponentially but followed the parameterization by Parker (1993). (4) Light limitation on photosynthesis was replaced with the Platt et al. (1980) functional form, which was also implemented in the newer version of NEMURO. (5) A constant C : Chl ratio was replaced with a variable C : Chl model according to the formulation by Li et al. (2010).

Neither the modified (Shropshire et al., 2020) nor the original (Kishi et al., 2007) NEMURO considered P and oxygen cycles. In this study, we introduced a P cycle into NEMURO, including three concentration-based state variables: phosphate (PO₄), particulate organic P (POP), and dissolved organic P (DOP). The P limitation on phytoplankton growth was introduced using the Michaelis–Menten formula. In NEMURO, N serves as the common “currency” when measuring the plankton concentration (mmol N m⁻³). In the river-dominated LaTex shelf, rivers supply inorganic and organic nutrients. In our model, riverine PO₄ (Fig. C1c), DOP, and POP were prescribed based on water quality measurements at river gages. When no measurement was available, the PO₄, DOP, and POP were approximated using total nitrate + nitrite (NO₃ + NO₂), dissolved organic N (DON), and particulate organic N (PON) measurements, respectively, via the Redfield ratio of P : N $=\mathrm{1}:\mathrm{16}$. We neglected the POP settling process but preserved these pools by introducing the stoichiometric ratio between P and N instead. In other words, the sinking process of POP is implicitly included by building linkages between PON and POP concentrations, as the sinking of PON is considered in the model. Governing equations for P state variables are given according to Eqs. (1)–(3). Please also refer to the appendices for more details on expressions of modified terms (Appendix A), state variables (Appendix Table B1), source and sink terms (Appendix Table B2), and values of parameters (Appendix Table B4).

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{ll}\frac{\mathrm{d}\left({\mathrm{PO}}_{\mathrm{4}}\right)}{\mathrm{d}t}& =\left(\mathrm{ResPSn}+\mathrm{ResPLn}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & +\left(\mathrm{DecP}\mathrm{2}\mathrm{N}+\mathrm{DecD}\mathrm{2}\mathrm{N}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & +\left(\mathrm{ExcZSn}+\mathrm{ExcZLn}+\mathrm{ExcZPn}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & -\left(\mathrm{GppPSn}+\mathrm{GppPLn}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\end{array}\text{(2)}& {\displaystyle }\begin{array}{ll}\frac{\mathrm{d}\left(\mathrm{DOP}\right)}{\mathrm{d}t}& =\left(\mathrm{DecP}\mathrm{2}\mathrm{D}-\mathrm{DecD}\mathrm{2}\mathrm{N}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & +\left(\mathrm{ExcPSn}+\mathrm{ExcPLn}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\end{array}\text{(3)}& {\displaystyle }\begin{array}{ll}\frac{\mathrm{d}\left(\mathrm{POP}\right)}{\mathrm{d}t}& =(\mathrm{MorPSn}+\mathrm{MorPLn}+\mathrm{MorZSn}\\ & +\mathrm{MorZLn}+\mathrm{MorZPn})\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & +\left(\mathrm{EgeZSn}+\mathrm{EgeZLn}+\mathrm{EgeZPn}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\\ & -\left(\mathrm{DecP}\mathrm{2}\mathrm{N}+\mathrm{DecP}\mathrm{2}\mathrm{D}\right)\cdot \mathrm{RPO}\mathrm{4}\mathrm{N}\end{array}\end{array}$$

We further adapted the oxygen cycle developed by Fennel et al. (2006, 2013) to NEMURO for hypoxia simulations. However, our model's biogeochemical processes are slightly different due to the different plankton functional groups considered. Sources for oxygen are contributed by the photosynthesis of two phytoplankton functional groups. In comparison, the sinks are attributed to respirations of two phytoplankton functional groups, metabolism of three zooplankton functional groups, light-dependent nitrification (Olson, 1981; Fennel et al., 2006), aerobic decomposition of particulate and dissolved organic matter (measured as PON, and DON, respectively), and SOC. Oxygen air–sea flux was estimated following parameterizations by Wanninkhof (1992). The biogeochemical dynamics of oxygen were adopted as follows (Eq. 4; also see detailed descriptions of variables and parameters in Appendix A–B).

$$\begin{array}{}\text{(4)}& \begin{array}{ll}\frac{\mathrm{d}\left(\mathrm{Oxyg}\right)}{\mathrm{d}t}& =\left({\mathrm{rOxNO}}_{\mathrm{3}}\cdot \mathrm{GppNPS}+{\mathrm{rOxNH}}_{\mathrm{4}}\cdot \mathrm{GppAPS}\right)\\ & +\left({\mathrm{rOxNO}}_{\mathrm{3}}\cdot \mathrm{GppNPL}+{\mathrm{rOxNH}}_{\mathrm{4}}\cdot \mathrm{GppAPL}\right)\\ & -\mathrm{ResPSn}\cdot [\mathrm{RnewS}\cdot {\mathrm{rOxNO}}_{\mathrm{3}}\\ & +\left(\mathrm{1}-\mathrm{RnewS}\right)\cdot {\mathrm{rOxNH}}_{\mathrm{4}}]\\ & -\mathrm{ResPLn}\cdot [\mathrm{RnewL}\cdot {\mathrm{rOxNO}}_{\mathrm{3}}\\ & +\left(\mathrm{1}-\mathrm{RnewL}\right)\cdot {\mathrm{rOxNH}}_{\mathrm{4}}]\\ & -{\mathrm{rOxNH}}_{\mathrm{4}}\cdot \left(\mathrm{ExcZSn}+\mathrm{ExcZLn}+\mathrm{ExcZPn}\right)\\ & -\mathrm{2}\cdot \mathrm{Nit}\cdot \mathrm{LgtlimN}\cdot \widehat{r}\\ & -{\mathrm{rOxNH}}_{\mathrm{4}}\cdot (\mathrm{DecD}\mathrm{2}\mathrm{N}+\mathrm{DecP}\mathrm{2}\mathrm{N})\cdot \widehat{r}\\ & -\mathrm{SOC}\cdot {\mathrm{THK}}_{\mathrm{bot}}\end{array}\end{array}$$

A PON_sed pool due to vertical sinking processes of PON was introduced for parameterization of SOC. The SOC scheme (Fennel et al., 2006) is known as the instantaneous consumption of DO. As soon as the PON falls into the sediment bed, PON will be decomposed instantaneously. This scheme tends to overestimate SOC at the peak of blooms and to underestimate SOC after blooms since the lag in SOC demand is neglected (Fennel et al., 2013). We considered such temporal delays in SOC by introducing a PON_sed pool. A portion of the PON ends with PON_sed, while the rest is buried (PON_burial) and removed from the system. The parameterization is shown in the following. (1) Organic matter settling down at the conceptual sediment layer is remineralized at a temperature-dependent aerobic remineralization rate, K_P2N. (2) Sediment oxygen is consumed only in the oxidation of sedimentary organic matter (represented by PON_sed) and the nitrification of ammonium to nitrate (Fennel et al., 2006). (3) Oxygen consumption at the conceptual sediment layer directly contributes to oxygen concentration decreases only at the bottom water column. (4) Sediment denitrification is linearly related to SOC according to observational-based estimates by Seitzinger and Giblin (1996), but the relationship was modified by Fennel et al. (2006) with a slightly smaller slope of denitrification on SOC rate, i.e.,

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathrm{denitrification}\left({\mathrm{mmol}\mathrm{N}\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}\right)\\ & =\mathrm{0.105}\times \mathrm{SOC}\left({\mathrm{mmol}\mathrm{O}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}\right).\end{array}\end{array}$$

(5) Aerobic decomposition of PON_sed, sediment nitrification, and denitrification follow chemical equations according to (Fennel et al., 2006)

$$\begin{array}{}\text{(R1)}& {\displaystyle }\begin{array}{rl}{\mathrm{C}}_{\mathrm{106}}{\mathrm{H}}_{\mathrm{263}}{\mathrm{O}}_{\mathrm{110}}{\mathrm{N}}_{\mathrm{16}}\mathrm{P}& +\mathrm{106}{\mathrm{O}}_{\mathrm{2}}\leftrightarrow \mathrm{106}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{16}{\mathrm{NH}}_{\mathrm{4}}\\ & +{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}+\mathrm{122}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\end{array}\text{(R2)}& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{3}}+\mathrm{2}\mathrm{H}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(R3)}& {\displaystyle }\begin{array}{rl}{\mathrm{C}}_{\mathrm{106}}{\mathrm{H}}_{\mathrm{263}}{\mathrm{O}}_{\mathrm{110}}{\mathrm{N}}_{\mathrm{16}}\mathrm{P}& +\mathrm{84.8}{\mathrm{HNO}}_{\mathrm{3}}\to \mathrm{106}{\mathrm{CO}}_{\mathrm{2}}\\ & +\mathrm{42.4}{\mathrm{N}}_{\mathrm{2}}+\mathrm{16}{\mathrm{NH}}_{\mathrm{3}}+{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\\ & +\mathrm{148.4}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}\end{array}$$

(6) Nitrate produced in sediments (Reaction R2) is used for denitrification (Reaction R3). The linear assumption of denitrification and SOC (Eq. 5) implicitly builds relationships among the aerobic decomposition of PON_sed (Reaction R1), sediment nitrification (Reaction R2), and denitrification (Reaction R3). Let us assume that the production rate of NH₄ by aerobic decomposition (Reaction R1) of organic matter is M mmol m⁻³ d⁻¹ and that the fraction of denitrification-produced CO₂ (Reaction R3) to the total CO₂ production (Reactions R1 and R3) is x. According to the linear assumption mentioned above, the consumption rate of NO₃ during denitrification (Reaction R3) is proportional to the total consumption rate of O₂ in the sediment (Reactions R1 and R2), yielding $\frac{\mathrm{84.8}Mx}{\mathrm{16}\left(\mathrm{1}-x\right)}=\mathrm{0.105}\times \left[\frac{\mathrm{106}M}{\mathrm{16}}+\frac{\mathrm{84.8}Mx}{\mathrm{8}\left(\mathrm{1}-x\right)}\right]$ and further x≈0.1425. The oxygen consumption rate (Eq. 6) and organic matter consumption rate (Eq. 7) due to the coupled aerobic decomposition, nitrification, and denitrification processes can be obtained by substituting the x value into the stoichiometric ratios according to Reactions (R1)–(R3).

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathrm{Oxyg}}_{\mathrm{consumption}}={\displaystyle \frac{\mathrm{106}M}{\mathrm{16}}}+{\displaystyle \frac{\mathrm{84.8}Mx}{\mathrm{8}\left(\mathrm{1}-x\right)}}=\mathrm{8.3865}M\text{(7)}& {\displaystyle }{\mathrm{OM}}_{\mathrm{consumption}}={\displaystyle \frac{M}{\mathrm{16}}}+{\displaystyle \frac{Mx}{\mathrm{16}\left(\mathrm{1}-x\right)}}=\mathrm{0.0729}M\end{array}$$

Accordingly, the SOC and consumption rate of PON_sed are given, respectively, as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{SOC}={\mathrm{Oxyg}}_{\mathrm{consumption}}\cdot {\mathrm{THK}}_{\mathrm{bot}}=\mathrm{8.3865}M\cdot {\mathrm{THK}}_{\mathrm{bot}},\text{(9)}& {\displaystyle }\begin{array}{rl}{\mathrm{PON}}_{{\mathrm{sed}}_{\mathrm{consumption}}}& =\mathrm{16}\cdot {\mathrm{OM}}_{\mathrm{consumption}}\cdot {\mathrm{THK}}_{\mathrm{bot}}\\ & =\mathrm{1.1662}M\cdot {\mathrm{THK}}_{\mathrm{bot}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\displaystyle }M={\displaystyle \frac{{\mathrm{PON}}_{\mathrm{sed}}\cdot \mathrm{VP}\mathrm{2}{\mathrm{N}}_{\mathrm{0}}\cdot \exp \left(K_{\mathrm{P}\mathrm{2}\mathrm{N}}\cdot \mathrm{TMP}\right)}{{\mathrm{THK}}_{\mathrm{bot}}}},\text{(11)}& {\displaystyle }{\mathrm{THK}}_{\mathrm{bot}}=\text{thickness of bottom water column}.\end{array}$$

We further added light inhibition to nitrification and aerobic decomposition. These parameterizations were applied following descriptions by Fennel et al. (2006, 2013). For the oxygen-dependent term, an oxygen threshold is specified below which no aerobic respiration or nitrification occurred. Detailed equations are listed in Appendix A. The structure of the newly modified NEMURO is shown in a schematic diagram in Fig. 1.

Figure 1 Schematic diagram of the modified NEMURO. Note that the P flow and the oxygen flow are two newly added flows to the original NEMURO.

2.3 Model setups

The coupled model was applied to the GoM using an Arakawa C-grid with a horizontal resolution of ∼ 5 km (Fig. 2a). There are 334 and 357 interior ρ points in the east–west and north–south directions, respectively. The model includes 36 σ layers vertically. The wetting and drying scheme (Warner et al., 2013) was implemented to provide a more accurate representation of shallow water. The computational time step (i.e., baroclinic time step) was set to 240 s, while the number of barotropic time steps between each baroclinic time step was set to 30. Model hindcast was carried out from 1 August 2006 to 26 August 2020, with the first 5 months as a spin-up period. Model historical and averaged results were output at a daily interval, while the historical fields were output at 00:00 UTC each day.

The physical model setups largely followed an earlier Gulf–COAWST application (Zang et al., 2018, 2019, 2020). Open boundaries were set at the south and east forced by daily water level, horizontal components of 3D current velocity, horizontal components of depth-integrated current velocity, 3D water salinity, and 3D water temperature derived from the Hybrid Coordinate Ocean Model (HYCOM) global analysis products (Bleck and Boudra, 1981; Bleck, 2002) with data assimilated via the Navy Coupled Ocean Data Assimilation system (Cummings, 2005; Cummings and Smedstad, 2013; Fox et al., 2002; Helber et al., 2013). For lateral boundary conditions, we utilized Chapman implicit for the free surface and water level (Chapman, 1985), Flather for depth-integrated momentum (Flather, 1976), gradient for mixing total kinetic energy, and mixed radiation nudging conditions for 3D momentum, temperature, and salinity (Marchesiello et al., 2001). The nudging time steps for the mixed radiation nudging condition were set to 1 d for inflows and 30 d for outflows. The boundary nudging technique was performed at the computational grids along the open boundary. The boundary condition types for passive biological and chemical tracers (i.e., PS, PL, ZS, ZL, ZP, NO₃, NH₄, PON, DON, Si(OH)₄, Opal, PO₄, POP, DOP, and Oxyg) were all prescribed as radiation.

Initial conditions for water level, horizontal components of 3D current velocity, horizontal components of depth-integrated current velocity, 3D water salinity, and 3D water temperature were provided by the same HYCOM products as well. Initial conditions for concentrations of NO₃, PO₄, and Si(OH)₄ were interpolated from measurements provided by the World Ocean Database (WOD, Boyer et al., 2018). Initial conditions for DO concentration were given by the World Ocean Atlas (WOA, Garcia et al., 2019). At the sediment layer, PON_sed, PON_burial, Opal_sed, and Opal_burial were initialized as 0.1 mmol m⁻³. Other biological and chemical tracers were initialized as 0.1 mmol m⁻³ due to the lack of observations.

Atmospheric forcings, including surface wind velocity at 10 m height above sea level, net longwave radiation flux, net shortwave radiation flux, precipitation rate, air temperature 2 m above sea level, sea surface air pressure, and relative humidity 2 m above sea level, were derived from the National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) 6-hourly products (for years prior to 2011, Saha et al., 2010) and NCEP CFS Version 2 (CFSv2) 6-hourly products (for years starting from 2011, Saha et al., 2011) with a horizontal resolution of ∼ 35 and ∼ 22 km, respectively. In our model, 63 rivers were considered to be horizontal point-source forcings along the coastal GoM. They were split into 280 point (red dots in Fig. 2a) sources transporting time-varying salinity (nearly zero), temperature, 3D horizontal momentum (based on the magnitude of river discharges), nutrients (NO₃, NH₄, PO₄, Si(OH)₄, PON, DON, POP, and DOP; Fig. C1), and DO to the computational domain. Locations of river point sources of the Mississippi and the Atchafalaya rivers are shown as red dots in Fig. 2b. For reconstructions of time series of river forcing terms, we composed measurements from various sources, including the US Geological Survey (USGS) National Water Information System (NWIS), National Oceanic and Atmospheric Administration (NOAA) Tides and Currents System (TCS), NOAA National Estuarine Research Reserve System (NERRS), and Mexico National Water Commission (CONAGUA, for rivers in Mexico's territory). Daily averaged river discharges were given based on measurements by USGS NWIS and CONAGUA. The magnitude of river discharges was multiplied by 1.4 to account for adjacent watershed areas and the lateral inflow of tributaries (Warner et al., 2005). River temperature and salinity time series were reconstructed from measurements by USGS NWIS, NOAA TCS, and NOAA NERRS. River nutrient concentrations were provided monthly by USGS NWIS and NOAA NERRS and were extended to daily time series with values in the corresponding months. Riverine DO concentration was set to be a constant (258 mmol m⁻³), assuming that riverine DO was saturated at 25 °C under 1 atm. In addition, tidal forcings were introduced in the hydrodynamic model, taking into account the influences of tidal elevations and tidal currents. There were 13 tidal constituents considered in the model including M2, S2, N2, K2, K1, O1, P1, Q1, MF, MM, M4, MS4, and MN4.

Figure 2 (a) Bathymetry of the entire domain of the Gulf–COAWST, (b) zoomed-in bathymetry plot of the northern Gulf of Mexico (nGoM), and (c) locations of observed inorganic nutrient and DO profiles derived from WOD, SEAMAP, and NOAA's shelf-wide cruises. In (a), locations of river point sources are denoted by red dots. In (b), only bathymetry between 6 and 50 m is mapped with colors; computational meshes are split by solid gray lines; main river channels are denoted by solid blue curves; locations of river point sources of the Mississippi and Atchafalaya rivers are indicated by red dots; sampling locations for SOC and overlying water respiration measurements by McCarthy et al. (2013) are denoted by dark yellow dots.

3 Biogeochemical model validations

3.1 Available measurements

In this section, biogeochemical model validations are conducted for surface inorganic nutrient concentration (i.e., NO₃, PO₄, and Si(OH)₄), types of limited nutrients, ratios of diatoms to total phytoplankton, SOC, DO concentration profiles, spatial distributions of bottom DO concentration, and temporal variability of the hypoxic area against multiple field and lab datasets. Validation of the hydrodynamic model can be found in Zang et al. (2019).

Inorganic nutrient concentrations from WOD and NOAA's shelf-wide cruises were used for model validation. WOD measurements cover the period from 11 January 2007 to 5 July 2009, while the shelf-wide records cover the 2007–2019 period. The types of limited nutrients across the LaTex shelf were discussed based on multiple bioassay studies (Turner and Rabalais, 2013; Quigg et al., 2011; Smith and Hitchcock, 1994; Sylvan et al., 2006, 2007; Zhao and Quigg, 2014; Nelson and Dortch, 1996). The diatom percentage of total phytoplankton was derived from measurements by Chakraborty and Lohrenz (2015) and Schaeffer et al. (2012). The SOC measurements were provided by an incubation study (McCarthy et al., 2013). Available DO concentration profiles were obtained from NOAA-supported midsummer shelf-wide cruises and the Summer Groundfish Survey in GoM supported by the Southeast Area Monitoring and Assessment Program (SEAMAP) conducted annually by the Gulf States Marine Fisheries Commission. The shelf-wide cruises provided 1450 measured profiles with 70 401 available records from 2007 to 2019. There were at least 83 DO profiles for each summer (June–August, except 2016) from the shelf-wide cruise observations. The selected SEAMAP DO dataset covers a time range from 2007 to 2019 with measurements including 1651 profiles with 94 200 sampled records. Locations of the selected profiles from different archives are shown in Fig. 2c. Summer measurements by the shelf-wide cruises were used to validate spatial patterns of bottom DO concentration and time series of summer hypoxic areas. Estimated hypoxic areas by the cruises are available from 2007 to 2020, with a range from 5480 to 22 720 km².

3.2 Surface nutrient concentration

One-to-one comparisons for surface nutrient concentration validation were seldom carried out in previous numerical studies, where spatially averaged or temporally averaged matrices were frequently validated. To provide a more detailed quantification of model performance in surface nutrients, we performed one-to-one differences between simulations and measurements at each sampling location on specific dates. Modeled results showed good agreement with the cruise measurements from both shelf-wide and WOD records (Fig. 3) in terms of magnitudes. A total of 86 % of surface NO₃ differences are within a range of ± 10 mmol m⁻³, with most biases ranging from −2.5 to 0 mmol m⁻³ (56 %, Fig. 3a). It indicates a slight underestimation, which is mostly found in the middle and western shelf (> 150 km from the Mississippi River mouth, Fig. 3b). Surface NO₃ biases exhibit a higher variance near the mouth than in other regions. A total of 92 % of surface PO₄ bias pairs are within ± 1 mmol m⁻³ (Fig. 3c), exhibiting a more even distribution pattern than the NO₃ differences. It results from the model underestimation in the middle and eastern shelf but overestimation in the west (Fig. 3d). A total of 88 % of surface Si(OH)₄ differences are within a range of ± 20 mmol m⁻³, with a slight underestimation (Fig. 3e). We found higher biases near the Mississippi (first to third quartiles within ± 8 mmol m⁻³ at 0–150 km) and the Atchafalaya (−5 to 7 mmol m⁻³ at 150–300 km) River mouths (Fig. 3f) than at the western shelf. Mean Mississippi and Atchafalaya riverine PO₄ concentrations were 2.7 ± 0.7 and 2.3 ± 0.7 mmol m⁻³, respectively, and mean riverine Si(OH)₄ concentrations were 118 ± 23 and 116 ± 21 mmol m⁻³, respectively. Thus, the nutrient concentration bias between simulations and observations is acceptable, considering the possible transient influence from the riverine nutrient loads during a survey.

Figure 3 Comparison of surface nutrient concentration between model hindcasts and cruise measurements (both shelf-wide and WOD) for (a)–(b) NO₃, (c)–(d) PO₄, and (e)–(f) Si(OH)₄. The left bar graphs illustrate the distribution of concentration differences by percentage within specific concentration ranges, while the right box charts show the first quartiles, third quartiles, and medians of the concentration differences against the distance to the Mississippi River mouth.

3.3 Nutrient limitation

Nutrient limitation could vary among different phytoplankton species with different efficiencies in nutrient uptakes. In our model, the Si limitation was modeled only for the PL growth. Depth-averaged nutrient limitation coefficients (see Eqs. A9–A10) over the surface at 1 m depth were compared to bioassay studies. When a modeled coefficient is lower than 0.75, the water body is defined to be limited by the corresponding nutrient for the corresponding phytoplankton group. A bioassay study by Turner and Rabalais (2013) demonstrated that N limitation was more common than P limitation along transects C and F in June and July 2010 (Fig. 4). All July samples were found to be N-limited, while only some June samples along transect C were found to be P-limited with the rest being N-limited. The model mostly captured the dominant N limitation pattern along both transects. As there was a lack of location information in this bioassay study, we could not pinpoint the location of the observed P limitation in Fig. 4. However, our model indicated that the P limitation was more common around the Mississippi River mouth for both phytoplankton groups. In June 2010, transect C, located at the boundary of the modeled N and P limitation, showed that the model could successfully capture the observed spatial pattern of nutrient limitation.

Dominant P limitation adjacent to the Mississippi River mouth was observed in other bioassay studies (e.g., Quigg et al., 2011; Smith and Hitchcock, 1994; Sylvan et al., 2006, 2007) and was also captured by the model, indicated by high percentage occurrences over the simulation period (2007–2020) (Fig. 5b, e). N limitation was mostly found in the shallow parts of the middle and western shelf during spring (Fig. 5a) and became more widespread offshore and eastward in July (Fig. 5d). This pattern was also seen in earlier bioassay estimates (e.g., Quigg et al., 2011; Sylvan et al., 2007; Zhao and Quigg, 2014). The Si limitation occurrence formed a distinct offshore gradient in spring (Fig. 5c). Bioassay studies have illustrated that Si limitation occurred in the eastern shelf during spring (e.g., Quigg et al., 2011; Nelson and Dortch, 1996; Smith and Hitchcock, 1994). The gradient tilted westward in July, indicating a potential oligotrophic water intrusion from deep waters when the circulation pattern changed during the summer months. However, there is a knowledge gap regarding Si limitation over the western shelf region, where no known bioassay studies have been conducted. We gather some clues from the Dortch and Whitledge (1992) study of spring 1988 and summer 1987 in the Mississippi plume (mostly east of 90° W with depth > 50 m), where they found that Si had higher potential as a limiting nutrient than N in summer in high-salinity waters. Salinity in the western shelf is usually high in July due to the changing predominant current system from westward to eastward or northeastward. The low-salinity and Si-rich plume waters can be replaced by deep waters with higher salinity and lower Si. We expect a more Si-limited environment in the western shelf than in other parts during July, which, however, needs further support from additional bioassay studies.

Figure 4 Simulated nutrient limitation coefficient patterns (color) and locations of two sampling transects from the 2010 midsummer shelf-wide cruise (black dots). Modeled nutrient limitation coefficients (for PS, left column; for PL, right column) are averaged over the surface at 1 m depth. A lower coefficient indicates the corresponding nutrient is more limited. Limited nutrient types were examined in the June and July water samples from transects C and F in Turner and Rabalais (2013).

Figure 5 Modeled nutrient limitation occurrences (in percentages) overlaid with locations of observed limited nutrients by bioassay studies in spring (a, b, c) and July (e, f, g). Modeled occurrences are obtained based on the entire simulation period (2007–2020). Solid gray lines indicate bathymetry of 10, 20, and 50 m, while dashed black lines represent the contour lines of 10 %, 50 %, and 70 %.

3.4 Diatom ratios

Cruise observations confirmed that diatoms are one of the dominant phytoplankton groups in the LaTex shelf (Schaeffer et al., 2012; Chakraborty and Lohrenz, 2015). When compared to the Schaeffer et al. (2012) measurements, vertical averages of PS and PL concentration over the surface at 0.5 m depth at the sampled points (black dots in Fig. C2) were extracted from the model hindcast. Statistics of modeled diatom ratios were derived from the daily ratios at the selected locations over the cruise months in 2008. When compared to the Chakraborty and Lohrenz (2015) measurements, we only calculated the modeled diatom ratios at the surface, middle, and bottom layers. Statistics of modeled ratios were given based on the daily ratios at these layers over the cruise regions (polygons shown in Fig. C2) and during cruise months in 2009 and 2010. The modeled ratios reasonably reproduced the measured ones in magnitudes, monthly variability, and cross-shelf variability (Table 1). During the cruise periods in 2008, the range of modeled diatom percentage (59 % to 87 %) matched the measurements well (71 % to 86 %) except for May 2008, when underestimations were found. In 2009, our model results agreed well with the measurements in inner-shelf waters but overestimated the measurements in the mid-shelf regions, especially in the summer and fall of 2009. The measured percentages exhibited salient monthly variations with higher values in winter and spring and lower ones in summer and fall. In the cross-shelf direction, the phytoplankton community shifted from highly diatom-dominated in the inner-shelf waters to less diatom-dominated in the mid-shelf waters, especially in summer. It should be noted that high uncertainty was found in the diatom ratio from both hindcast and measurements (comparable standard deviation against mean values). Therefore, model–measurement biases are expected when comparing statistics derived from a whole month (model hindcast) and a few days (cruise measurements). Then, the biases should be acceptable as the magnitudes of modeled and measured statistics are closed.

Table 1 Comparison of simulated (mean ± 1 SD) and measured (mean ± 1 SD in parentheses) diatom percentage of the total phytoplankton. Note that the statistics for the simulated percentages were calculated based on concentration values over the cruise months and over regions that cover the cruise sampling locations (Fig. C2). The measured percentages by Schaeffer et al. (2012) (for measurements in 2008) were calculated based on biovolume values, while those by Chakraborty and Lohrenz (2015) (for measurements in 2009 and 2010) were given by chlorophyll a attributed to different phytoplankton groups.

3.5 SOC rates

Modeled SOC rates were compared against a laboratory incubation study by McCarthy et al. (2013) at five shelf sites (for the location see the Fig. 1 in that paper) using sediment and water samples collected during six cruises (i.e., July 2008, September 2008, January 2009, August 2009, May 2010, and May 2011). The modeled SOC was averaged over the cruise months for four shelf sites (i.e., F5, C6, B7, and MRM; Fig. 2b). Our model captured the SOC magnitude well. The model generally overestimated the SOC at all sites except for May 2010 at site C6 and August 2009 at site MRM (Fig. 6). The largest overestimations were found in September 2008 when measurements were carried out shortly after hurricanes Gustav and Ike. These measurements tended to provide a low SOC but a high water column respiration, possibly induced by the mixing incurred by storms. Note that the model results shown in Fig. 6 were averaged over an entire month because no exact cruise date information was reported in McCarthy et al. (2013).

Figure 6 Comparison of modeled and measured SOC (unit: µmol m⁻² h⁻¹) at four LaTex shelf sites (dark yellow dots in Fig. 2b). Note that the measurements are provided by the McCarthy et al. (2013) incubation study and the modeled SOC for each sampled site is averaged over the specific months.

3.6 DO profiles

Both the shelf-wide and SEAMAP cruise studies provide high-resolution measurements of DO profiles in the vertical direction, with the observed layers ranging from the surface to bottom. The number of observed layers is close to or even more than that of the modeled layers. Therefore, the observed DO profiles were interpolated to the modeled layers using the nearest interpolation method for the one-to-one comparisons between modeled and observed DO profiles. Mean, median, and 25–75th percentile ranges of the model–observation differences were derived and compared against normalized depths ranging from −1 (bottom) to 0 (surface) (Fig. 7). Most of the biases were within ± 1 mg L⁻¹, indicating a robust model performance in reproducing DO profiles. We noticed that the model tended to overestimate the shelf-wide observed DO by more than 1 mg L⁻¹ but less than 2 mg L⁻¹ on average over the upper layers in shallow waters (Fig. 7a). When validating against the SEAMAP profiles, a wider range of biases was also found at near-surface layers of the shallower water (Fig. 7d) than in deeper waters (Fig. 7e and f). On the one hand, in shallow water, cruise measurements seldom resolved the vertical layers finer than the model where 36 layers were designed, which introduced biases when interpolating the measured profiles to the modeled layers. On the other hand, ROMS tends to overmix the water column in shallow water regardless of the vertical mixing parameterizations chosen (Robertson and Hartlipp, 2017). Despite the slight overestimations of DO profiles, our model results performed better than those of previous numerical studies. For example, DO concentration biases against profile measurements in Yu et al. (2015) were mostly within 2 mg L⁻¹.

Figure 7 Concentration difference statistics of DO profiles between model hindcasts and measurements by (a–c) NOAA's shelf-wide cruises and (d–f) SEAMAP. The statistics are derived from one-to-one differences between hindcasts and measurements for specific sampling locations and dates. The normalized depths of 0 and −1 represent the surface and bottom, respectively. The total counts (n) of profiles within different depth ranges are shown in each panel.

3.7 Spatial distributions of bottom DO and temporal variability of hypoxic area

As the annual NOAA shelf-wide cruises were conducted from the eastern shelf to the west in the summer, the model-simulated bottom DO was resampled following the cruise periods. For example, if the westernmost location of the cruise is 90° W on day 1, the simulated bottom DO concentration over the east of 90° W on that day is extracted. On the following day, if the westernmost location of the cruise is 91° W, the simulation between 91 and 90° W on day 2 is extracted, and so forth. All the extracted frames were blended to reconstruct the spatial distribution of simulated bottom DO concentration during the summer cruise period. Simulated results outside the LaTex shelf and over the deep (> 50 m) and shallow (< 6 m) water regions were excluded since observations were unavailable. Model results showed good agreement with the observations in terms of interannual variability and spatial extent of bottom hypoxic waters (Fig. 8). The spatial distribution of the hypoxic regions varied over different summers. For example, the hypoxic area was small and primarily restricted to nearshore (< 20 m) regions during the summers of 2007, 2009, 2012, 2014, and 2018. The size of the hypoxic zone was more prominent and extended offshore in 2008, 2011, 2013, and 2019. The spatial dispersion of hypoxic waters occurred mostly over the west of the LaTex shelf, where bathymetry gradients were gentle. Over the eastern shelf, the hypoxic water was mostly constrained within a narrow belt. These results suggest that the hypoxia development on the LaTex shelf was complex and generally followed the bathymetry and distances from the major river mouths.

The daily time series of the size of the hypoxic zone was calculated over the LaTex shelf (6–50 m; Fig. 9). There was good agreement between simulated hypoxia zone size and that captured by the shelf-wide cruises in terms of variability and magnitude. The overall correlation coefficient (CC) was 0.69 over the 99 % significant level (Table 2). The 10-year running CCs ranged from 0.66 to 0.76, surpassing at least the 95 % significance threshold. Underestimations were found in 2007, 2008, and 2017 with a root mean squared error (RMSE) of 1693 km², while there are overestimates in other summers of interest with RMSE = 8084 km². The model performed an apparent overestimation for 2019 summer. Nevertheless, biases in other summers were acceptable, considering the relative sporadic converge of cruise data.

Figure 8 Modeled summer bottom DO concentration (colored patches) and NOAA's summer shelf-wide hypoxia observations (black dots and open circles). The black dots and the open circles are indicators of observed bottom hypoxia and normoxia, respectively. The solid gray lines indicate bathymetry of 10, 20, 50, and 100 m.

Figure 9 Comparison of the hypoxic area (in km²) between model simulations and shelf-wide cruise observations from 2007 to 2020 (except 2016). The gray patches denote the cruise periods, while the solid black lines represent the measured hypoxic area.

Table 2 The overall (2007–2020) and 10-year running correlation coefficients (CCs) of summer hypoxic area between model simulations and shelf-wide measurements. Note that the comparison in 2016 is excluded due to the lack of measurement. Superscripts ^* and ${}^{**}$ indicate that the corresponding CCs are above the 95 % and 99 % significant levels, respectively.

4 Results and discussion

4.1 Nutrient limitation

In this study, the riverine nutrient loads from the Mississippi and Atchafalaya rivers were calculated based on measurements from the USGS NWIS. During the investigated period (2007–2020), the riverine N : P ratio was higher than 16:1 during spring and reached its minimum in midsummer to early fall (Fig. 10a). It indicated that P limitation in the shelf could be more severe in spring than in midsummer and early fall (also seeing Fig. 5). Most riverine N : Si ratios fluctuated between 0.5 and 1 and were slightly higher in late spring and summer than in other seasons (Fig. 10b). The riverine N and Si loads were at a similar level when compared to the Redfield ratio of N : Si = 1:1. However, recent studies have pointed out that marine diatoms require a lower N : P : Si ratio ($\mathrm{16}:\mathrm{1}:\mathrm{20}$, Billen and Garnier, 2007; Royer, 2020), indicating that N may be more excessive over Si than previously thought. Riverine Si : P ratios were much higher than 16:1 and 20:1, suggesting that the major river systems transported excessive Si over P to the LaTex shelf. From the perspective of riverine supply, the plume's extent appeared to be more constrained by P availability (see Figs. 4–5) than by N and Si. The limitation effects of N and Si might be relatively similar, given that the N : Si ratio was around 16:20. However, the nutrient limitation is also related to the phytoplankton assimilation efficiency for nutrients (half-saturation coefficients for nutrient uptakes) and the water exchanges between the shelf and the adjacent waters.

Figure 10 Daily time series of ratios of nutrient loads from the Mississippi and Atchafalaya rivers. The dashed black lines denote the nutrient ratios of 16:1, 1:1, and 16:1 in (a), (b), and (c), respectively. The gray patches indicate the late spring and summer (May–August) each year. The capitalized letters M, J, S, and D on the x axis denote the first day of March, June, September, and December, respectively.

The half-saturation coefficient for phytoplankton nutrient uptake is a critical factor associated with nutrient limitation. In our model, PL was parameterized to be more competitive than PS in nutritious waters with a higher half-saturation coefficient. The half-saturation coefficients for NO₃ and NH₄ used in this model study (Table B4) followed the parameterization in Shropshire et al. (2020). The half-saturation coefficients for PO₄ were set as 0.03125 mmol P m⁻³for PS and 0.1875 mmol P m⁻³ for PL, according to the Redfield stoichiometry of N : P = 16:1. This parameterization method was also applied in Laurent et al. (2012) for discussion of P limitation effects in the LaTex shelf. The half-saturation coefficient for Si(OH)₄ ($K_{{\mathrm{SiOH}}_{\mathrm{4}}}$) was set to be 6.0 mmol Si m⁻³, mirroring the choice in Shropshire et al. (2020), although there was no discussion on how this parameter was determined. Uptake kinetic studies for different marine diatom species have suggested a wide range of $K_{{\mathrm{SiOH}}_{\mathrm{4}}}$ from 0.8 to 17.4 mmol Si m⁻³ (Table 3). The average, median, and first and third quartile of the measured coefficients in Table 3 were 5.9, 4.5, 2.3, and 7.0 mmol Si m⁻³, respectively. We opted for the average over the median coefficient in our model, considering the PL group to be a representative marine diatom assemblage. However, the $K_{{\mathrm{SiOH}}_{\mathrm{4}}}$ for a diatom assemblage may shift given changing ambient silicate concentration. For example, as pointed out by Nelson and Dortch (1996), $K_{{\mathrm{SiOH}}_{\mathrm{4}}}$ for the sampled phytoplankton assemblage (dominated by diatom species) remained low from 0.48 to 1.71 mmol Si m⁻³ when the ambient silicate concentration was low between 0.13 and 0.41 mmol Si m⁻³, but it increased to 5.29 mmol Si m⁻³ when the ambient silicate concentration was 4.72 mmol Si m⁻³. Along the Mississippi and Atchafalaya River plumes, which deliver silicate-rich waters to the shelf (average concentrations are 118 ± 23 and 116 ± 21 mmol m⁻³, respectively), the silicate concentration remains high, suggesting a high half-saturation coefficient. We acknowledge that a constant half-saturation coefficient cannot fully capture the dynamics of silicate and diatom outside the plumes, as indicated by Nelson and Dortch (1996). Further investigations and improvements in model parameterization for the dependency of $K_{{\mathrm{SiOH}}_{\mathrm{4}}}$ on silicate concentration are needed in future studies.

Table 3 Half-saturation coefficient (unit: mmol Si m⁻³) for silicate uptake by different diatom species according to multiple uptake kinetic studies.

The changing coastal wind and current systems during summer can lead to significant changes in nutrient distribution, alternating the growth of phytoplankton and summer hypoxia development. Here, we show three snapshots in August 2019 (Fig. 11) when seasonal hypoxia reached its maximum (Fig. 9) to demonstrate the highly varying shelf hydrodynamics and the resultant nutrient dispersion patterns. During spring, the westward alongshore current system dominated the LaTex shelf, while in summer, currents shifted eastward and southward, forming a clockwise circulation in the middle and western shelf (Fig. 11a). This shift not only pushed the river plume eastward but also allowed water intrusion from the west and deep gulf. Waters from the outer shelf were typically high in salinity and low in nutrient content with higher N : Si and lower Si : P ratios than local waters (Fig. 11c–e). Although silicate concentration remained high and was usually excessive in the plume area, the intrusion of deep gulf waters led to an enlarging Si limitation domain in the western LaTex shelf (Fig. 11f–g). The PL concentration and primary production (PS + PL) (Fig. 11h–j) in the western shelf decreased pronouncedly after the intrusion of Si-limited waters. Pronounced declines in PON_sed concentration (Fig. 11k) in the shallow western shelf were also detected 5 d after the primary production decreased. The SOC was expected to decrease, which could relieve the summer bottom hypoxia in the shallow western shelf.

We also note that the upwelling system along the nearshore far western shelf (> 95° W) and the direct transport of PON from the west could affect the evolution of bottom hypoxia on the LaTex shelf. In the Northern Hemisphere, the clockwise circulation system was favorable for the development of coastal upwelling systems, which induced cooling at the surface along the coast (Fig. 11b) and led to elevated concentrations of surface inorganic nitrogen, phosphate, and silicate along the nearshore western shelf. Total surface primary production remained high roughly along the 20 m isobath, where the water column PON concentration was also elevated. The clockwise circulation system carried the PON offshore and northeastward to the LaTex shelf, inducing an increase in the PON_sed pool (around 28° N; Fig. 11k) and SOC. The high alongshore production was limited by N rather than Si or P. However, the N limitation band narrowed around the coastal upwelling zones. Such patterns – including low-Si water intrusion, eastward transport of PON, and a narrow N limitation band in the upwelling zone – were also found in other summer snapshots when the current system changed (e.g., Fig. C3).

Previous bioassay studies suggested the potential Si limitation on the LaTex shelf (Quigg et al., 2011; Nelson and Dortch, 1996; Smith and Hitchcock, 1994; Lohrenz et al., 1999). However, N and P limitations were reported more frequently than Si limitations along the shelf. Part of the reason was that samples collected in previous studies were mainly from the eastern shelf, where N and P typically appeared to be limited. Our understanding of potential nutrient limitations, particularly in the western shelf during the recent decade, still needs to be completed. Nevertheless, this lack of in situ data should not hinder model developments, as indirect evidence supports the potential Si limitation in the western shelf, especially during the summer. For instance, a recent study using in situ incubations and laboratory experiments showed that the oligotrophic open gulf, generally low in N, could also be Si-limited, as indicated by lower maximum growth rates of diatoms compared to other culture and field measurements (Yingling et al., 2022). Additionally, earlier concentration measurements (Dortch and Whitledge, 1992) showed that Si limitation sometimes overwhelmed the N limitation in the deep gulf waters (depth > 50 m). Water exchanges between the LaTex shelf and adjacent deep waters become more pronounced in summer with changes in wind and current systems. The intrusion of low-Si waters can promote the development and expansion of Si limitation, which in turn affects the phytoplankton community and oxygen dynamics. Therefore, the accuracy of the boundary conditions along the LaTex shelf is crucial in biogeochemical modeling. Indeed, earlier numerical studies (e.g., Fennel et al., 2013) emphasized the significance of the correct physical boundary conditions for hypoxia modeling. Our results further illustrate that biogeochemical boundary conditions, such as nutrient concentrations, are as critical as river forcings in influencing the shelf's nutrient distribution, plankton, and oxygen dynamics. These effects have yet to be addressed in numerical studies of the LaTex shelf.

Figure 11 Summer snapshots of (a) sea surface salinity (overlaid with surface current velocity), (b) surface temperature (° C), (c) surface total inorganic nitrogen concentration (mmol N m⁻³), (d) surface phosphate concentration (mmol P m⁻³), (e) surface silicate concentration (mmol Si m⁻³), (f–g) surface nutrient limitation coefficients, (h–i) surface phytoplankton concentration (mmol N m⁻³), and (k) PON_sed concentration (mmol N m⁻³) with a 5 d lag in the nGoM. The nutrient, phytoplankton, and PON_sed concentrations are displayed on the log₁₀ scale.

4.2 Plankton community interactions

On the LaTex shelf (Fig. 2b colored area), total production, primarily supported by the primary production (Fig. 12a), exhibited a bi-peak pattern in spring and summer (e.g., 2007, 2009, 2010, 2014, 2015, 2016, 2017, 2019, and 2020) with both peaks being of similar magnitude. This pattern was hardly captured by numerical models featuring a less complex plankton community (e.g., Fennel et al., 2011) and was seldom reported or discussed even in model simulations where this pattern appeared (see comparisons of modeled and satellite-derived chlorophyll a concentration in Gomez et al., 2018). Satellite-derived chlorophyll a concentration from multiple products, averaged over the LaTex shelf, also showed a bi-peak pattern from March to August (Fig. 12a), closely resembling the pattern observed in our hindcast primary production. A cruise study conducted in March, May, and July 2004 similarly depicted a higher chlorophyll a peak in May and a lower one in July (Quigg et al., 2011). The bi-peak pattern shown was attributed to the negative correlation between PS and PL, where a decrease in PS typically coincided with an increase in PL, and vice versa (Fig. 12b). For example, the peaks in primary production and chlorophyll a, observed from March to May 2019, coincided with the transition from a PS peak to a PL peak. The secondary peak, observed from June to July 2019, was attributed to sustained high PS biomass.

Figure 12 Daily time series of (a) PS + PL and PS + PL + ZS + ZL + ZP biomass (represented by mmol N) as well as (b) PS and PL separated and integrated over the LaTex shelf (Fig. 2b colored area). (a) Monthly time series of regionally averaged (over the LaTex shelf) chlorophyll a concentration (in mg m⁻³) derived from multiple satellite products. The gray patches indicate the late spring and summer (May–August) period of each year. The capitalized letters M, J, S, and D on the x axis denote the first day of March, June, September, and December, respectively.

Competition for nutrients between PS and PL (bottom-up) and grazing pressure from zooplankton (top-down) jointly contribute to the differing fluctuation patterns of PS and PL and the bi-peak total primary production pattern. However, their effects are mostly nonlinear and are not straightforward to explain. We sampled six snapshots around the primary production peaks in the spring (early April) and summer (mid-June) of 2019 to illustrate the responses of both phytoplankton groups to the changing nutrient environments and grazing pressure. Analysis was based on depth averages within the surface at 1 m depth (Figs. 13–14).

In April 2019, a consistent westward current system dominated in the LaTex shelf, corresponding to an east–west elongated river plume region, as indicated by the low-sea-surface-salinity band (Fig. 13a). The spatial pattern of total primary production (PS + PL) followed the plume, within which the PS concentration increased, and PL concentration decreased westward (Fig. 13e–g). These patterns were associated with the nutrient distribution on the shelf (Fig. 13b–d). Inorganic nutrients were abundant around the riverine outlets and diluted and consumed westward following the currents. PL, having a greater half-saturation constant for nutrients than PS, typically achieved higher growth efficiency or reached the maximum growth rate more easily than PS when background nutrients were abundant. By contrast, PS could outcompete PL when nutrient supplies were low. In addition, a downwelling system was established along the shallow coast in the middle and western shelf, leading to decreased nutrient concentrations and allowing PS to outcompete PL. The grazing pressure from the zooplankton group appeared to be minor and did not significantly affect the distribution of PS and PL during these days (Fig. 13h–j).

Pronounced bottom-up and top-down effects on the primary production were found around the biomass peak in June 2019, coinciding with a shift in the coastal current system to a northward direction (Fig. 14a). The northward currents not only constrained the river plume but also introduced oligotrophic deep water, as evidenced by the high surface salinity, to the inner shelf. Note that the discharges of the Mississippi and Atchafalaya rivers remained high from May to July 2019 (Fig. C1). A distinct difference in the patterns of PS and PL was observed between 89 and 93° W and between 93 and 97° W (Fig. 14f–g). In the former region, where constrained river plumes and oligotrophic water intrusions were detected, PS exhibited a higher nutrient uptake efficiency than PL. In contrast, PL concentration was slightly higher than PS concentration in the latter regions, where the plume was pushed offshore. However, two areas of low PS concentration and corresponding high PL concentration were identified between 93 and 96° W, nearshore stretching from southwest to northeast, and between 91 and 92° W, stretching from nearshore to offshore. In these regions, the concentration of ZS, which grazes on PS only (Fig. 14h), was high, exerting strong grazing pressure on PS but inversely allowing PL to bloom (Fig. 14h).

The results indicated that the responses of PS, PL, and PS + PL to the riverine nutrient loads were nonlinear due to the mixing among the waters on the shelf, from the river, and intruding from the deep ocean. The riverine nutrient supplies were much greater in June 2019 than in March–April 2019 (Fig. C1). Higher primary production and PL concentration in June would have been expected if a nutrient-based linear relationship had been applied. However, as shown in the model and the satellite products, primary production was higher in April than in June. This indicates that variations of phytoplankton concentration are affected not only by riverine nutrient inputs but also the current system, which limits the expansion of river plumes, pronounced upwelling or downwelling, and water exchanges with the oligotrophic open ocean. In the April and June 2019 snapshots, mesoscale eddies were found south of the Mississippi River outlets. The intensity and impact area of the June eddy were greater than those of the April eddy, causing a more pronounced northward flow and more constrained river plumes along the shelf in June. These eddy systems are known as loop current eddy (LCE) systems, which can propagate westward and interact with the LaTex shelf waters after the detachment from the GoM Loop Current (LC). A recent study indicated that LCE has distinct bio-optical properties (e.g., temperature, salinity, density, DO concentration, and chlorophyll a concentration) from the surrounding waters, highlighting the importance of open-ocean dynamics to the shelf biogeochemical processes (Zhang et al., 2023). Another recent study analyzing water samples from the LaTex shelf emphasized the significant impact of mesoscale circulation features on the summer planktonic community composition (Anglès et al., 2019). This study revealed that between 20 and 25 June 2013, diatoms proliferated on the western shelf, where upwelling was detected, whereas the flagellate group dominated within the river plumes. From 18 to 23 June 2014, diatoms and flagellates bloomed in proximity to the Mississippi River and Atchafalaya River outlets, respectively. In contrast, blooms on the western shelf were characterized by a mixture of the two phytoplankton groups. Similar patterns were observed in our model results, as depicted in Figs. C4–C5.

In addition to the impacts of upwelling and LCE systems, direct advection of river outflow waters by coastal currents was also found to be significant for phytoplankton community composition, carbon export, and the associated bottom DO conditions based on other field studies in the nGoM (Chakraborty and Lohrenz, 2015) and northeastern GoM (Qian et al., 2003). Our results suggest that the grazing pressure exerted by zooplankton groups can be variable, manifesting as significant in some instances while remaining minimal in others. Laboratory experiments on surface water samples collected around the Mississippi River outlets in May 1993 suggested significant grazing pressures by microzooplankton on phytoplankton growth (Strom and Strom, 1996). However, no salient grazer impact was found on phytoplankton growth according to bioassay studies on the water samples collected around the plumes in April and August 2012 (Zhao and Quigg, 2014). Other unmodeled factors can also affect shelf primary production. For example, a reduction of chlorophyll a between 2011 and 2014 detected in the nGoM was attributed to the Deepwater Horizon oil spill disaster in 2010 (Li et al., 2019). Incorporating a complex community into the model to address the nonlinear interactions among different plankton groups enhances our understanding of the primary production variability and associated DO dynamics on the LaTex shelf (e.g., the bi-peak patterns that were seldom discussed before).

Figure 13 Snapshots of (a) sea surface salinity (overlaid with surface current velocity), (b) surface total inorganic nitrogen concentration (mmol N m⁻³), (c) surface phosphate concentration (mmol P m⁻³), (d) surface silicate concentration (mmol Si m⁻³), (e–g) surface phytoplankton concentration (mmol N m⁻³), and (h–j) surface zooplankton concentration (mmol N m⁻³). The nutrient and plankton concentrations are displayed on the log₁₀ scale.

Figure 14 Same as Fig. 13, but for snapshots from 14 to 16 June 2019.

4.3 A re-examination of LaTex shelf DO dynamics

In this section, we specify the bottom waters as the layers within 2 m above the sea floor, while the upper waters represent all layers above this 2 m bottom layer. The purpose is to understand the contributions of different processes, including water column biochemistry, air–sea flux (in upper layers), SOC (in bottom layers), and water transports (advection + diffusion), to the daily variations of DO in the LaTex shelf during summers (May–August) of 2007–2020.

In the upper LaTex shelf, daily DO changes were primarily driven by shelf physics and local water column biochemistry (Fig. 15a), as reflected by their significant contributions to the variability and magnitude. The advection and diffusion terms together explained the greatest spatiotemporal variability of total DO changes. The ranges of the first and the third quartiles were closely shown in the total rate of changes (−124 to 107 mmol O₂ m⁻² d⁻¹) and changes by water transports (−117 to 72 mmol O₂ m⁻² d⁻¹). Detailed separation of the water transport terms indicated that horizontal advection of DO contributed the most to the variability of the physical terms. The water column biochemistry contributed the second most to total DO variability, with a wide range of first and third quartiles (−41 to 96 mmol O₂ m⁻² d⁻¹). The phytoplankton groups contributed positively to the upper DO pool, with the majority of the contribution from the PS group. PS biomass was usually higher than PL biomass in summer when the allocation of nutrients was more favorable for the growth of PS. The net DO changes by water column biochemistry could be negative, indicating net metabolism, which was also reported by previous field studies demonstrating consistent net water column heterotrophy across the Louisiana shelf (e.g., Murrell et al., 2013). The air–sea interactions contributed negatively to the total DO changes and accounted for the least contribution. This indicated that the upper LaTex shelf was mostly a source of oxygen to the atmosphere during summer.

In the bottom layers, the DO variability was controlled by SOC and water transports (Fig. 15b). SOC was a steady (narrow range of quartiles) but major DO loss term (median −32 mmol O₂ m⁻² d⁻¹, first quartile −45 mmol O₂ m⁻² d⁻¹, and third quartile −24 mmol O₂ m⁻² d⁻¹), driving the total rate of changes in DO to be negative at most shelf grids during summer (median −8 mmol O₂ m⁻² d⁻¹, first quartile −32 mmol O₂ m⁻² d⁻¹, and third quartile 11 mmol O₂ m⁻² d⁻¹). The advection and diffusion terms together acted as a major source of DO in the bottom layers (median 21 mmol O₂ m⁻² d⁻¹, first quartile 7 mmol O₂ m⁻² d⁻¹, and third quartile 48 mmol O₂ m⁻² d⁻¹). However, they hardly offset the DO loss due to SOC. Such a positive contribution to DO by physical transports was mainly a result of steady and strong net DO supplies through vertical diffusion, as the variability and magnitude of DO changes due to total advection were less pronounced than those due to vertical diffusion. The vertical diffusion of DO is influenced by both water stratification and vertical DO concentration gradient. Water stratification results from multiple processes, including river plume dynamics, tidal dynamics, wind patterns, and surface heating and cooling, and has been identified as an important indicator of bottom DO supply (Hetland and DiMarco, 2008; Bianchi et al., 2010; Fennel et al., 2011, 2013, 2016; Justić and Wang, 2014; Wang and Justić, 2009; Feng et al., 2014; Yu et al., 2015; Laurent et al., 2018). The variation of the vertical gradient was more related to the DO dynamics in the upper layers than in the bottom, as the DO variability is more pronounced in the upper layers (wider range in total rate of changes). Thus, while SOC and water stratification play crucial roles in DO changes in the bottom layers, DO changes in the upper shelf can affect the bottom DO through vertical diffusion.

Figure 15 Depth-integrated rate of changes in DO due to different modeled processes in (a) the upper layers and (b) the bottom layers. The total rate of changes is the summation of DO sources and sinks by three groups of contributors (water column biochemistry, DO transports, and air–sea flux in upper layers or SOC in bottom layers) separated by vertical gray lines. In each group, DO changes by specific processes are illustrated by black boxes. Boxes represent the first and third quartiles, with lower and upper whiskers extending to the lowest and highest values within the 1.5 interquartile range of the first and third quartiles, respectively. The median is indicated by a black line in the middle of the boxes. Statistics are summarized from the summer (May–August) records of 2007–2020 for all grid cells in the LaTex shelf.

The interactions within the plankton community (e.g., competition for nutrients and grazing pressure), which led to biomass differences, also resulted in different DO patterns at the bottom layer. Such impacts became more apparent when the DO contribution by water biochemistry outweighed that from transport processes in the upper ocean. For illustration, three summer snapshots of 14–16 June 2019 (Figs. 16 and 17) were sampled when widespread bottom hypoxia was detected. The water column biochemical processes contributed more than 50 % of total DO changes in most computational cells in the upper layers (Fig. 16a). First of all, the DO contribution by phytoplankton, zooplankton, and microbes exhibited distinct spatiotemporal patterns, complicating the net DO changes in the upper layers. Generally, the PS and PL groups enhanced DO levels, whereas zooplankton and microbes tended to deplete DO. During 14 June 2019, the DO losses by biochemical processes (Fig. 16b) in the shallow western shelf were mostly attributed to high ZS metabolism (Fig. 16e). The net DO gains between 91.5 and 92.5° W reflected high PL concentrations (Fig. 14g) and the associated high DO supplies (Fig. 16d). The scattered DO losses over the shelf were primarily due to the homogenously high DO consumption by microbes (Fig. 16h). During 15 and 16 June 2019, when DO supplies by PS and PL (Fig. 16c–d) increased, net DO gains predominated in the shelf (Fig. 16b). However, the net DO gains in the western (> 92.5° W) and eastern (< 91.5° W) shelf were mainly contributed by PS, while those in the middle shelf were by PL.

At the same time, changes in upper DO could affect the bottom DO through vertical diffusion, the spatial patterns (mostly positive; Fig. 17b) and daily variability of which aligned with biochemical DO alterations in the upper layers (Fig. 16b). However, water column stratification, as indicated by the potential energy anomaly (PEA; Fig. 17a), resulted in noticeable spatial disparities in the vertical diffusion of DO. On 15 June 2019, for example, the effects of vertical diffusion were weakened in areas that featured strong stratification, as evidenced by high PEA values. In contrast, in regions of weak stratification, such as the shallow waters between 90.5 and 92.5° W, vertical diffusion was markedly stronger. During the sampled period, among various factors (i.e., total advection, horizontal diffusion, water column biochemistry, and SOC), the vertical diffusion term contributed the most to the total rate of changes in bottom DO, especially over the middle shallow shelf. As the rates of changes were daily-averaged and the bottom DO concentration was sampled at 00:00 UTC on each sampled day (Fig. 17i), the elevated bottom DO level and relief of bottom hypoxia in the shallow middle shelf on 16 June 2019 were mainly due to the significant vertical diffusion on the preceding day, driven by high PL-supported DO sources and weak water stratification. Thus, through the interactions within the community in the upper ocean and DO diffusion processes between the upper and bottom layers, the influence of planktonic community complexity on the bottom DO dynamics and the hypoxia evolution is evident.

The influence of SOC and water stratification on bottom hypoxia in the LaTex shelf has been well-documented. Yet, the role of planktonic community complexity has received scant attention in prior numerical and observational studies. This study devoted considerable effort to validating various factors from nutrient dynamics (concentration and limitation types) to phytoplankton composition (diatom ratio and temporal variations in total primary production) and oxygen variables (SOC, DO profiles, and hypoxia patterns). Our findings illustrated how both bottom-up mechanisms (phytoplankton competition for nutrients) and top-down effects (zooplankton grazing on phytoplankton) shape plankton composition, thereby influencing DO levels in the upper water column and affecting subsequent changes in bottom DO and hypoxia patterns through physical transports (e.g., vertical diffusion). The insights obtained suggest that the impacts of planktonic community complexity on bottom DO and hypoxia patterns could be of high importance.

Nonetheless, incorporating a more complex plankton community in the model requires reasonable parameterizations for different groups to represent their interactions. The large number of parameters can sometimes hamper the reliability of a biogeochemical model due to the lack of support from in situ observations or laboratory experiments. This is also a critical reason why prevailing lower-trophic biogeochemical models are often “oversimplified”. Even in complex models, the number of plankton functional groups considered needs to be constrained to avoid over-parameterization. For example, there are two phytoplankton and two zooplankton functional groups in the PISCES (Aumont and Bopp, 2006) and CoSiNE models (Chai et al., 2002), three phytoplankton and two zooplankton functional groups in the PlankTOM5 model (Buitenhuis et al., 2010), and three phytoplankton groups and one zooplankton functional group in the CCSM-BEC model (Moore et al., 2004).

Figure 16 Snapshots of DO contribution by the (a) water column biochemical processes (percentages) in the upper layers, DO gain and loss rates (mmol m⁻² d⁻¹) due to (b) water column biochemical processes, (c) PS, (d) PL, (e) ZS, (f) ZL, (g) ZP, and (h) microbes in the upper layers. The percentage contribution is related to the sum of absolute DO changes due to water column biochemical processes, water transports (advections and diffusions), and air–sea fluxes in the upper layers. The solid black lines in (a) indicate the −50 % and 50 % contour lines.

Figure 17 Snapshots of (a) potential energy anomaly (PEA; J m⁻³), DO gain and loss rates (mmol m⁻² d⁻¹) due to (b) vertical diffusion (vdiff), (c) horizontal advection (hadv), (d) vertical advection (vadv), (e) horizontal diffusion (hdiff), (f) water column biochemical processes in the bottom layers, (g) SOC, (h) total bottom DO gain and loss rates (mmol m⁻² d⁻¹), and (i) bottom DO concentration (mg L⁻¹). Rate snapshots are daily averages, while snapshots of state variables (i.e., PEA and bottom DO concentration) are extracted at 00:00 UTC on each sampled day.

5 Conclusions

In this study, we modified a three-dimensional coupled hydrodynamic–biogeochemical model (NEMURO) and adapted it to the GoM to investigate the mechanisms of bottom DO variability in the LaTex Shelf from 2007 to 2020. In addition to N and Si, a P flow was embedded into NEMURO to account for the impacts of P limitation on phytoplankton growth rates. Drawing upon the SOC scheme of the instantaneous remineralization developed by Fennel et al. (2006), a pool of sedimentary PON was added to capture temporal delays in SOC relative to the peak of plankton blooms. The model reproduced the surface inorganic nutrient concentration (i.e., nitrate, phosphate, and silicate), nutrient limitation patterns, the ratio of diatom to total phytoplankton, and the magnitude of SOC well. The model's robustness in DO simulation was affirmed via comparison of the DO profiles against cruise observations from two different databases, comparison of spatial distributions of bottom DO, and time series of the hypoxic area against shelf-wide cruise observations.

Model results revealed that the changing dominant current system in summer can significantly alter the distribution of shelf nutrients and types of nutrient limitations. While N limitation and P limitation dominate the Mississippi and Atchafalaya River plume area, Si limitation becomes pronounced as the coastal current system shifts from westward to eastward or northward, facilitating the intrusion of low-Si waters from the west and the deep gulf. This effect, particularly evident on the western shelf, has rarely been addressed in previous studies on nutrient limitation. Model results also indicated that under a westward background current system, upwellings can enhance nearshore surface nutrient content, with the two modeled phytoplankton functional groups, PS and PL, exhibiting distinct responses to the redistribution of surface nutrients.

Our findings underscore the importance of incorporating complex community dynamics and sophisticated nonlinear interactions into biogeochemical models to capture the variability in primary production on the LaTex Shelf. The model identified a bi-peak production pattern in spring and early summer, aligning with satellite-derived chlorophyll a variations – a pattern not commonly reported in earlier research. We linked this bi-peak pattern to plankton community interactions, including both bottom-up and top-down effects, as demonstrated in the sampled spring and summer snapshots. Changes in nutrient distribution arising from interactions between the LaTex shelf and its adjacent waters, the passages of LCE, the formation of upwelling or downwelling systems, and variations in river plume patterns are crucial in influencing plankton interactions, highlighting the important role of open-ocean dynamics and boundary conditions along the LaTex shelf in LaTex biogeochemical modeling.

While the effects of SOC and water stratification on bottom hypoxia are well-recognized, our study illuminates how plankton composition, influenced by bottom-up and top-down effects, can affect DO levels in the upper water column and lead to changes in bottom DO and hypoxia patterns through physical transport processes, such as vertical diffusion. These insights suggest the potential impacts of planktonic community complexity on bottom DO and hypoxia patterns, emphasizing the need for future in situ and modeling efforts.

Appendix A: Expressions of process terms modified in this study

Detailed descriptions of related terms and parameters are listed in Appendix B.

A1 Update gross primary production of PS and PL due to the additional phosphate limitation

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathrm{GppPSn}=\mathrm{GppNPS}+\mathrm{GppAPS},\text{(A2)}& {\displaystyle }\mathrm{GppPLn}=\mathrm{GppNPL}+\mathrm{GppAPL},\end{array}$$

where

$$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}\mathrm{GppNPS}=& \mathrm{PSn}{V}_{\mathrm{maxS}}\exp \left(K_{\mathrm{GppS}}\mathrm{TMP}\right)\\ & \left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{PS}}}{V_{\mathrm{maxS}}}}{I}_{\mathrm{PS}}\right)\right]\\ & \exp \left(-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{PS}}}{V_{\mathrm{maxS}}}}{I}_{\mathrm{PS}}\right)\mathrm{NutlimPS}\mathrm{RnewS},\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}\mathrm{GppAPS}=& \mathrm{PSn}{V}_{\mathrm{maxS}}\exp \left(K_{\mathrm{GppS}}\mathrm{TMP}\right)\\ & \left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{PS}}}{V_{\mathrm{maxS}}}}{I}_{\mathrm{PS}}\right)\right]\\ & \exp \left(-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{PS}}}{V_{\mathrm{maxS}}}}{I}_{\mathrm{PS}}\right)\mathrm{NutlimPS}\\ & (\mathrm{1}-\mathrm{RnewS}),\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}\mathrm{GppNPL}=& \mathrm{PLn}{V}_{\mathrm{maxL}}\exp \left(K_{\mathrm{GppL}}\mathrm{TMP}\right)\\ & \left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{PL}}}{V_{\mathrm{maxL}}}}{I}_{\mathrm{PL}}\right)\right]\\ & \exp \left(-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{PL}}}{V_{\mathrm{maxL}}}}{I}_{\mathrm{PL}}\right)\mathrm{NutlimPL}\mathrm{RnewL},\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}\mathrm{GppAPL}=& \mathrm{PLn}{V}_{\mathrm{maxL}}\exp \left(K_{\mathrm{GppL}}\mathrm{TMP}\right)\\ & \left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{PL}}}{V_{\mathrm{maxL}}}}{I}_{\mathrm{PL}}\right)\right]\\ & \exp \left(-{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{PL}}}{V_{\mathrm{maxL}}}}{I}_{\mathrm{PL}}\right)\mathrm{NutlimPL}\\ & \left(\mathrm{1}-\mathrm{RnewL}\right),\end{array}\text{(A7)}& {\displaystyle }\begin{array}{rl}\mathrm{RnewS}=& {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{S}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{S}}}\right)}}\\ & {\displaystyle \frac{\mathrm{1}}{\frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{S}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{S}}}\right)}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{{\mathrm{NH}}_{\mathrm{4}}+K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{S}}}}},\end{array}\text{(A8)}& {\displaystyle }\begin{array}{rl}\mathrm{RnewL}=& {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{L}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{L}}}\right)}}\\ & {\displaystyle \frac{\mathrm{1}}{\frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{L}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{L}}}\right)}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{{\mathrm{NH}}_{\mathrm{4}}+K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{L}}}}},\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}\mathrm{NutlimPS}=& min({\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{S}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{S}}}\right)}}\\ & +{\displaystyle \frac{{\mathrm{NH}}_{\mathrm{4}}}{{\mathrm{NH}}_{\mathrm{4}}+K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{S}}}},{\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K_{{\mathrm{PO}}_{\mathrm{4}}\mathrm{S}}}}),\end{array}\text{(A10)}& {\displaystyle }\begin{array}{rl}\mathrm{NutlimPL}=& min({\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{\left({\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}\mathrm{L}}\right)\left(\mathrm{1}+\frac{{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{L}}}\right)}}\\ & +{\displaystyle \frac{{\mathrm{NH}}_{\mathrm{4}}}{{\mathrm{NH}}_{\mathrm{4}}+K_{{\mathrm{NH}}_{\mathrm{4}}\mathrm{L}}}},{\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K_{{\mathrm{PO}}_{\mathrm{4}}\mathrm{L}}}},\\ & {\displaystyle \frac{{\mathrm{SiOH}}_{\mathrm{4}}}{{\mathrm{SiOH}}_{\mathrm{4}}+K_{{\mathrm{SiOH}}_{\mathrm{4}}\mathrm{L}}}}),\end{array}\text{(A11)}& {\displaystyle }\begin{array}{rl}{I}_{\mathrm{PS}}& =\mathrm{PAR}\mathrm{frac}\exp \mathit{\{}z\mathrm{AttSW}+\mathrm{AttPS}{\int }_z^{\mathrm{0}}[\mathrm{PSn}\left(\mathit{\zeta }\right)\\ & +\mathrm{PLn}\left(\mathit{\zeta }\right)\left]\mathrm{d}\mathit{\zeta }\mathit{\right\}},\end{array}\text{(A12)}& {\displaystyle }\begin{array}{rl}{I}_{\mathrm{PL}}& =\mathrm{PAR}\mathrm{frac}\exp \mathit{\{}z\mathrm{AttSW}+\mathrm{AttPL}{\int }_z^{\mathrm{0}}[\mathrm{PSn}\left(\mathit{\zeta }\right)\\ & +\mathrm{PLn}\left(\mathit{\zeta }\right)\left]\mathrm{d}\mathit{\zeta }\mathit{\right\}}.\end{array}\end{array}$$

A2 Update aerobic decomposition from PON to NH₄ and from DON to NH₄ due to the introduction of oxygen dependency

$$\begin{array}{}\text{(A13)}& {\displaystyle }\mathrm{DecP}\mathrm{2}\mathrm{N}=\mathrm{PON}\mathrm{VP}\mathrm{2}{\mathrm{N}}_{\mathrm{0}}\exp \left(K_{\mathrm{P}\mathrm{2}\mathrm{N}}\mathrm{TMP}\right)s\widehat{r},\text{(A14)}& {\displaystyle }\mathrm{DecD}\mathrm{2}\mathrm{N}=\mathrm{PON}\mathrm{VD}\mathrm{2}{\mathrm{N}}_{\mathrm{0}}\exp \left(K_{\mathrm{D}\mathrm{2}\mathrm{N}}\mathrm{TMP}\right)\widehat{r},\end{array}$$

where

$$\begin{array}{}\text{(A15)}& \widehat{r}=max\left[{\displaystyle \frac{max\left(\mathrm{0},\mathrm{Oxyg}-{\mathrm{Oxyg}}_{\mathrm{th}}\right)}{K_{\mathrm{Oxyg}}+\mathrm{Oxyg}-{\mathrm{Oxyg}}_{\mathrm{th}}}},\mathrm{0}\right].\end{array}$$

A3 Update water column nitrification due to the introduction of oxygen dependency and light limitation

$$\begin{array}{}\text{(A16)}& \mathrm{Nit}={\mathrm{Nit}}_{\mathrm{0}}\exp \left(K_{\mathrm{Nit}}\mathrm{TMP}\right)\mathrm{LgtlimN}\widehat{r},\end{array}$$

where

$$\begin{array}{}\text{(A17)}& {\displaystyle }\mathrm{LgtlimN}=\mathrm{1}-max\left(\mathrm{0},{\displaystyle \frac{I_{\mathrm{N}}-I_{\mathrm{0}}}{I_{\mathrm{N}}-I_{\mathrm{0}}+k_{\mathrm{I}}}}\right),\text{(A18)}& {\displaystyle }\begin{array}{rl}{I}_{\mathrm{N}}=& \mathrm{PAR}\mathrm{frac}\exp \mathit{\{}z\mathrm{AttSW}+max\left(\mathrm{AttPS},\mathrm{AttPL}\right)\\ & {\int }_z^{\mathrm{0}}\left[\mathrm{PSn}\left(\mathit{\zeta }\right)+\mathrm{PLn}\left(\mathit{\zeta }\right)\right]\mathrm{d}\mathit{\zeta }\mathit{\}}.\end{array}\end{array}$$

A4 Additional SOC term

$$\begin{array}{}\text{(A19)}& \mathrm{SOC}=\mathrm{8.3865}{\mathrm{PON}}_{\mathrm{sed}}\mathrm{VP}\mathrm{2}{\mathrm{N}}_{\mathrm{0}}\exp \left(K_{\mathrm{P}\mathrm{2}\mathrm{N}}\mathrm{TMP}\right)\end{array}$$

Appendix B: Descriptions of terms and parameters

Table B1 Descriptions of state variables.

Table B2 Descriptions of related terms involved in the phosphorus cycle and nutrient limitation. Superscripts ^* and + denote the fact that the mathematic expressions of corresponding terms are the same as those in Kishi et al. (2007) and Shropshire et al. (2020), respectively. Expressions of terms with no superscript are updated and reported in Appendix A.

Table B3 Descriptions of other variables.

Table B4 Descriptions and values of all model parameters. Superscripts S, L, F06, and F13 denote the fact that the corresponding parameters follow Shropshire et al. (2020), Laurent et al. (2012), Fennel et al. (2006), and Fennel et al. (2013), respectively. Superscript ^* indicates that the corresponding parameters are from this study.

Appendix C: Supporting figures

Figure C1 Daily time series (2007–2020) of river discharges of freshwater, nitrate, phosphate, and silicate from the Mississippi and Atchafalaya rivers.

Figure C2 The model computational meshes over which the regionally averaged diatom ratios are calculated for validation purposes. Black dots indicate the sampling locations in Schaeffer et al. (2012), while the regions restricted by two black polygons are two regions (i.e., inner shelf and mid-shelf) where samples were collected in the Chakraborty and Lohrenz (2015) study.

Figure C3 Summer snapshots of (a) sea surface salinity (overlaid with surface current velocity), (b) surface temperature (°C), (c) surface total inorganic nitrogen concentration (mmol N m⁻³), (d) surface phosphate concentration (mmol P m⁻³), (e) surface silicate concentration (mmol Si m⁻³), (f–g) surface nutrient limitation coefficients, (h–i) surface phytoplankton concentration (mmol N m⁻³), and (k) PON_sed concentration (mmol N m⁻³) with a 5 d lag in the nGoM. The nutrient, phytoplankton, and PON_sed concentrations are displayed on the log₁₀ scale.

Figure C4 Snapshots of (a) sea surface salinity (overlaid with surface current velocity), (b) surface total inorganic nitrogen concentration (mmol N m⁻³), (c) surface phosphate concentration (mmol P m⁻³), (d) surface silicate concentration (mmol Si m⁻³), (e–g) surface phytoplankton concentration (mmol N m⁻³), and (h–j) surface zooplankton concentration (mmol N m⁻³). The nutrient and plankton concentrations are displayed on the log₁₀ scale.

Figure C5 Same as Fig. C4, but for snapshots from 21 to 23 June 2014.