Quantifying the impacts of stratification and nutrient loading on hypoxia in the northern Gulf of Mexico.

Stratification and nutrient loading are two primary factors leading to hypoxia in coastal systems. However, where these factors are temporally correlated, it can be difficult to isolate and quantify their individual impacts. This study provides a novel solution to this problem by determining the effect of stratification based on its spatial relationship with bottom-water dissolved oxygen (BWDO) concentration using a geostatistical regression. Ten years (1998-2007) of midsummer Gulf of Mexico BWDO measurements are modeled using stratification metrics along with trends based on spatial coordinates and bathymetry, which together explain 27-61% of the spatial variability in BWDO for individual years. Because stratification effects explain only a portion of the year-to-year variability in mean BWDO; the remaining variability is explained by other factors, with May nitrate plus nitrite river concentration the most important. Overall, 82% of the year-to-year variability in mean BWDO is explained. The results suggest that while both stratification and nutrients play important roles in determining the annual extent of midsummer hypoxia, reducing nutrient inputs alone will substantially reduce the average extent.

Introduction

The regular, seasonal occurrence of bottom-water hypoxia in the northern Gulf of Mexico has been well documented over the last 25 years.[1] Hypoxic extent, defined operationally as the area where bottom-water dissolved oxygen (BWDO) concentrations are below 2 mg L–1, has been determined primarily through midsummer monitoring cruises conducted regularly since 1985. These cruises demonstrate that the areal extent of hypoxia has increased over time, corresponding with a period of elevated nutrient loading from the Mississippi River basin.[1] This observation, along with other evidence of increasing hypoxia,[1,2] has motivated development of nutrient reductions strategies for the basin.[3,4]

The relationship between nutrient loads and the interannual (i.e., year-to-year) variability of the midsummer hypoxic area has been explored through a range of computational models. Scavia et al.[5] use a one-dimensional, simple mechanistic model that explains 45% of the interannual variability as a function of the May–June nitrogen load. These results are generally consistent with the statistical regressions of Bianchi et al.[6] who find that the May–June nitrogen load explains 47% of the variability, and of Greene et al.[7] who find that the May nitrate plus nitrite (NO) load explains 42% of this variability. In addition, studies by Turner et al.[8,9] and others[7,10] have shown that model performance is improved by also accounting for the long-term cumulative effects of nutrient loading.

It is widely understood, however, that water column stratification also affects the temporal variability of hypoxia. Stratification, which inhibits the reoxygenation of bottom waters, can be particularly strong in the northern Gulf due to warmer surface waters significantly freshened by discharges from the Mississippi and Atchafalaya Rivers overlying denser waters derived from the deep shelf. Wiseman et al.,[11] while acknowledging the importance of nutrient loading, determined that the majority of the interannual variability in the hypoxic area (for a nine-year period) could be explained using the mean river flow from the preceding eleven months. Using a more expansive data set, Bianchi et al.[6] demonstrate that 41% of the interannual variability in hypoxic area can be explained in terms of the May–June river flow.

The fact that both nutrients and flow can explain a large portion of the variability in hypoxic extent is not surprising. Bianchi et al.[6] report that these variables are highly correlated, with May–June flow accounting for 95% of the variability in May–June nitrogen load. Therefore, use of either of these variables masks the mechanistic effect of the other, leaving room for debate regarding whether eutrophication (via nutrient loading) or stratification (via freshwater flow) is the primary control on the interannual variability of hypoxia. Accordingly, Hetland and DiMarco[12] suggest that it will be necessary to separate the physical and biological causes of hypoxia to develop models with greater predictive capability; and a recent scientific assessment[2] stresses the need to include physical factors within hypoxia models. Moreover, a recent study by Murphy et al.[13] demonstrates the importance of both nutrients and stratification for predicting hypoxia in the Chesapeake Bay. This need likely extends to other systems, as Diaz and Rosenberg[14] report that there are over 400 coastal hypoxic zones around the world, and that they are controlled by both biological and physical factors (generally eutrophication and stratification, respectively). Because stratification is primarily a natural phenomenon, while nutrient loading is often an anthropogenic watershed pollution issue, the ability to distinguish between these factors is critical for managing these complex systems.

