Biogeosciences Simulating the effects of phosphorus limitation in the Mississippi and Atchafalaya River plumes

The continental shelf of the northern Gulf of Mexico receives high dissolved inorganic nitrogen and phosphorus loads from the Mississippi and Atchafalaya rivers. The nutrient load results in high primary production in the river plumes and contributes to the development of hypoxia on the Louisiana shelf in summer. While phytoplankton growth is considered to be typically nitrogen-limited in marine waters, phosphorus limitation has been observed in this region during periods of peak river discharge in spring and early summer. Here we investigate the presence, spatio-temporal distribution and implications of phosphorus limitation in the plume region using a circulation model of the northern Gulf of Mexico coupled to a multi-nutrient ecosystem model. Results from a 7-yr simulation (2001–2007) compare well with several sources of observations and suggest that phosphorus limitation develops every year between the Mississippi and Atchafalaya deltas. Model simulations show that phosphorus limitation results in a delay and westward shift of a fraction of river-stimulated primary production. The consequence is a reduced flux of particulate organic matter to the sediment near the Mississippi delta, but slightly enhanced fluxes west of Atchafalaya Bay. Simulations with altered river phosphate concentrations ( ±50 %) show that significant variation in the spatial extent of phosphorus limitation ( ±40 % in July) results from changes in phosphate load.


