Multivariable statistical regression models of the areal extent of hypoxia over the Texas–Louisiana continental shelf

Observations of the areal extent of seasonal hypoxia over the Texas–Louisiana continental shelf from 1985 to 2010 are correlated with a variety of physical and biogeochemical forcing mechanisms. Significant correlation is found between hypoxic area and both nitrogen load (r2 = 0.24) and east–west wind speed (r2 = 0.16). There is also a significant increasing trend in the areal extent of hypoxia in time; a linearly increasing trend over the entire record (r2 = 0.17), a step increase in area for the years 1994 and beyond (r2 = 0.21), and a step increase for 1993 and beyond (r2 = 0.29) were all found to be significantly correlated with area. The year 1988, often included in other studies, was found to be a statistical outlier, in that the statistical regression properties are strongly modified when this year is included. The exclusion of any other year does not have as great an effect as excluding 1988 from the record. The year 1989 is also excluded, as this year had no full shelf survey, for a total of 24 years of data for the record. Multivariable regression models using all possible combinations of the forcing variables considered were calculated. The best performing models included east–west wind, either a linear trend in time or step in time (1994 and beyond), and either nitrogen load or river discharge combined with nitrogen concentration. The range of adjusted correlation coefficients ranged from r2 = 0.47 to 0.67. The best model (east–west wind, a step increase in time 1994 and beyond, river discharge, and nitrogen concentration) has a standard error of 3008 km2.