In this study, we determine the effect of stratification on the variability of hypoxia through a geostatistical regression (GR)[15,16] for BWDO. In the model, stratification is represented by metrics derived from salinity and temperature profiles, with regression coefficients determined based on how stratification explains the spatial variability in BWDO. While a portion of the interannual variability in BWDO is also explained by the stratification effect, the remaining portion must be explained by other factors, such as Mississippi River nutrient inputs. We evaluate the efficacy of this methodology using 10 recent years of midsummer cruise data (1998–2007). As such, this study generally focuses on the interannual variability of the system within its current state; and it does not focus on long-term trends in the system.

Materials and Methods

Data Description

We use data from midsummer hypoxia monitoring cruises conducted between 1998 and 2007 along the Louisiana-Texas shelf. These cruises are performed by the Louisiana Universities Marine Consortium (LUMCON), and the data were retrieved from the National Ocean Data Center.[17] To help ensure a consistent spatial envelope for this study, only locations sampled during at least 9 of the 10 years were included in the analysis (Figure 1), resulting in 61–64 monitoring locations for each year. Sampling locations were geo-referenced using the UTM Zone15 projection, and water depths were determined from a 3-arc-second digital elevation model (DEM) obtained from NOAA.[18] We use the dissolved oxygen, salinity, and temperature profile data collected at these locations. Typically, data were collected by two different instruments at each monitoring site: a Hydrolab and a Sea-Bird profiler.[19] Based on a comparison with the DEM, the Hydrolab and Sea-Bird typically reached to within zero and one meters of the sea floor, respectively, with minor variability. For this study, the BWDO values were taken from whichever instrument reached the greatest recorded depth (typically the Hydrolab). To determine salinity and temperature profiles, Sea-Bird data, which have better vertical resolution, were chosen preferentially over Hydrolab data; and when the Hydrolab reached a greater depth, we appended these additional measurements to the Sea-Bird profile.

Figure 1

Map of shelf bathymetry and monitoring sites used in study.

The stratification intensity at each sampling location was quantified using metrics derived from the salinity and temperature profiles. Because sampling intervals in the raw profile data were not uniform, the profiles were first resampled by linear interpolation to 0.2 m resolution. The difference (ΔS) between the 25th and 100th percentiles of salinity (Slo and Shi, respectively) was calculated for each location, and then the maximum salinity gradient (Sg) over a 0.4 m vertical interval, within the region of ΔS, was determined. The thickness (m) of the subhalocline (Hs), defined as the region below Sg, was also calculated. Figure S1 (Supporting Information (SI)) provides an illustration of these metrics. Use of the 25th and 100th percentiles for defining Slo and Shi was not arbitrary; a range of percentile values was considered, and the values that optimized the GR (based on BIC score, described subsequently in Section 2.4) were selected. For temperature profiles, an identical analysis was performed, and the zeroth and 75th percentiles were found to be optimal for defining Tlo, Thi, and ΔT. The maximum temperature gradient and subthermocline thickness, Tg and Ht, respectively, were also calculated. Because temperature tends to decrease with depth, while salinity increases, the 25th percentile of salinity and the 75th percentile of temperature both generally correspond to the 25th percentile of depth, suggesting that near-surface stratification conditions are less important for predicting BWDO.

The Mississippi and Atchafalaya River outfalls (Figure 1) provide the vast majority of fresh water and nutrients to the shelf. The Atchafalaya River is important because approximately 30% of the Mississippi River flow (and load) enters the Atchafalaya River near Simmesport, Louisiana.[5] Monthly flow and loading data were retrieved from USGS,[20] and we used the sum of the two river inputs. The USGS uses both an adjusted maximum likelihood estimator (AMLE) and a composite (COMP) method to estimate loads.[7] We used AMLE data because they generally performed better in the modeling analyses. However, the two data sets are highly correlated, and the choice does not greatly affect the results.

Geostatistical Model Formulation