Introduction
The Mississippi-Atchafalaya river basin is the world's thirdlargest river basin, over 3.2 × 10 6 km 2 in size and constitutes nearly 80 % of the United States freshwater discharge into the Gulf of Mexico (Dunn, 1996).The Mississippi and Atchafalaya rivers, which contribute 2/3 and 1/3 of the discharge respectively, are the main sources of freshwater to the Louisiana shelf and represent about 95 % of its total nitrogen (N) and phosphorus (P) load (Dunn, 1996).Changes in agricultural practices in the Mississippi-Atchafalaya river basin -in particular an increase in the use of fertilizers -have tripled dissolved inorganic N concentration (DIN, mainly nitrate + nitrite) since the 1960s in the lower Mississippi River (Turner and Rabalais, 1991;Goolsby et al., 2001).This increase in nutrient load has resulted in eutrophication, with high chlorophyll concentration and high primary production (Turner and Rabalais, 1994;Lohrenz et al., 1997), and contributes to the development of hypoxic bottom waters on the Louisiana shelf in summer (Rabalais et al., 2002;Greene et al., 2009).Total P (TP) in the river has also increased, but dissolved inorganic P (DIP) has decreased since the 1980s (Lohrenz et al., 2008) modifying the stoichiometric balance of DIN and DIP inputs.
Primary production on the Louisiana shelf is usually Nlimited during periods of low river flow in late summer and fall (Rabalais et al., 2002), whereas light limitation occurs in late fall and winter and near the river delta (Lohrenz et al., 1999;Fennel et al., 2011).The DIN to DIP ratio (DIN : DIP) in the lower Mississippi River varies seasonally, reaching values well above the Redfield N : P ratio of 16 : 1 during the period of peak discharge in spring and early summer when DIN > 100 mmol m −3 .
These high N : P ratios should lead to P limitation of primary production assuming that nutrients are removed near Redfield stoichiometry.DIP is considered limiting when its concentrations reach values of less than 0.2 mmol m −3 (Dortch and Whitledge, 1992).It has been suggested that TP rather than DIP should be used in the calculation of N : P (Rabalais et al., 2002;Dodds, 2006;Turner et al., 2007), which would indicate near-Redfield stoichiometry in the lower Mississippi River.However, while dissolved organic P (DOP) may be an important source of P on the Louisiana shelf (Dagg et al., 2007), DIP is the only form of P that is readily available to phytoplankton.Moreover, the DIN : DIP observations are consistent with physiological measurements and nutrient addition bioassays carried out near the Mississippi delta region (Dortch and Whitledge, 1992;Smith and Hitchcock, 1994) and across the Louisiana shelf (Sylvan et al., 2006(Sylvan et al., , 2007(Sylvan et al., , 2011;;Quigg et al., 2011) all indicating that P limitation occurs during spring and early summer.The period of P limitation coincides with the peak in primary production (Lohrenz et al., 1997) that is thought to contribute to the development of summer hypoxia (Rabalais et al., 2002).Phosphorus limitation may therefore have a significant role in the biogeochemistry of the Louisiana shelf, leading to a delay in the assimilation of riverine DIN by phytoplankton and spreading the river-induced enhancement of primary production over a larger shelf area (Quigg et al., 2011).This would modify the distribution and magnitude of particulate organic matter (POM) fluxes to the sediment and could subsequently affect the development and location of summer hypoxia.Clearly, the effects of P limitation on the biogeochemistry of the Louisiana shelf and the consequences of changing DIN : DIP stoichiometry in the Mississippi and Atchafalaya rivers need to be investigated.
Recently, an N-based model of the lower trophic ecosystem has been coupled to a realistic 3-dimensional circulation model of the northern Gulf of Mexico (Fennel et al., 2011).
The model successfully simulates the seasonal cycle and the spatial distribution of nitrate and phytoplankton across the Louisiana shelf.A 15-yr simulation has provided insights into the temporal and spatial patterns of primary production and phytoplankton loss terms in the region and was used to explain the relationship between primary production and N loads (Fennel et al., 2011).The model does not include DIP and therefore assumes that primary production is limited by light and N only.Here, we extend the model of Fennel et al. (2011) to explicitly simulate the dynamics of DIP.A 7-yr simulation (2001)(2002)(2003)(2004)(2005)(2006)(2007) is analyzed to assess if and how P limitation affects patterns of primary production across the Louisiana shelf.The results compare well with spatial and temporal patterns of P and N limitation observed during the nutrient limitation studies of Sylvan et al. (2006Sylvan et al. ( , 2007Sylvan et al. ( , 2011) ) and Quigg et al. (2011) and corroborate their hypothesis of delayed primary production in the Mississippi delta resulting from P limitation.This effect is analyzed for several zones across the shelf, and nutrient dynamics in the Mississippi and Atchafalaya River plumes are compared.Finally, sensitivity experiments are presented to assess the effect of altered DIP loads on the extent of P limitation.

Circulation model
The Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008) was configured to simulate water circulation on the northern Gulf of Mexico shelf near the Mississippi and Atchafalaya river outflow region.The model grid covers the continental shelf region from 94.6 • W to 87.8 • W with a variable spatial resolution (Fig. 1).The horizontal resolution varies between ∼20 km in the southwestern region to up to 1 km near the Mississippi River delta.The grid has 20 vertical layers with increased resolution near the surface and bottom.The setup and validation of the circulation model is described in detail in Hetland andDiMarco (2008, 2012).
The model uses a fourth-order scheme for the horizontal advection of tracers and a third-order upwind scheme for the advection of momentum, with conservative parabolic splines to calculate vertical gradients.Vertical mixing is parameterized using the Mellor-Yamada level 2.5 turbulent closure scheme (Mellor and Yamada, 1982;Galperin et al., 1988).Climatologies from the World Ocean Database (Boyer et al., 2005) are used to prescribe temperature and salinity at the open boundaries.Atmospheric forcing is specified using 3-h near-surface winds from the NCEP North American Regional Reanalysis (NARR) high-resolution climate data set (Mesinger et al., 2006) and surface heat and freshwater flux climatologies of da Silva et al. (1994a,b).Daily freshwater fluxes from the Mississippi and Atchafalaya rivers are based on measurements of freshwater transport by the US Army Corps of Engineers at Tabert Landing and Simmesport, respectively.

Biological model
The circulation model is coupled to a modified version of the N cycle model developed by Fennel et al. (2006Fennel et al. ( , 2008)).The original model describes the pelagic N cycle with seven compartments representing phytoplankton (Phy), chlorophyll (Chl), zooplankton (Zoo), nitrate (NO 3 ), ammonium (NH 4 ) and detrital particulate organic N (PON), which is divided into a pool of small detritus (SDet, particles size < 10 µm) that receives losses from phytoplankton and zooplankton and a pool of large detritus (LDet) that represents aggregates of phytoplankton and small detritus.A salinitydependent light attenuation coefficient accounts for increased light attenuation in the river plumes due to colored dissolved organic matter and suspended terrigenous sediment.A detailed description and validation of the N cycle model for the northern Gulf of Mexico domain is available in Fennel et al. (2011).For completeness' sake the model equations are presented in the Appendix and parameter values given in Table A1.
For this study a new state variable for DIP was added to the N cycle model.It is assumed to represent orthophosphate (hereafter referred to as DIP).The four pathways that control the dynamics of DIP are uptake for phytoplankton growth, excretion by zooplankton due to feeding and basal metabolism, remineralization of POM in the water column, and remineralization of POM in the sediment.The release of P sorbed to particulate matter when Mississippi waters enter the Gulf of Mexico (Fox et al., 1985) or the use of DOP as an alternate P source for phytoplankton growth are not accounted for in the model.A schematic of the modified model structure is given in Fig. 2.
As for ammonium and nitrate, the uptake of DIP by phytoplankton is assumed to follow Michaelis-Menten dynamics.The degree of nutrient limitation is determined by the most limiting element, either N corresponding to DIN = NO 3 + NH 4 or P corresponding to DIP, using the limitation factors L N = L NO 3 + L NH 4 for N and L P for P, which are calculated as: (1) (2) where k NO 3 , k NH 4 and k DIP (mmol m −3 ) are the halfsaturation concentrations for nitrate, ammonium and DIP uptake, respectively.Their values are set to k NO 3 = k NH 4 = 0.5 mmol N m −3 (Fennel et al., 2006(Fennel et al., , 2011) )  0.03 mmol P m −3 .The range of observed values for k DIP is large (Dortch and Whitledge, 1992).Here, a conservative value was chosen for k DIP that assumes nutrient uptake occurs in Redfield stoichiometry (i.e.k DIP = k NO 3 /16).We consider hereafter the following types of limitation to phytoplankton growth: 1. N limited: if L N < L P and L N < 0.75 2. P limited: if L P < L N and L P < 0.75

No nutrient limitation otherwise
The same criteria were used to calculate limiting factors from the observations, allowing for a direct comparison between model results and observations.The specific growth rate of phytoplankton (µ; d −1 ) is then expressed as a function of light (E; W m −2 ), temperature (T ; • C) and the most limiting nutrient such that: where µ max (T ) is the temperature-dependent maximum growth rate of phytoplankton (see Eq. A9) and f (E) is the light limitation factor (see Eq. A10).Model phytoplankton, zooplankton and detritus pools are assumed to be in Redfield stoichiometry and from assimilation of phytoplankton.Water column remineralization of detrital particulate organic matter (small and large pools) into DIP is a linear function of particulate matter concentration.The time rate of change of DIP due to biological processes in the water column is: where l BM (d −1 ) is the rate of excretion by zooplankton due to basal metabolism, l E (d −1 ) is the maximum rate of excretion by zooplankton due to assimilation (Leonard et al., 1999), k P ((mmol N m −3 ) 2 ) is the half saturation constant for zooplankton grazing on phytoplankton, β is the zooplankton assimilation efficiency and r SD (d −1 ) and r LD (d −1 ) are the remineralization rates of small and large detritus, respectively.
All sinking particulate organic P (phytoplankton and detritus) is instantaneously remineralized into DIP when reaching the sediment-water interface.This is analogous to the treatment of N at the bottom, except that a fraction of PON reaching the sediment-water interface is lost through denitrification (see Eqs. A13-A14).This fraction is fixed and was determined empirically from a relationship between sediment denitrification (representing all processes of N 2 gas production) and oxygen consumption in a range of aquatic environments (Seitzinger and Giblin, 1996;Fennel et al., 2009).The sediment-water interface parameterization assumes that denitrification occurs through coupled nitrification-denitrification only.A detailed description of the calculation is presented in Fennel et al. (2006).

Simulations
The model was run for the period 2001-2007 using monthly measurements of nutrient loading for the Mississippi and Atchafalaya rivers (Fig. 3, Aulenbach et al., 2007).This simulation is referred to as the control run, and its results are compared with satellite chlorophyll estimates and in situ nutrient observations.The control simulation is then compared to one with the original N-based version of the ecosystem without P dynamics (Fennel et al., 2011), but otherwise identical forcing, initial and boundary conditions and parameter values (referred to as the N-only simulation).Two additional model runs were carried out to investigate the sensitivity of the system to variations in DIP load.In these runs, the model setup remains the same as in the control run except that DIP concentrations in the Mississippi and Atchafalaya rivers are increased or decreased by 50 % compared to the control run.This range of variation has been used in a recent modeling study of the effect of nutrient loading in the Mississippi River plume (Eldridge and Roelke, 2010) and corresponds to the range of variability associated with nutrient reduction strategies (Eldridge and Roelke, 2010, and references therein).
For model analysis, five geographical zones were defined (Mississippi delta, Mississippi intermediate, far-field, Atchafalaya delta and Atchafalaya intermediate) and are shown in Fig. 1.The first three sub-regions correspond to an ecological gradient associated with the Mississippi River plume (Rowe and Chapman, 2002) and were used in the recent modeling study of Fennel et al. (2011).The Mississippi delta and Intermediate regions were originally defined by Lohrenz et al. (1997).The Atchafalaya sub-regions were defined for the present study to account for the ecological gradient associated with the Atchafalaya River plume.Correlations were computed between simulation results and satellite observations (SeaWiFS) of monthly surface chlorophyll concentrations over the whole domain (for both, the control run and the N-only run).Root mean square errors (RMSEs) were computed for the spatially averaged monthly mean surface chlorophyll values in the five subregions for both, the control run and the N-only run.In addition, RM-SEs and bias were calculated between simulated surface DIP concentrations and observations from three sources: Sylvan (Sylvan et al., 2006(Sylvan et al., , 2007(Sylvan et al., , 2011)), EPA (Lehrter et al., 2009(Lehrter et al., , 2012) ) and LUMCON (Rabalais et al., 1999(Rabalais et al., , 2007)).  1, including the median (black line), the range between the 25th and 75th percentile (dark grey area) and the range between the minimum and maximum value (light grey area).Also shown are observations from Sylvan et al. (2006Sylvan et al. ( , 2007Sylvan et al. ( , 2011, circles), circles), LUMCON (Rabalais et al., 1999(Rabalais et al., , 2007, squares) , squares) and EPA (Lehrter et al., 2009(Lehrter et al., , 2012, triangles), triangles).Observations are averaged by month with error bars indicating the range between the 25th and 75th percentile.

Surface chlorophyll and DIP concentrations
Time series of monthly mean simulated and observed surface chlorophyll averaged over the five areas defined in Fig. 1 are shown in Fig. 4. Selected climatological monthly mean chlorophyll fields over the entire model domain are presented in Fig. 5. Simulated chlorophyll follows a pronounced seasonal cycle in the two delta regions and in the Atchafalaya intermediate region, reaching annual maxima of 10-15 mg m −3 between June and August (Fig. 4a, b) after the peak in river discharge that occurs in late spring (Fig. 3).An annual minimum of 2-4 mg m −3 is found during the late fall and winter.The amplitude of the annual chlorophyll cycle is smaller in the Mississippi intermediate region (Fig. 4c) and almost disappears in the far-field region (Fig. 4e), where chlorophyll concentrations remain low throughout the year (< 1 mg m −3 for most of the time).Overall, chlorophyll decreases westward of the river sources and southward from the Louisiana coast.
The spatial distribution of simulated surface chlorophyll on the Louisiana Shelf is in agreement with satellite observations as shown for the period of interest (April to September) in Fig. 5.This is indicated by correlation coefficients of 0.68 on average and a maximum correlation of 0.74 in May (Fig. 5c).The time series comparison shows that the largest differences between the control and the N-only simulations occur in the Mississippi delta region (Fig. 4a).We consider a change in RMSE > 10 % to be significant.In the Mississippi delta region, error statistics indicate a significant improvement of model results in the control simulation over the N-only simulation (Table 1).RMSEs are improved for every year of the simulation except 2002.On average, RMSEs decrease by 19 % in this region (Table 1).Model results are also improved significantly in the Mississippi intermediate region where RMSEs are lower for every year except 2003 and decrease by 31 % on average (Table 1).The changes in error statistics are not significant in the other regions.
In the Atchafalaya delta and intermediate regions, which are both directly influenced by the Atchafalaya river inflow, the model systematically underestimates chlorophyll concentration in winter and early spring.This may be due in part to unresolved processes affecting chlorophyll concentration in this region, such as inputs of dissolved organic matter from wetlands adjacent to Atchafalaya Bay, and partly due to an overestimation of satellite chlorophyll related to the presence of colored dissolved organic matter in this area.
The seasonal cycles of DIP and DIN concentrations are out of phase in the Mississippi River.DIP concentration is at its annual maximum in late summer, whereas DIN concentration is at its maximum in spring (Fig. 3).However, their total loads are dominated by discharge, which is at its annual maximum in spring.Surface DIP concentrations follow the annual cycle of DIP load in the delta regions, with annual maximum concentrations of 1-1.5 mmol P m −3 for the Mississippi delta (Fig. 6a) and 2-2.5 mmol P m −3 for the Atchafalaya delta (Fig. 6b) in winter and spring, and with annual minimum concentrations of less than 0.2 mmol P m −3 in summer (Fig. 6a, b) after the peak in discharge and after the annual peak in phytoplankton biomass.A similar seasonal pattern, but of lower magnitude, is found in the intermediate regions (Fig. 6c, d).In the Atchafalaya intermediate region, DIP concentrations are generally higher than in the Mississippi intermediate region.This difference is important in determining the pattern of P limitation on the shelf.In the farfield region, DIP concentrations are low throughout the year, with a minimum in early summer (Fig. 6e), but are usually not entirely depleted.
Surface DIP observations from the three data sets fall within the range of variability found in the simulation for all data sources in the Mississippi delta (Fig. 6a), Atchafalaya delta (Fig. 6b) and Atchafalaya intermediate (Fig. 6d) regions, and with the Sylvan and EPA data sets in the Mississippi intermediate (Fig. 6c) and far-field (Fig. 6e) regions.In the Mississippi intermediate (Fig. 6c) and the far-field (Fig. 6e) region the LUMCON observations agree with the model in 2001, but not from 2002 onward.In particular during 2002 in the Mississippi intermediate region the observed DIP concentrations are significantly higher even though the same stations were sampled in all years.RMSE and bias between the simulation and the three data sets are given in Table 2. RMSEs are smallest for the EPA data set with practically no bias.For the Sylvan data set the RMSEs are higher and the model has a slight positive bias (i.e. the model predicts larger DIP concentrations than observed).The RMSEs are largest for the LUMCON data set with a pronounced negative bias (i.e. the model underestimates DIP observations).Below we use all three data sources to calculate observed nutrient limitation factors on the Louisiana shelf.

Nutrient limitation on the Louisiana shelf
Nutrient limitation on the Louisiana shelf is explored using the limitation factors L N and L P (Eqs.1-3) and the criteria defined for N and P limitation.The degree of nutrient limitation is given by the size of the limitation factors, with smaller values indicating stronger limitation and values of L above 0.75 indicating the absence of nutrient limitation.
The spatial distribution of nutrient limitation during the eight cruises described in Sylvan et al. (2006Sylvan et al. ( , 2007Sylvan et al. ( , 2011) ) and Quigg et al. (2011) is presented in Fig. 7. Data from the first four cruises illustrate the seasonal evolution of nutrient limitation during 2001 (Fig. 7a-d).Nutrients are not limiting along the Louisiana coast in March (L > 0.9).This period corresponds to the beginning of the phytoplankton bloom and high discharge with high nutrient loading from the Mississippi River (Fig. 3).Nutrient concentrations are high but phytoplankton biomass is still close to its annual minimum (Fig. 4).Nutrient limitation is found downstream, mainly from N in the far-field and from P in deeper areas south of Atchafalaya and Terrebonne Bays (Fig. 7a).In May, P limitation develops at the edge of the Mississippi River plume and near Terrebonne Bay (L P ∼ 0.5, Fig. 7b), which corresponds to the Mississippi intermediate region (Fig. 1).This is explained by the high DIN : DIP ratio in the Mississippi River during the preceding discharge period (March-May).Nitrogen limitation (L N < 0.2) occurs in the deeper shelf regions and offshore.Phosphorus limitation is still present in the Mississippi intermediate region in July and extends even further offshore, reaching maximum strength in the deeper shelf area south of Atchafalaya Bay (L P < 0.2, Fig. 7c).Phosphorus limitation has its maximum extent at this time.With the exception of the Mississippi and Atchafalaya River plumes, which are not nutrient-limited, the entire shelf is N-limited in September (Fig. 7d).
Similar patterns are found in 2002 (Fig. 7h) and 2004 (Fig. 7e-g).In March 2004, nutrients are not limiting in the nearshore areas, except for a region of weak P limitation within the Mississippi intermediate region (Fig. 7e).Nitrogen limitation occurs outside the river plume at this time.In May 2004, nutrient limitation is well developed with an extended region of P limitation at moderate levels (0.6 < L P < 0.75) near Terrebonne Bay and westward, and a region of strong N limitation (L N < 0.2) on deeper shelf regions and in the far-field (Fig. 7f).Again, the P-limited region increases in size in July 2004, covering the same area as in July 2001, except for an extended area without nutrient limitation in the Atchafalaya River plume (Fig. 7g).In July 2002, the P-limited area is smaller than in 2001 and 2004 (Fig. 7h).
The simulated patterns of nutrient limitation agree well with observations from the Sylvan data set in March (Fig. 7a), July (Fig. 7c) and September 2001 (Fig. 7d), but the model underestimates P limitation in the Mississippi and Atchafalaya delta regions in May 2001 (Fig. 7b) when observed DIP concentrations are below detection limits (DIP < 0.03 mmol m −3 ).The agreement is also good in 2002 (Fig. 7h) and in 2004 (Fig. 7e-g), with some discrepancy in the extent of the P-limited area.The main discrepancies with the observations occur in March 2004, when the model simulates the observed area of N limitation correctly but underestimates the magnitude of P limitation near Terrebonne Bay (Fig. 7e), and in July 2004, when the simulation indicates P limitation near Terrebonne Bay while weak N limitation is observed (Fig. 7g).Aside from these discrepancies, the spatial and temporal distribution of P limitation is well represented by the model.
Table 2. RMSE and bias (both in units of mmol P m −3 ) between simulated surface DIP concentrations and observations from the three data sources.For the calculation, simulated data were spatially averaged for days when observed data are available (observations were spatially averaged as well) and then error statistics are computed on these averages.Bias is calculated as model-observations; thus a positive bias indicates that the model overestimates the observations.N is the number of observations available for each calculation.The seasonal cycle and spatial gradient of nutrient limitation over the simulation period is summarized in Fig. 8.In the Mississippi and Atchafalaya delta regions (Fig. 8a, b) DIN and DIP concentrations are always high due to river inflow and thus nutrient limitation does not occur.The study of Fennel et al. (2011) using the N-only model indicates that light, rather than nutrients, is the dominant limiting factor in the Mississippi delta region.In contrast, N is almost always strongly limiting in the more oceanic far-field region (L N < 0.5, Fig. 8e).In the Mississippi and Atchafalaya intermediate regions nutrient limitation develops seasonally (Fig. 8c, d).The Atchafalaya intermediate region is Nlimited from August to September but not nutrient-limited otherwise (Fig. 8d).P limitation does not occur in this region and N limitation did not occur at all in summer 2001 and 2004.For those two years, the June DIN load from the Atchafalaya River was more than double the June load of the other simulated years (Fig. 3).This late input of DIN has alleviated N limitation westward 91 • W during summer, resulting in the high annual peak in chlorophyll (Fig. 4d).

Sylvan
In the Mississippi intermediate region, P limitation develops between May and July (Fig. 8c) following the annual peak discharge.This is also illustrated on the spatial maps presented above (Fig. 7).In late spring/early summer, waters transported from the Mississippi delta reach the Mississippi intermediate region where they are depleted in DIP, initiating a period of P limitation.This is followed by a period of strong N limitation (0.2 < L N < 0.6) from late August to November, when river inflow reaches its annual minimum (Fig. 3).At this time, the region is more influenced by oceanic waters where N limitation is dominant.Nutrient limitation does not occur in late fall and winter (Fig. 8c), but light limitation is important during this period (Fennel et al., 2011).
Phytoplankton growth is limited by either N or P but the possibility of N + P co-limitation was investigated by modifying the nutrient limitation criteria.Results are presented in the online supplementary material and show that N + P www.biogeosciences.net/9/4707/2012/Biogeosciences, 9, 4707-4723, 2012 co-limitation seldom occurs with the modified limitation criteria.Therefore, considering N + P co-limitation does not change the results described above or our conclusions in any qualitative way.

Effects of P limitation
The effect of P limitation on primary production and POM depositional fluxes is estimated by comparing the control simulation against the N-only simulation.In both simulations the mean rate of water column-integrated primary production is highest between May and July and lowest between November and January (Fig. 9a-e).In the control run, the highest rates of depth-integrated primary production occur in June with an average of 20.6 mmol N m −2 d −1 in the Mississippi delta region (Fig. 9a) and 19.1 mmol N m −2 d −1 in the Atchafalaya intermediate region (Fig. 9b).Primary production in the Atchafalaya delta region is lower than in the intermediate region due to light limitation near the Atchafalaya River outflow.These rates are lower than previously reported (Lohrenz et al., 1997;Quigg et al., 2011), although the local maximum daily rate of primary production (79.2 mmol N m −2 d −1 in July within the Mississippi delta region) is within the range of primary production values reported for the region (65.4-144.7 mmol N m −2 d −1 , Lohrenz et al., 1990Lohrenz et al., , 1999;;Quigg et al., 2011).The comparison between the control and N-only simulations indicates a significant decrease of primary production due to P limitation in spring/early summer in the Mississippi delta (26 % decrease in June, Fig. 9a) and Mississippi intermediate (12 % decrease in May, Fig. 9c) regions, but an increase during summer and early fall in the Mississippi intermediate region and westward by up to 18 % in the intermediate regions (Fig. 9c, d) and 27 % in the far-field region (Fig. 9e).In other words, a fraction of the nutrient-stimulated primary production is shifted downstream in space and delayed in time (Fig. 9f).This redistribution of primary production influences the depositional flux to the sediment.In the control run, the depositional flux is largest in June and July in the delta and intermediate regions (Fig. 10), with average values of 5.
These values are at the lower end of observed depositional fluxes which vary between 4.3 mmol N m −2 d −1 and 19.3 mmol N m −2 d −1 at the base of the photic layer in the Mississippi delta region (Redalje et al., 1994).Similar to primary production, the depositional flux is significantly lower in spring and early summer in the control run compared to the N-only simulation (Fig. 10f), mainly in the Mississippi delta (33 % decrease in June/July, Fig. 10a,f) and Mississippi intermediate (28 % decrease in May, Fig. 10c, f) regions, and is higher by up to 30 % in August and September in the Mississippi intermediate region and westward (Fig. 10f).
In the model, organic matter that is deposited to the sediment is remineralized instantaneously, as described in Sect.2.2.Organic P is restored to the bottom water as DIP, while a fraction of the organic N is assumed to be denitrified (the remainder is restored to bottom waters as ammonium).Denitrification rates vary over the course of the year.Multi-year monthly mean rates range from 1.6 to 5.5 mmol N m −2 d −1 in the delta regions, from 0.6 to 4.2 mmol N m −2 d −1 in the intermediate regions and from 0.5 to 1.3 mmol N m −2 d −1 in the far-field region.These rates are similar to the denitrification rates of Lehrter et al. (2012), who measured rates from 0.9 to 2.8 mmol N m −2 d −1 on the Louisiana Shelf.In the Atchafalaya intermediate region this N removal amounts to 37 % of primary production in June, but only to 21 % in the Mississippi intermediate region.

Consequences of altered DIP load
The sensitivity of nutrient limitation patterns to changes in the N : P ratio of the river nutrient load was evaluated by varying river DIP concentrations by ±50 %.In these scenarios, the N : P ratio peaks around 38 (for increased DIP) and 150 (for decreased DIP).The latter is an order of magnitude higher than the Redfield ratio.P limitation is most strongly affected in May and July during the peak of primary production (Fig. 11a).In the control simulation, the extent of the Plimited area increases from an average of 1.94 × 10 4 km 2 in May to 3.38 × 10 4 km 2 in July.A 50 % increase in river DIP concentration significantly reduces the area of P limitation to 1.16 × 10 4 km 2 (40 % decrease) in May and 1.95 × 10 4 km 2 (42 % decrease) in July (Fig. 11a).This reduction is compensated by an increase of the N-limited area that expands in July (Fig. 11b).Conversely, a 50 % decrease in river DIP concentration extends the P-limited area to 3.03 × 10 4 km 2 (56 % increase) in May and 5.11 × 10 4 km 2 (51 % increase) in July (Fig. 11a) and reduces the N-limited area (Fig. 11b).P limitation lasts until September in this case.
Interannual variations in river DIN concentrations, in the timing of water discharge and in circulation patterns on the shelf induce significant interannual variability in the extent of P limitation and its response to altered river DIP concentrations (Fig. 11a).When river DIN concentrations are low (2002 and2003, Fig. 3), the spatial extent of P limitation is at its minimum and becomes very small with added DIP (Fig. 11a).In the years when DIN load is highest in June (2001 and2004, Fig. 3), the P-limited area is at its maximum during July (Fig. 11a).

Discussion
Variations in primary production in the surface waters of the northern Gulf of Mexico are primarily driven by N delivery from the Mississippi and Atchafalaya rivers (Rabalais et al., 2002).The annual cycle of N and P loads controls the patterns of nutrient limitation over the shelf, while light limitation due to the presence of suspended terrigenous sediments and chromophoric dissolved organic matter is the most important factor limiting primary production in the delta regions of the Mississippi and Atchafalaya rivers (Lohrenz et al., 1990(Lohrenz et al., , 1999;;Fennel et al., 2011;Quigg et al., 2011).Realistic model simulations presented here demonstrate that in the intermediate salinity region between the light-limited and the N-limited regions there is a large zone of P limitation (3.4 × 10 4 km 2 on average) during the annual peak of primary production between May and July.This result is consistent with the high alkaline phosphatase activity and P stress measured in May and July at several locations within the Mississippi delta and intermediate regions (Sylvan et al., 2006(Sylvan et al., , 2007(Sylvan et al., , 2011;;Quigg et al., 2011).Moreover, the improved agreement with observations when including P limitation in the model demonstrates that multi-nutrient interactions are an important component of phytoplankton dynamics on the Louisiana Shelf.These results substantiate the framework of resource limitation recently proposed by Quigg et al. (2011) that relates DIN concentration (decreasing downstream of the Mississippi River) to resource limitation.Significant interannual variability in the intensity of P limitation (Fig. 8) results from variability in the timing and magnitude of river nutrient inputs.
P limitation has been suggested to induce a delay in the uptake of N in the Mississippi River plume (Quigg et al., 2011), an idea that is confirmed by the simulations presented here.The simulations demonstrate how P limitation alters the distribution of primary production and depositional fluxes on the Louisiana shelf, spreading the effect of allochthonous nu-trient over a larger area.A similar effect has been described for the Baltic Sea (Granéli et al., 1990).
The simulation also suggests that P limitation is more pronounced in the Mississippi intermediate region than in the Atchafalaya intermediate region.In the shallow Atchafalaya delta region bioavailable N is removed more efficiently by sediment denitrification (37 % of primary production in June) than in the deeper Mississippi delta region (only 21 % of primary production in June) where a larger fraction of remineralization occurs in the water column.Sediment denitrification removes bioavailable N but does not affect DIP (Caraco et al., 1990;Blomqvist et al., 2004), resulting in a net decrease in the DIN : DIP ratio ultimately eliminating P limitation.In the shallow Atchafalaya regions nutrients remineralized in the sediments are also more readily available to primary producers than in the deeper Mississippi regions.Thus the role of P limitation varies between the plumes of the Mississippi and Atchafalaya rivers.
It has been speculated that spatial shifts in POM deposition resulting from P limitation may result in a positive feedback on the development of hypoxic areas (Paerl et al., 2004;Conley et al., 2009).Hypoxic bottom waters are typically found in summer on the Louisiana shelf (Rabalais et al., 2002), during the peak of P limitation in surface waters between the Mississippi and Atchafalaya delta regions, which is also the main location of hypoxia (Rabalais and Turner, 2006).Phosphorus limitation may therefore enhance hypoxia in the region west of 91 • W. Conversely, P limitation may dilute the effect of eutrophication on bottom water oxygen (Quigg et al., 2011), thus reducing the extent of the hypoxic zone.Further investigation of this feedback is warranted.
An uncertainty in drawing conclusions about the real system from the model simulations presented here is the role of DOP fluxes from the river and DIP fluxes from the sediment.DOP may be a significant source of P for phytoplankton when DIP is depleted (Dyhrman et al., 2007), and could alleviate P limitation on the Louisiana shelf (Dagg et al., 2007).The use of DOP by phytoplankton has been observed in bioassays of the northern Gulf of Mexico (Quigg et al., 2011), but the availability of DOP on the Louisiana shelf be-tween May and July is not well known.Limited measurements indicate the presence of DOP (0.07-0.35 mmol m −3 ) in the Atchafalaya River plume where DIP is depleted (Pakulski et al., 2000;Cai and Guo, 2009).Moreover, the model does not account for adsorption and desorption of P in the sediment, a process that can affect the timing and magnitude of DIP fluxes between sediment and overlaying water column (Sundby et al., 1992;Conley, 2000), in particular due to the release of DIP under hypoxic conditions (Conley et al., 2002).Further investigation is required to estimate DIP fluxes at the sediment-water interface and to improve our understanding of the use of DOP by phytoplankton, in order to better evaluate the extent of P limitation on the Louisiana shelf.

Conclusions
Results from a 7-yr simulation of a multi-nutrient physicalbiological model of the northern Gulf of Mexico indicate that primary production is limited by P between the Mississippi and Atchafalaya river deltas from May to July.This is in agreement with several independent lines of empirical evidence.Model simulations illustrate that, as previously hypothesized, P limitation caps the uptake of allochthonous N, resulting in a delay and westward shift of a fraction of riverstimulated primary production.The consequence is a reduced POM flux to the sediment near the Mississippi delta, but enhanced fluxes in the Atchafalaya and far-field regions.

Ecosystem model equations
The biological model used in this study is the N cycle model developed by Fennel et al. (2006Fennel et al. ( , 2008) ) and was modified to include P. The model has 8 state variables that represent phytoplankton (Phy), chlorophyll (Chl), zooplankton (Zoo), nitrate (NO 3 ), ammonium (NH 4 ), DIP, small detritus (SDet) and large detritus (LDet).The bottom boundary condition is prescribed as instant mineralization of POM with a denitrification pathway.The time rates of change of the state variables due to biological processes are modeled as follows:  Parameters and parameter values are presented in Table A1.
The specific biological processes are described below.
The maximum growth rate of phytoplankton (µ max ; d −1 ) is modulated by temperature (T ) following the formulation of Eppley (1972): where µ 0 is the phytoplankton growth rate at 0 • C.
The chlorophyll content in a phytoplankton cell varies with time as the cell acclimates to the changes in light and nutrient conditions, following the model of Geider et al. (1996Geider et al. ( , 1997)).Only a fraction of phytoplankton growth (ρ chl ) is dedicated to chlorophyll synthesis: where θ max is the maximum chlorophyll to carbon ratio in a phytoplankton cell.The ratio ρ chl represents the ratio between actual and maximum possible photosynthesis.The rate of nitrification (n; d −1 ) is inhibited by light using the relationship: where n max is the maximum rate of nitrification, k E is the light intensity for half-saturated inhibition of nitrification and E 0 the threshold for light inhibition of nitrification (Olson, 1981).Particulate organic matter reaching the bottom boundary at depth H is remineralized instantly into ammonium and DIP.A fraction of the remineralization is assumed to occur through denitrification, resulting in a loss of N from the nutrient pool.DIP and the remainder of ammonium are returned to the overlying water.

Fig. 1 .
Fig. 1.Model domain, including grid and bathymetry (indicated in meters with the contour lines).The black boxes indicate the location of the five areas used during the analysis.

Fig. 4 .
Fig. 4. Comparison between monthly averaged surface chlorophyll concentrations in the simulations (control: red, N-only: grey) and from SeaWiFS (blue) for the five regions shown in Fig. 1.The light blue areas and light red vertical lines represent one standard deviation of the data from SeaWiFS and the control run, respectively.
Fig. 5. (A, B) Monthly averaged (2001-2007) surface chlorophyll concentration in the control run (A) and from SeaWiFS (B) in May, July and September.The black boxes indicate the five areas described in Fig. 1. (C) 2-dimensional histograms showing the comparison between simulated (control run) and SeaWiFS surface chlorophyll concentration from April to September.Their correlation is indicated in each panel.The correlations between the N-only simulation and SeaWiFS are indicated in parenthesis.The color scale represents the number of data pairs in each bin.The thin black line represents the 1:1 line.

Fig. 7 .
Fig. 7. Simulated (colour maps) and observed (circles) nutrient limitation on the Louisiana shelf in March, May, July and September 2001 (A)-(D), in March, May and July 2004 (E)-(G) and in June-July 2002 (H).The colour-coded scale represents N limitation (blue), P limitation (red) and no nutrient limitation (white).

Fig. 8 .
Fig. 8. Time series of monthly mean, area-averaged type and magnitude of N (open circles) and P (closed circles) limitation factors for the five regions described in Fig. 1.Squares indicate the absence of nutrient limitation.

Fig. 9 .
Fig. 9. Annual cycle of water column-integrated primary production for the five regions shown in Fig. 1. (A)-(E) Primary production in the control (black lines) and N-only (grey lines) simulations.Monthly means were calculated for the period 2001-2007.(F) Relative change in primary production due to the occurrence of P limitation.Anomalies are calculated by subtracting results from the control run (with DIP) from results of the N-only model.Positive values indicate an enhancement of primary production in the control run.

Fig. 10 .
Fig. 10.Annual cycle of depositional flux at the sediment-water interface for the five regions shown in Fig. 1. (A)-(E) Depositional flux at the sediment-water interface in the control (black lines) and N-only (grey lines) simulations.Monthly means were calculated for the period 2001-2007.(F) Relative change in depositional flux due to the occurrence of P limitation.Anomalies are calculated by substracting results from the control run (with DIP) from results of the N-only model.Positive values indicate an enhancement of depositional flux in the control run.

Fig. 11 .
Fig. 11.Total area of surface P (A) and N (B) limitation (L < 0.75) in March, May, July and September for the control simulation (grey) and for the model experiments with increased (red) and decreased (blue) river DIP.Filled bars indicate the monthly average for each year (from left to right, 2001-2007) and open bars indicate the average calculated over the whole simulation.Total model area is 14.6 × 10 4 km 2 .

Remineralization NO 3 NH 4 Phy NH4 DIP DIP SDet WATER COLUMN SEDIMENT PHOSPHORUS MINERALIZATION LDet Zoo EXCRETION N2 Denitrification
and k DIP =

Table 1 .
Annual RMSE for monthly-averaged surface chlorophyll between the control simulation and SeaWiFS time series presented in Fig.4.Values for the N-only simulation are in parenthesis.The total average of RMSEs is presented in the last column, bold values indicating a significant difference (> 10 %) in average RMSE between the models.

material related to this article is available
The bottom boundary conditions are therefore: Phy Phy z=H + w SDet SDet| z=H Phy Phy z=H + w SDet SDet| z=H+ w LDet LDet| z=H (A14)where z is the thickness of the bottom most grid box, w Phy , w SDet and w LDet are the sinking velocities of phytoplankton and small and large detritus, respectively.r ox represents the yield of POM oxidation to ammonium in the sediment.