Introduction
The areal extent of seasonal hypoxia over the Texas-Louisiana continental shelf has been monitored every year since 1985 as part of an annual shelf-wide survey in late July (Rabalais et al 2007). Measured hypoxic areas ranged from 4400 to 22 000 km 2 from 2000 to 2010, with an average value of 15 370 km 2 during this period. The Mississippi River/Gulf of Mexico Watershed Nutrient Task Force Action Plan (2008) calls for nutrient reductions in the Mississippi River watershed such that the average areal extent of hypoxia is reduced to a running 5 year mean of 5000 km 2 . Further details on this objective, and a general description of the present knowledge of the factors influencing northern Gulf of Mexico hypoxia are found in a report by the EPA's Science Advisory Board Hypoxia Advisory Panel (2008).
There is a significant correlation between the total nitrogen load delivered to the Gulf of Mexico via the Mississippi and Atchafalaya Rivers and the areal extent of hypoxia (Justić et al 2007, Bianchi et al 2010 andreferences therein). This relationship has been exploited by a number of statistical regression models geared at predicting the extent of seasonal hypoxia, as measured in July, given the estimated nitrogen load in May, see table 1. These statistical models have, in turn, been used to estimate the required reductions in nitrogen load that would achieve mandated levels of average hypoxic areal extent.
A statistical relationship between the areal extent of hypoxia and the Mississippi streamflow-averaged from  Wiseman et al (1997) August-June average streamflow r 2 = 0.6 (N = 9) Scavia et al (2003) May-June average total nitrogen load r 2 = 0.45 (N = 18) Turner et al (2006) June NO 2+3 load r 2 = 0.53 a (N = 20) Greene et al (2009) MayNO 2+3 load r 2 = 0.42 (N = 22) This study May NO 2+3 load r 2 = 0.24 b (N = 24) a Higher correlations (r 2 = 0.60) were observed with respect to phosphorous load. b Excludes data from the years 1988 and 1989. August to June, the eleven months preceding the hypoxic area measurement in late July-was first noted by Wiseman et al (1997). They reported a correlation of approximately r 2 = 0.6. 3 At that time, only nine years were considered, including the drought year 1988 and the flood year 1993. A linear relationship between nitrogen load and hypoxic area has most recently been verified statistically by Greene et al (2009), with the highest correlations between May nitrogen load and observed (late July) hypoxic area. They also report a number of other important facts regarding the temporal character of the nitrogen load-there has been no significant trend in streamflow or nitrogen concentration. This statistical stationarity stands in contrast to an apparent increase in the hypoxic area per unit nitrogen load post-1993, attributed to an increased sensitivity of the system to nitrogen loading also noted by Turner et al (2008) and Liu et al (2010).
A series of papers (Scavia et al 2003, Scavia and Donnelly 2007, Donner and Scavia 2007) use a Streeter-Phelps model to relate nitrogen load to hypoxic area. The Streeter-Phelps model used by Scavia et al (2003) has four parameters that need to be specified or fit to obtain the relationship between nitrogen load and hypoxic area. The specification of these parameters in the model do not correspond directly to processes acting in the ocean. As with all linear models, every parameter that is used to better fit the model to observations reduces the degrees of freedom of the model by one. In essence, however, the Streeter-Phelps model is very similar to the univariate linear regression models that correlate hypoxic area to nitrogen load.
There have also been a number of studies that examine the changes in hypoxic area with time. Turner et al (2006) relate the areal extent of hypoxia to a linear regression based on both nitrogen load and time, and suggest this represents an increase in the sensitivity of the system to nitrogen load. The addition of time as a variable in the regression significantly increases the ability of the model to reproduce prior observations. Greene et al (2009) relate the trend of increasing area with time to an 'epoch', where area measurements taken post-1993 have a mean hypoxic area that is 6450 km 2 greater compared to pre-1993 measurements. Liu et al (2010) add yet another variable to the Streeter-Phelps model that represents the sensitivity of the hypoxic area to a given nitrogen load.
Wind has been observed to be an important influence on hypoxia in other estuarine systems (Møhlenberg 1999, Scully 2010a, 2010b. Turner et al (2008) exclude a number of years from their calculations because of severe storm activity at or just prior to the annual shelf-wide area survey. Variations in wind are also an important driver of seasonal circulation in the northern Gulf of Mexico (e.g., Cochrane and Kelly 1986, Cho et al 1998, Morey et al 2003, Nowlin et al 2005, Morey et al 2005. At shorter timescales, winds can shift the position of the Mississippi and Atchafalaya river plumes (e.g., García Berdeal et al 2002, Hetland 2005, Hetland and DiMarco 2008, excite inertial oscillations, particularly in summer (Zhang et al 2010(Zhang et al , 2009, and enhance mixing. As we will demonstrate below, non-storm seasonal wind patterns predict a significant fraction of the interannual variability in hypoxic area. The goals of this letter are to examine the relationship of hypoxic area to wind, nitrogen load, and time using a multivariable framework. To date, most studies examining statistical relationships have examined only univariate statistics, with the exception of the regression models presented by Turner et al (2006), Turner et al (2008) and Liu et al (2010) who examine the influence of two variables concurrently: nitrogen load and time, and Greene et al (2009), who present an stepwise multiple regression analysis on the effect many different variables on hypoxic area. This work builds upon Greene et al (2009), in that we do not consider all of the possible variables, but select a small set of variables that represent different physical or biogeochemical processes and are known to affect seasonal hypoxia. In this letter, we build on all of these previous studies by allowing a multivariable model framework to select the best combination of parameters. A Bayesian Information Criterion is used to evaluate the improvement to the regression such that the optimal model configuration may be determined. We also reexamine some previously published statistical relationships in the context of identified outliers within the data set.

Methods
Images of the areal extent of hypoxia and figures of the hypoxic area timeseries are available online 4 . We used the timeseries from 1985, with the exception of 1988 and 1989. The 1989 estimate of the hypoxic area was made from other field data during an August cruise. 1988 was a drought year, with an abbreviated cruise. While 1989 is not included in many other studies (e.g., Greene et al 2009), 1988 has always been included. As discussed in more detail below, including 1988 has a profound effect on many of the statistical relationships and therefore must be considered an outlier if the prior assumptions of the linear regression model are to be met.
Nutrient and streamflow data are available from the United States Geological Survey online 5 We used the AMLE load estimates because the data are available for all years since 1985. The composite data (suggested to be more accurate by Greene et al 2009) are missing three years. Since there is a strong relationship (r 2 = 0.84, p = 5 × 10 −12 ) and no bias between the two estimates, we thought it more important to include as many years as possible in the calculation. We computed nitrogen concentration as a ratio of the load to the basin discharge. We use the May mean nitrogen load and discharge for two reasons. First, the mean nitrogen load for May and May-June have both been demonstrated to have the highest correlation with hypoxic area (Greene et al 2009, Bianchi et al 2010. May nitrogen load is chosen as the primary nutrient variable here, as it is used for predictions of late July hypoxic area in June when the data become availableexcluding June allows the prediction to be made one month earlier. We do not consider other nutrients in this study. The statistical relationships between wind and seasonal hypoxia have not been previously investigated. Wind data observations are available from the National Data Buoy Center online 6 . We used the continuous observations from Station BURL1-Southwest Pass, LA, filling missing wind velocity data with hourly data and bivariate linear models for each of two remote NDBC stations: 42040-Mobile AL South, and 42007-Biloxi MS. We first aggregate the station observations to daily mean values and then select and aggregate a time window covering June 15 through July 15 for the eastward component of the wind, U , the northward component of the wind, V , the eastward component of the wind stress, U |U |, and the mean wind power magnitude, ( √ U 2 + V 2 ) 3 . Both wind and wind stress along-shore (east-west) wind are major drivers of along-shore currents (e.g., Cochrane andKelly 1986, Cho et al 1998), and the magnitude of the average wind stress has been shown to strongly influence the position of the Mississippi/Atchafalaya River plume (Hetland, unpublished). Thus, wind may influence the stratification envelope (the extent of the stratified region, see Hetland and DiMarco 2008) and may thereby influence the areal extent of hypoxia by increasing or decreasing the area where stratification is sufficient for hypoxia to develop. The wind power is related to vertical mixing by the wind. Although it is well known that strong winds may ventilate hypoxic bottom waters, no significant relationship between wind power and hypoxic area was found. This may be because the strong winds associated with vertical mixing events happen on short timescales that are not well captured by monthly means.
Sea surface temperature (SST) anomalies have been shown to have an impact on brown shrimp populations in the northern Gulf of Mexico (Li and Clarke 2005), and thus it is possible that other biogeochemical processes may be affected. The effect of sea surface temperature on seasonal hypoxia has also not been previously investigated. Global sea surface temperatures and anomalies are available from the NOAA extended reconstruction sea surface temperature (ERSST) online 7 . We aggregated these by month and after screening the ERSST monthly averages for univariate correlation with hypoxic area, we present results for July near the location 92 • W, 28 • N. Below, we show that no significant relationship (for July or other months) was found between SST and hypoxic area.
We fit univariate linear models relating the observed hypoxic area to the above parameters for years excluding 1988 and 1989 and tabulated model fit statistics including the standard error of the residuals (ŝ), coefficient of determination (R 2 ), adjusted R 2 (R 2 adj ), Bayesian Information Criterion (BIC), and Pearson correlation coefficient (r ). In addition, we compare the BIC statistics to the BIC of the null or intercept model with δBIC; univariate relationships with δBIC < 0 are considered to perform better than the null model. Since the univariate models have the same complexity, or number of terms, these rankings are all monotonically related to the mean error.
We fit all subsets of the above parameters as multivariable linear models, ranking them by R 2 , R 2 adj and δBIC to identify the best subsets of predictor variables. Within each subset size class, the rankings are monotonically related to each other, however, the complexity penalty differs for each statistic. R 2 adj scales R 2 by a factor of (n − 1)/(n − p − 1) to account for the degrees of freedom consumed by the number of fitted model parameters ( p) and the sample size (n). The Bayesian Information Criterion, assuming the model errors are normally distributed, is BIC = n ln(σ 2 ) + p ln(n), where,σ 2 , is the variance of the model errors. For the selected models, examinations of the distributions indicated that the model errors were not significantly non-normal, indicating that a linear model is a reasonable choice.
We applied the multivariable linear models to predict the detectable effects of a significant nitrogen load reduction on the hypoxic area after Greene et al (2009). Rather than a Monte Carlo simulation to predict detection time, we used a two-sample, one-sided t-test, t crit(0.05,n 1 +n 2 −2) = A 1 −A 0 s model , to determine the critical marginal hypoxic area reduction (A 1 − A 0 ). We estimate the requisite marginal reduction in nitrogen load or concentration by dividing by the corresponding model parameter.
The data used in this study are assembled in table 2. The table includes 26 years from 1985 to 2010, however 1988 and 1989 are excluded from the model fitting. The 1989 cruise was truncated, and is flagged as no-data 8 . The year 1988 is demonstrated to be an outlier in the analysis below.

Results
Investigations into the influence of each data point on the regression statistics show that the year 1988 has a unique and profound impact on the regression and associated correlation coefficient (figure 1). 1988 had the lowest flow, lowest nitrogen load, and the smallest hypoxic area during the  period of shelf-wide hypoxic area surveys; it was also a year of strong downwelling westward winds (i.e., winds unfavorable for hypoxia formation). Calculation of the correlation coefficients excluding different years indicates that excluding 1988 always results in an extreme value of the correlation coefficient, and is often well outside of the range of correlation coefficients obtained when excluding other years. This tendency for the year 1988 to more dramatically influence the correlation coefficient is stronger when the datum goes against the tendency of the regression for a particular property (e.g., summer mean east-west wind strength, figure 1, fourth panel). Including 1988 in a regression increases the correlation coefficient between hypoxic area and both nitrogen load and streamflow by r 0.1, nitrogen concentration by r 0.15, and reduces the magnitude of summer mean east-west wind stress correlation by r 0.25. The survey period in 1988 was also later than for any other year, Aug 12-16, as opposed to all other surveys which all ended before July 30. 9 Typically, cruises start July 21-23 and end July 25-29. For these reasons, and although this year was included in other studies, we exclude 1988 from our calculations. 1989 was also a partial, late season survey, and is typically excluded from regression calculations.
The univariate model fitness statistics in table 3 indicate that river discharge, nitrogen concentration, the northward component of wind (V ), wind power, and SST anomaly are not individually predictive. All time-related variables, Year, post-1992Year, post- (i.e., 1993Year, post- and beyond), and post-1993Year, post- (i.e., 1994 and beyond) are significant, with the post-1992 and -1993 providing a slightly lower standard error than year. The strengths and weaknesses of each of these temporal variables is discussed further below. We find, as do many other studies, that nitrogen load is also significantly correlated with Figure 1. The correlation coefficient for the regression between hypoxic area and a particular property (year, nitrogen load, streamflow, east-west wind, and nitrogen concentration) is shown in each panel. The panels show the correlation coefficient calculated by excluding a given year from the calculation. Outliers indicate that a given year has a disproportionate effect on the correlation coefficient between hypoxic area and that particular property. The thick gray line at ±0.344 indicates the one-sided 95% significance limit for 22 degrees of freedom. The year 1988, highlighted by the open circle, is identified as an outlier. hypoxic area. Of particular note here is that river discharge is not significantly correlated with hypoxic area, although this correlation was estimated to be significant by Wiseman et al (1997), Greene et al (2009), andBianchi et al (2010). This is because of the exclusion of 1988; if 1988 is included the river discharge regression becomes significant in our calculations as well. Univariate regressions confound or mask the correlations between possible predictor variables.
Multivariable regression models on the 11 predictor variables in table 3 using all possible (2 11 = 2048) subsets were fitted and compared using r 2 , r 2 adj , and BIC-metrics that measure model skill relative to model complexity as described in the methods section. Results for the top models in each Figure 2. The 'best' models for various subset sizes as ranked by Bayesian information criterion (δBIC) is shown. Thus, the best univariate model (subset size of one) with the lowest δBIC is NO x load (L), the best two-variable model (subset size of two) is a model that depends on both epoch and NO x load (E-L), and so on. δBIC increases with model error and model complexity, with better fitting, simpler models having lower δBIC. Thus, lower δBIC indicates a model that better fits the observations. Here, δBIC is defined as the improvements in BIC relative to the null model of hypoxic area = 14 km 2 . Not shown on this graph are the scores for the 2 11 − 7 = 2041 models of similar complexity but higher residual error, i.e., higher δBIC scores.
subset size class are shown in figure 2. The parameters and metrics of four selected models, (the best 3-and 4-term models, and the models which use 'Year' instead of the post-1993 epoch), are shown in table 4. Hindcasts of hypoxic area 1985-2010 using these models are shown in figure 3. These four models appear to give qualitatively similar predictions. However, the models using the post-1993 epoch outperform the models which include Year, while the four-term models using river discharge and nitrogen concentration outperform the three-term models using nitrogen load. It is interesting to note that the 1988 and 1989 partial estimates of hypoxic area fall within the normal range of this model. For the case of the year 1988, which was excluded from training the models due to its high influence on the model parameters, the abbreviated survey, and the late survey time, this seems to suggest that it was not an unusually small value for hypoxia, but rather the unusually strong westward wind and unusually weak river discharge may explain the unusually low value.
The calculation of all subsets of multivariable regression selects a different set of optimal model variables than might be suggested by the univariate regressions. For example, although post-1992 was the variable with the highest correlation with hypoxic area, post-1993 (epoch) is the variable selected in all of the subsets larger than two. This is because 1993 had a very large hypoxic area, which causes the post-1992 epoch to fit the data better when considered alone. Without considering nutrient load, the additional point in the step function in 1993 (the post-1992 epoch) better fits the large area in that year. However, 1993 also had a very streamflow and nutrient load, so these other variables can explain the large hypoxic area in  1993 when they are also included in the regression without requiring the step function to begin in 1993. That is, when other variables are included, the post-1993 epoch better fits the data. Also, although eastward wind stress gives a higher univariate correlation with hypoxic area than eastward wind, it is eastward wind that is selected in the optimal model. Finally, neither river discharge nor nutrient concentration have significant (univariate) correlations with hypoxic area, however they appear in combination in the overall best-fit model of epoch, discharge, nitrogen concentration and wind. Of course, the product of mean discharge and nutrient concentration is nutrient load, which has a relatively strong univariate correlation with hypoxic area. We have calculated the amount of nitrogen reduction required in terms of nitrogen load (table 5) or nitrogen concentration (table 6) to achieve a statistically significant hypoxic area reduction using a two-sample, one-sided t-test for a number of different time horizons. In all cases, larger reductions in nitrogen load or concentration will result in a detectable reduction in hypoxic area sooner. In addition to the calculations presented above, we also ran a series of regressions using the log of the area (not shown), with similar results to the linear response model, but the residuals appeared to be non-uniformly distributed and the correlation coefficients were smaller. The log-transformed models do have the conceptual advantage of a zero lower bound in hypoxic area, such that low discharges or nitrogen loads cannot produce negative hypoxic areas. Finally, we performed preliminary investigations on possible memory effects in the system (not shown). There is no significant correlation between areal extent of hypoxia in a given year and either previous year's nitrogen load (integrated over any time window within that year) or the previous year's hypoxic area. Thus, as suggested by Wiseman et al (1997), Turner et al (2008) and Liu et al (2010), memory effects may exist, but the processes are more complicated than can be explained by a simple regression model.

Discussion
Using the all-subsets BIC rankings in figure 2 suggests that the best model is a four-term model containing the eastwards component of wind, the post-1993 Epoch, the Mississippi River basin discharge, and the concentration of NO x . The more parsimonious three-term model including wind, epoch, and nitrogen loading is nearly as predictive, as the product of river discharge and nitrogen concentration equals the nitrogen load. However, it is interesting to note that the model improves when the discharge and concentration are considered separately. The fact that model results are more accurate using both discharge and concentration as inputs implies that there are likely two effects of the discharge on hypoxic area: one the combined effect of discharge and concentration on nitrogen load influencing eutrophication, the other the effect of discharge alone influencing stratification. However, because discharge and nitrogen load are highly correlated, it is not possible to separate these effects in a statistically significant way using univariate regression models.
There is a significant increasing temporal trend in hypoxic area, either in the form of a linear increase in area, or a step increase after 1993. Using epoch, a post-1993 step function, to represent the increase in time results in a more accurate hindcast of observed hypoxic areal extent in the multivariable regressions when nutrient load is also considered. However, epoch does not predict hypoxic area significantly better in a statistical sense than using a linear temporal trend, and so at this point we are unable to exclude either class of temporal trend from consideration. Given that it is impossible to statistically distinguish between a step and a linear trend in hypoxic area, there are practical reasons to prefer using a step to reproduce temporal trends. A linear increase in hypoxic area with time impairs the ability of the model to be used for prescribing nitrogen reduction targets. For example, the regression model by Turner et al (2006) 10 is neutral, in that zero nitrogen load would result in zero hypoxic area in 1990; in 2010, the minimum area, with zero nitrogen load, is 13 086 km 2 . Thus, it is impossible to achieve the mandated reduction to 5000 km 2 in the future within the framework the Turner et al (2006) regression model. The best approach for forecast models based on historical regressions may be to include a (nonlinear) sensitivity factor, as done by Liu et al (2010). This approach was not followed here, because we have chosen to only consider linear models.
The year 1988 was identified as an outlier that was seen to affect the statistical regression models much more than any other year. It is possible that the late August survey was during a period when hypoxic area was in decline, and the full extent of the maximum area was not measured during the abbreviated survey. It is also possible that the extremely low flow during 1988 shifted the system to a different state, such that linear relationships between forcing variables and hypoxic area are no longer valid. One possible mechanism that could be responsible for such a switch is the combination of weak flow and strong downwelling winds that also occurred that year, which confined the Mississippi and Atchafalaya river plumes to very near the coast, inshore of the region of the shelf where hypoxia is most often found. Wind effects on buoyant plumes have been the subject of a number of previous studies (Fong andGeyer 2001, Hetland 2005), and the Mississippi/Atchafalaya plume system is known to respond strongly to wind forcing (Hetland andDiMarco 2008, Wang andJustic 2009).
The multivariable regression models can be used to estimate nutrient reductions required to achieve the stated goal of a 5 year running mean hypoxic area of 5000 km 2 or lower. In order to do this, we will use the two best-fit models using epoch in table 4, using mean values for the variables not associated with nitrogen. As discussed above, the two models with a linear trend in time (with a year variable) are not useful when predicting future areas. For example, the UYDC model (see table 4  present state. It is not presently clear what caused this shift in the system, or if it will return to its previous state in the future. Thus, nutrient reduction estimates must consider both states to obtain a range of nutrient reduction estimates. The UEDC model requires a 68% reduction in nutrients in the present epoch state (E = 1), a 29% reduction in the previous epoch state (E = 0), to achieve a hypoxic area of 5000 km 2 or lower. This model estimate has the advantage that river discharge (an unmanageable variable) keeps its mean value, and nitrogen is altered independently. It is not possible to achieve a hypoxic area of 5000 km 2 or lower using the UEN model in the present epoch state (E = 1), for the previous epoch state (E = 0) a 54% reduction would be required. Care should be taken when interpreting these results, however, as the nutrient reduction scenarios are well outside of previous system states. The standard deviation of nitrogen concentration is 18% of the mean, nitrogen load 28% of the mean. Thus, all of these nutrient reduction estimates represent shifts of two standard deviations or more in the nitrogen variables. Such a strong perturbation will very likely induce a nonlinear system response that is unpredictable using a linear model.
There is an important question of whether future nitrogen reductions will result in changes to the shelf ecosystem that will be observable in a statistical sense. We find that it would take approximately 5 years to say with some certainty that a 25% reduction in nitrogen load resulted in a statistically significant reduction in hypoxic area (see table 5), similar to Greene et al (2009) who suggest that a 30% reduction in nitrogen load would be observable with 50% certainty after 5 years. Interestingly, for a model that separates streamflow and nitrogen concentration, much smaller reductions in nitrogen concentration are required to create a statistically significant reduction in area (see table 6). Note that our calculations require knowing the system response to the other forcing variables. As mentioned above, these nutrient reductions exceed the natural variability within the system over the observational record, and it is not clear if the response of the system to all of the variables will remain the same with such a large perturbation to even just one of the forcing variables. The estimates of nitrogen loading reductions needed to achieve the target hypoxia are made assuming that the relationships observed in the data hold at nitrogen concentration much lower than the range of observations. The actual hypoxic area associated with the reduced nitrogen loading could potentially be much different from the prediction.

Conclusions
Statistical regression analysis shows that there are four key forcing variables in predicting the areal extent of seasonal hypoxia on the Texas-Louisiana continental shelf: average May Mississippi/Atchafalaya riverine nitrogen concentration, average May combined river discharge, mid-June to mid-July average east-west wind speed, and a shift in the system post-1993. Alternative models include a three-variable system that substitutes nitrogen load for discharge and concentration, and substituting a linear trend for the post-1993 temporal step.
Automated comparison of many regression models enables identification of key forcing variables from a larger set of interacting candidate variables. This procedure is extensible to examine other candidate forcing variables and their relative influence on hypoxic area. Modeling the effects of exogenous forcing variables on hypoxic area helps build understanding of the effects of policy-controllable variables and extends the applicability of models to more varied conditions. Previous models of areal extent have used nitrogen load, which is the product of discharge and concentration; we find this simplified model is similar to the four-variable model that includes discharge and concentration as separate terms. Our results indicate that larger per cent reductions are required to achieve an average hypoxic area of 5000 km 2 or less when nitrogen concentration is considered alone; our estimates of required reductions in nitrogen load are in line with previous studies. The model that includes concentration and discharge as separate terms should be preferred for nutrient reduction scenario predictions, as discharge is not a manageable quantity. However, our results also suggest that it will be easier to measure statistically significant reductions in hypoxic area with a smaller per cent change in concentration. Thus, although our study suggests that large nitrogen reductions may be required to achieve an average hypoxic area of 5000 km 2 or less, smaller reductions will still have an observable impact on the extent of seasonal hypoxia.
Previous statistical models have also considered temporal trends in hypoxic area, either as a linear trend (Turner et al 2006), as a single step (Greene et al 2009), or as a series of steps (Turner et al 2008, Liu et al 2010. We do not believe that there are sufficient data to support (or reject) the hypothesis that multiple shifts have occurred in the system. We favor a single step in time, since this model may be used to investigate nitrogen load reduction strategies. However, as the reasons for the temporal trend in hypoxia are still uncertain, and the exact character of this trend is unknown, care must be taken when interpreting any predictions of hypoxic area.
A new finding presented in this letter is that wind has a significant effect on estimations of hypoxic areal extentexplaining a large percentage of the interannual variance in hypoxic area, independent of riverine processes. The relationship between eastward winds and decreasing hypoxic area may be explained through the effect of the wind on the Mississippi/Atchafalaya river plume system. In nonsummer, seasonal winds are downwelling, and tend to force the freshwater plume downcoast, toward the west. In summer, seasonal winds switch to upwelling favorable. These winds push the freshwater plumes upcoast, toward the east (Morey et al 2003(Morey et al , 2005. Stronger upwelling winds tend to push the plume further to the east, decreasing the area that is stratified by the plumes. Hetland and DiMarco (2008) refer to this as the stratification envelope, as the areal extent of hypoxia cannot exceed the area of the shelf that contains stratification great enough to support the formation of hypoxia. Thus, although the hypoxic potential, the maximum potential areal extent of hypoxia caused by a given nitrogen load, may be large, the actual observed area is constrained by the physical boundaries associated with the stratification envelope.
Our analysis identifies the hypoxic area observed in August 1988 as an outlier. Most previous statistical models of areal extent of hypoxia have included 1988 in their regressions. The consequences of excluding 1988 from the statistical analysis include significant shifts in the regression coefficients-a decrease in the correlation between hypoxic area and either nitrogen load or discharge, and a substantial increase in the correlation between hypoxic area and eastwest wind stress. No other year had such a profound effect on the regression statistics. Years with hurricanes that have been excluded from other studies show no discernible effect on the regression statistics. It is reasonable for process studies to exclude hurricane years as a method to reduce the noise variance, and to highlight the statistical relationships in the processes being investigated. However, we believe, as do Greene et al (2009), that these hurricane years must be included in any statistical model of northern Gulf of Mexico hypoxic area used for management purposes, because they are a natural and important part of the system variability.
The standard deviation in predicted area is greater than 3000 km 2 in all of the regression models identified as the best choices. These four models explain roughly half to threequarters of the observed variance in hypoxic area. Thus, although these models are able to explain a large portion of the interannual variability in hypoxia, as with many natural systems, there are still many aspects of the variation that we are not able to predict with simple models.