Our geostatistical modeling approach takes into account the spatial correlation of the dependent variable (i.e., BWDO), which improves estimates of the model parameters relative to a regression based on the assumption of independent and identically distributed residuals, and allows for a more realistic assessment of model uncertainty.[16] The effectiveness of GR methods for modeling environmental phenomena has been demonstrated in studies of rainfall,[21] snow depth,[22] water quality,[23] and biospheric CO2 exchange.[15] With the exception of the last example, these works focused primarily on spatial interpolation, and less on inference of causal factors. As such, our work adds to the limited number of studies demonstrating how geostatistical modeling can be useful for confirming and quantifying causal relationships within environmental systems.

In the GR model, the dependent variable, z, is expressed as a combination of deterministic and stochastic components.[16] The deterministic component, Xβ, is the portion of z that can be expressed as a function of predictor variables and annual constants (i.e., intercepts). The stochastic component is the portion of z not accounted for by the deterministic component, and which can include both spatially correlated and uncorrelated parts, η and ε, respectively. Overall, the GR is represented aswhere z is an n × 1 vector of BWDO measurements taken at different locations and times. The k × 1 vector β includes the parameters of the deterministic component, which can be divided into annual intercepts (β) and regression coefficients (β) (eq 2). Correspondingly, the n × k matrix X includes the annual classifiers (X) and predictor variables (X) for each observation. The annual classifiers are binary values that bin samples by year. Each predictor variable is normalized to a mean of zero and variance of one over the 10 year period.

As discussed in the Introduction, a portion of the interannual variability in mean BWDO is accounted for through stratification metrics (which are included as predictor variables). The regression coefficients for these metrics are fit to the spatial variability in BWDO, based on the assumption that sites which are strongly stratified will have lower BWDO due to lower rates of reoxygenation. They are not fit to the interannual variability in BWDO because this could confound the role of stratification with nutrient loading, as they are temporally related, both being dependent on river flow. This approach does not, however, preclude the stratification metrics from explaining a portion of the interannual variability. As the average intensity of stratification varies from year to year, the predicted impact of stratification on BWDO varies proportionally, as a function of the regression coefficients. This is reasonable because years of intense stratification would be expected to have less reoxygenation and thus lower BWDO. The remaining interannual variability is accounted for primarily through the annual intercepts, which will be discussed further in Section 2.6.

The stochastic component, η + ε, was found to be well represented by the commonly used exponential covariance model with a nugget effect,[16] as follows:where the covariance (Q) between two observations (z and z) is modeled as a function of the distance between the observation locations (h) and the covariance model parameters σε2, ση2, and r. Here, r is a range parameter, and 3r is the effective range of spatial correlation. The parameter σε2 represents the variance of the portion of the residuals that is not spatially correlated (microvariability and measurement error), and ση2 represents the variance of the portion that is spatially correlated. Spatial correlations are expected to be stronger in the east–west direction than in the north–south direction due to along-shore currents.[24] To account for this phenomenon, h is scaled by an anisotropy ratio, α, which is the ratio of east–west to north–south correlation ranges.

Geostatistical Parameter Estimation

We use restricted maximum likelihood (REML)[25] to estimate the covariance model parameters and the anisotropy ratio. This method is recommended for models with strong spatial correlation and/or a large numbers of predictor variables,[26] and it has been applied in previous geostatistical studies.[15,27] The parameters are optimized by minimizing L (eq 4) with respect to σε2, ση2, r and α, which define Q (an n × n covariance matrix with elements determined from eq 3). Because the covariance between the stochastic components of the observations collected in different years is assumed to be zero, Q becomes block diagonal. I is the identity matrix.

Once Q is optimized, the deterministic model parameters are calculated in a straightforward manner:Thus, a complete set of model parameters is uniquely defined.

Geostatistical Variable Selection

The GR can be formulated for any subset of the available predictor variables. Candidate variables include the salinity stratification metrics (ln[ΔS], Shi, Sg, and Hs) and temperature stratification metrics (ln[ΔT], Tlo, Tg, and Ht). Note that Slo and Thi were omitted to avoid issues of colinearity; and the salinity and temperature differences were log-transformed to account for an apparent nonlinearity in their relationships with BWDO. Variables for water column depth, UTM easting, and UTM northing (D, E, N) were also included, along with their squares (D2, E2, N2) to allow for potential spatial and bathymetric trends. Overall, there are fourteen candidate predictor variables, yielding over 16 000 combinations of 1–14 variables to be considered. The challenge is to determine an optimal combination of predictor variables that best describes the phenomenon of interest, without resulting in overparameterization of the model.[15] In general, this is achieved by including only those variables that have a sufficiently high degree of explanatory ability, as defined by a statistical criterion. In this study, we used the Bayesian information criterion (BIC),[28] where models with different sets of predictor variables (X) are compared based on their BIC scores, which quantify the explanatory power of a model relative to its complexity, with the lowest score indicating the model with optimal balance. Assuming the residuals follow a normal distribution, the BIC score is calculated as follows:[15]

The GR model should be optimized in terms of both its covariance parameters (Section 2.3) and its selection of predictor variables. Because the covariance parameters of the stochastic component depend on the selected predictor variables and vice versa, an iterative approach was used.

Annual Mean Values

The interannual variability of midsummer hypoxic intensity can be assessed by comparing the mean BWDO concentrations between years. Annual mean BWDO values were determined using geostatistical kriging of the mean,[29] which assigns lower weights to observations from clustered sampling locations. The mean BWDO for year i is calculated aswith weights, λ, corresponding to observations, z. The weights are determined from the system of linear equations presented in eq 8, where Q  is the n × n covariance matrix for the observations from year i, with elements determined from eq 3. The scalar, ν , is a Lagrange multiplier, its square root being the standard error of the estimated mean value.

The annual mean impacts of the predictor variables (on BWDO) are also determined geostatistically. For a given year i, the mean impact of each predictor variable is calculated by applying the vector of location-specific weights, λ, to the vector of location-specific impacts, x β, where j refers to the jth predictor variable.

Model for the Annual Intercepts

Linear regression models were developed to explore potential relationships between the annual intercepts and nutrient and flow data. The response variable is the set of 10 annual intercepts, β̂, determined from the GR. Candidate predictor variables include the April, May, and June monthly loads and concentrations for NO, total Kjeldahl nitrogen (TKN), and total phosphorus, as well as monthly flows. These concentrations are flow-weighted averages, calculated by dividing monthly load by monthly flow. To account for potential changes in the system’s susceptibility to hypoxia, a linear trend across years is also included as a candidate variable.

Because of the large number of candidate variables (p = 22), many of which are correlated, and relatively small sample size (n = 10), the “elastic-net” method[30] was used to identify important predictor variables. This method constrains the magnitude of the model parameters (and shrinks some parameters to zero, thereby eliminating them) to prevent problems commonly associated with the overfitting of models.[31] The elastic-net method was implemented using the “glmnet” package for the R statistical computing program.[32] The degree to which model parameters are constrained is determined through a 10-fold cross validation routine included in this package. We used the elastic-net as an exploratory tool and to select one or more primary variables from which to create a more parsimonious regression model.

Hypoxic Area Model and Nutrient Reduction Scenarios

The information developed above can be used to predict how the intensity of midsummer hypoxia will respond (within the system’s current state) to changes in nutrient loading. However, while this study focuses on BWDO concentration, policy goals have centered on the areal extent of hypoxia.[3,4] Thus, we develop a regression model for midsummer hypoxic area (as determined by LUMCON[1,33]) using the annual mean stratification effects (Section 2.5) and one or more nutrient variables (Section 2.6). Using this model, we can explore the effect of hypothetical nutrient reduction scenarios over the 10 year study period.

The goal of the Intergovernmental Task Force Action Plan[3,4] is to reduce the five-year running average of hypoxic area to less than 5000 km2. Thus, we develop modeling results for all possible consecutive five-year periods across the 10 year study period. The consecutive five-year period which requires the greatest nutrient reduction to achieve the Action Plan goal is defined as the “critical” period.

Because of the non-negativity constraint on hypoxic area (any predictions of negative hypoxic area are treated as zero), probabilistic results are best determined through Monte Carlo simulation (eq 9).[34] Here, ỹ is a 5 × 1 vector of simulated hypoxic areas for a given five-year period, ỹ is the predicted (i.e., mean) response, u is a vector of independent samples from a standard normal distribution, and C(Σ̃) is the upper-triangular matrix obtained by Cholesky decomposition of the covariance matrix for the predictions (eq 10). In eq 10, W is the matrix of predictor variables used to develop the regression (corresponding to the observed conditions), W̃ is a matrix of predictors with reduced nutrient levels, σ2 is the variance of the model residuals, and I is the identity matrix.

For each nutrient reduction scenario, 10 000 simulations of ỹ were generated in accordance with eq 9, using different random draws of the vector u. Any simulations of negative hypoxic area were replaced with a value of zero, and the mean of each ỹ was then determined, resulting in 10 000 simulations of five-year average hypoxic area. Summary statistics for each nutrient reduction scenario were then calculated based on this ensemble of results.

Results and Discussion

Deterministic Model Parameters

The deterministic component of the GR (eq 2) was optimized in terms of the regression coefficients and annual intercepts. The regression coefficients (β̂) indicate how BWDO is related to the various predictor variables; all parameters are significantly different from zero (p < 0.05, Table 1). Because the predictor variables were normalized, the magnitudes of the parameter values demonstrate their relative impacts on BWDO. The regression coefficients for northing (N) and depth (D) suggest that BWDO concentrations tend to be lower at locations that are shallower and farther north (i.e., closer to shore). This is likely because the near-shore waters mix less with the more oxygenated and less eutrophic deep-shelf waters. As expected, the regression coefficients for the stratification metrics all indicate BWDO tends to be lower in areas of intense water column stratification. Annual intercepts (β̂) ranged from 1.66 to 3.50 mg L–1, and are discussed further in Section 3.2. For the 10 year study period, the deterministic component of the model (Xβ̂) explains 52% of the variability in BWDO across space and time. For individual years, the deterministic component of the model explains from 27% to 61% of the spatial variability in BWDO. The remaining spatial variability may be due largely to the effects of varying coastal current patterns that can influence transport of fresh water, nutrients, and organic matter. The stochastic portion of the model, which accounts for this remaining variability, has substantial spatial correlation, consistent with current patterns acting across large spatial scales. Additional model details including the covariance parameters, a test of the linearity assumption, and site specific results are provided in the SI.

Table 1

Regression Coefficients (β̂) with Standard Errors (σβ̂) for GRa

variable	β̂p	σβ̂
northing, N	–0.77	0.11
depth, D	0.83	0.12
salinity diff., ln(ΔS)	–0.43	0.11
max. salinity, Shi	–0.45	0.10
salinity gradient, Sg	–0.21	0.07
temp. diff. ln(ΔT)	–0.62	0.13
min temp. Tlo	0.69	0.15

Units are mg L–1.

Analysis of Interannual Variability

The results of the GR can be summarized on an annually averaged basis to analyze the interannual variability of the factors affecting mean BWDO concentration. Figure 2 illustrates how mean BWDO is affected by northing and bathymetry; salinity stratification (ln(ΔS), Shi, and Sg); temperature stratification (ln(ΔT) and Tlo); and annual intercepts. The combination of these effects can be used to exactly predict the estimated mean BWDO (eq 7). For illustrative purposes, however, each factor is presented as a BWDO depletion, calculated by determining each effect relative to its year of minimum impact; for example, temperature stratification was least severe in 2006 and most severe in 2004. For each factor, the standard deviation of the 10 annual impacts is calculated, providing a metric for assessing the degree to which each factor contributes to the interannual variability of BWDO, with higher standard deviations indicating greater contributions. From the figure, it is clear that stratification has a substantial impact on the interannual variability in BWDO depletion, and the standard deviation of the net stratification effect (salinity plus temperature) is 0.51 mg L–1 BWDO. The standard deviation of the northing/bathymetry effect is only 0.05 mg L–1 BWDO, which is expected because there is little interannual variation in the sampling locations used in this study. The greatest portion of the interannual variability is accounted for through the annual intercepts, with a standard deviation of 0.71 mg L–1 BWDO. While these intercepts do not have intrinsic explanatory value, they can be modeled in terms of other factors, such as nutrient inputs.

Figure 2

BWDO depletions attributed to different factors from GR model; each factor presented relative to its year of minimum impact.

The annual intercepts (β̂) were first related to candidate predictor variables through development of an elastic-net model (Section 2.6). The resulting model (eq 11) includes six selected variables. Here, variables enclosed in brackets represent nutrient concentrations, while unbracketed variables represent loads, and subscripts represent months. Because variables were normalized, the regression coefficients indicate their relative importance.

The elastic-net model indicates May NO concentration has the largest impact, followed by May NO load, which is discussed more below. The trend with time (Year) also has a relatively large effect, suggesting that hypoxia is becoming more intense irrespective of seasonal nutrient inputs and stratification. This is consistent with previous studies indicating increasing hypoxia irrespective of seasonal nutrient loading[7−10] or stratification,[35] but this is the first study to suggest this trend while considering both factors in combination. June NO concentration and April TKN and NO loads are also included, but appear to have relatively small effects. The selection of TKN only for the month of April is consistent with this more refractory nitrogen fraction requiring additional time to become biologically available. The elastic-net model did not select any of the variables for TP or flow, suggesting that these factors are less important for predicting midsummer hypoxia. Overall, the elastic-net model explains 85% of the variability in the annual intercepts. Figure S4 (SI) graphically illustrates how the various factors in eq 11 affect the annual intercepts.

We also developed a simple linear regression (SLR) for the annual intercepts using May NO concentration alone (eq 12). The NO concentrations were normalized so that the model coefficients are comparable to those in eq 11. This model explains 76% of the variability in the annual intercepts, indicating that its performance is similar to the elastic-net model, despite its relative simplicity.

The use of May (or May–June) nitrate data to predict hypoxia has precedent in multiple previous studies,[5−10] and is reasonable given the known processes of phytoplankton growth and bacterial decomposition that link nitrogen inputs to BWDO depletion over time.[1,11] However, while previous studies generally indicate that load is the best single nutrient predictor of hypoxia, this study indicates concentration. (Greene et al.[7] also consider models using flow and concentration, but the confounding correlation between flow and load was not addressed.) To explore this distinction further, Figure 3 illustrates how May NO loads and concentrations correlate with mean BWDO (eq 7) and the annual intercepts. As shown, mean BWDO is somewhat correlated with NO load (Figure 3a), and as expected, the strength of this correlation (R2 = 0.41, p = 0.045) is similar to that of previous studies relating nitrogen loads to hypoxic area.[5−10] The correlation between the annual intercepts and May NO concentrations (Figure 3d), corresponding to the SLR (eq 12), is much stronger (R2 = 0.76, p = 0.001). This improvement can be attributed to the fact that here we use nutrients to predict only the variability in hypoxia that remains after accounting for stratification (i.e., the annual intercepts), rather than the total variability in hypoxia (i.e., mean BWDO or hypoxic area). If we do not remove stratification effects, nutrient load is a better predictor of hypoxia than concentration (Figure 3a versus 3b). This is as expected given load is correlated with flow, and higher flow tends to increase stratification, making load partially reflective of both stratification and nutrient concentration.

Figure 3

Annual mean BWDO (top) and annual intercepts from GR (bottom) versus May NO loads (left) and concentrations (right), with 95% confidence intervals.

While these results suggest that nutrient concentration is the best predictor of hypoxia (after accounting for stratification), nutrient load may also be of some predictive importance. The elastic-net model (eq 11) did include multiple nutrient loading terms, though they generally had less impact than the concentration terms. Similarly, Figure 3c indicates that there is some correlation between the annual intercepts and May NO load, though markedly less than with concentration.

From a mechanistic perspective, nutrient load is an intuitive predictor of hypoxia based on the assumption that it is proportional to the primary production and oxygen demand that develops within the system. However, the importance of concentration (as suggested by this study) indicates that the amount of oxygen demand generated on the shelf is also affected by the degree to which that load is diluted by flow. The exact reasons for this are not readily apparent and could benefit from further research. We note, however, that loads are highly correlated with flows, and high flows can intensify coastal current velocities,[36] thus reducing the time available to produce and settle organic matter within the study area. Also, because biochemical process rates are typically concentration dependent, dilution can reduce the rates of oxygen demand formation and nutrient recycling.

The approach outlined in this study made it possible to isolate the effects of stratification through GR, and then determine the effects of nutrients based on the interannual variability in BWDO that remained. Overall, the results suggest that both river nutrient concentration and stratification play large and comparable roles in determining the interannual variability of hypoxia in the northern Gulf of Mexico. The standard deviation of the net stratification effect is 0.51 mg L–1 BWDO, while the standard deviation of the May NO effect (from the SLR) is 0.62 mg L–1 BWDO, indicating nutrients may play a slightly larger role. When the annual intercepts predicted by the SLR are used in combination with the net stratification and northing/bathymetry effects, the overall model explains a large majority (82%) of the interannual variability in mean BWDO.

Because of its explanatory power, this approach provides a more precise way to evaluate the response of Gulf hypoxia to changes in Mississippi River nutrients. Scavia et al.[5] note that when hypoxia is modeled using nutrient loads alone, the interannual variability in oceanographic conditions can mask the effects of changes in nutrient loading, at least in the short term. This can confound management of the system by obscuring the effect of watershed management practices aimed at curtailing nutrient inputs. However, by using the approach outlined here, we separate the effects of stratification from the effects of nutrients, allowing for a more detailed assessment of how each is impacting hypoxia over time. Although the intensity of stratification can be difficult to predict in advance because it varies over short time scales due to wind forcing,[36] it can likely be approximated using river flow and climate data. For example, the May river flow alone explains 27% of the variability in the net stratification effect determined by this study. The stratification effect can also be determined following each midsummer monitoring cruise using the measured salinity and temperature profiles, as was done here.

While May NO concentration and load were not well correlated over the period of this study (R2 = 0.14); longer-term trends in nitrate concentration and load have been shown to be parallel.[37] Thus, if the existing Action Plan is successful at reducing NO load, it will likely also be successful in reducing concentration. So, while we think that concentration should receive more attention in future studies and in evaluating year-to-year changes in hypoxic severity, we do not think this emphasis on concentration is in discordance with the long-term Action Plan goal of reducing nutrient loads.

Nutrient Reduction Scenarios

Stratification metrics and river nutrient concentrations were also used to develop a regression model for hypoxic area (y, km2), as described in Section 2.7. In the model (eq 13), S is the aggregated BWDO depletion due to salinity and temperature stratification (as illustrated in Figure 2); and [NOMay] is the May NO concentration. Both of the variables were normalized so that the regression coefficients indicate their relative importance. Overall, the model explains 79% of the variability in observed hypoxic area.

We used the model to estimate the average hypoxic area for different nutrient reduction scenarios, along with associated uncertainty determined through Monte Carlo simulation (eq 9); and to estimate the percent nutrient reduction required to consistently (over the 10 year study period) reduce the five-year running average of hypoxic area to less than 5000 km2 in accordance with Action Plan goals.[3,4] Figure 4 illustrates how nutrient reductions will affect the five-year running average of hypoxic area for the critical period (1998–2002 is the critical period because it requires the greatest nutrient reduction to achieve the Action Plan goal). The model suggests that a 42% reduction in May NO concentration is required, with a 90% confidence interval of 29–62% reduction. The non-negativity constraint (Section 2.7) had only a minor impact, resulting in a recommended nutrient reduction of 42% instead of 41%. These results support multiple previous modeling studies suggesting nutrient reductions of around 40–45%,[5,7,38] but are somewhat lower than the 71% reduction recommended by Liu et al.[10]

Figure 4

Predicted average hypoxic area for the critical five-year period (1998–2002) under different levels of nutrient reduction, with 90% confidence intervals.

Note that these predictions assume the linear relationships in eq 13 remain valid outside of the range of observed conditions. If the effects of stratification become less severe under scenarios of reduced nutrient inputs, this would make the estimated loading reductions somewhat conservative. It should also be noted that this study used recent data (1998–2007), so the results generally reflect the system at its current state. If the system becomes more or less susceptible to hypoxia over time, as a result of long-term trends in nutrient loading, climate change, or other factors, then the reductions required to meet the Action Plan goal would also change.