a large 
 hypoxic zone has formed nearly every year in the northern 
 gulf of mexico over at least the last few decades .
it is the largest 
 human - caused coastal hypoxic zone in the western hemisphere and one 
 of the largest worldwide . because of its negative 
 ecological impacts and potential fisheries impacts ,
the hypoxic zone has received attention from various stakeholders 
 and policy makers who have developed plans to reduce its average size .
the hypoxic zone is generally largest in summer due to the 
 seasonal 
 cycles of phytoplankton production and water column stratification , 
 which are the two primary factors leading to dissolved oxygen ( do ) 
 depletion .
both factors are related to the outflow from the 
 mississippi river basin , which typically peaks in march  may . 

nutrient loads from the basin stimulate phytoplankton production in 
 offshore waters , and much of the resulting organic matter eventually 
 settles along the louisiana  texas shelf .
aerobic bacteria decompose 
 the organic matter , consuming large quantities of do . at the same 
 time , freshwater outflows from the basin promote stratification through 
 a fresher surface layer , which warms in spring and summer , and overlays 
 a colder , saltier water layer , inhibiting reoxygenation of bottom 
 waters . the hypoxic extent , operationally 
 defined as the region ( area or 
 volume ) where
do concentrations are below 2 mg l , is often assessed using data from  shelfwide  sampling 
 cruises .
the cruises are considered  quasi - synoptic  
 because changes in weather conditions ( e.g. , tropical storms ) preceding 
 or during the cruises can disrupt typical seasonal patterns in hypoxia . 

cruises to document the occurrence of hypoxia , as well as related 
 physical and biological parameters , have been performed by the louisiana 
 universities marine consortium ( lumcon ) beginning in 1985 .
the estimated bottom - water hypoxic area has been determined by interpolating between sampling locations 
 and hand - contouring parallel to isobaths over a calibrated ( planimeter ) 
 grid .
these estimates have been used in multiple modeling studies 
 linking hypoxia to nutrient loads from the mississippi river basin 
 and other environmental factors and for setting policy goals to reduce the severity 
 of hypoxia . however , these estimates are generally conservative 
 in that the most inshore , most offshore , and most western extents 
 of hypoxia are not always captured because of logistical constraints .
also , these estimates do not quantify the uncertainty inherent in 
 making estimates from limited observations .
quantified uncertainties are useful for assessing the adequacy of 
 existing sampling programs and for improving models that link hypoxic 
 extent to environmental causes and effects .
the thickness and volume of hypoxia are also important , as 
 they relate to how hypoxia affects pelagic organisms .
in addition , hypoxic volume , rather than area , should be more closely 
 related to the total oxygen deficit of the system , potentially making 
 it a more useful metric for biogeochemical modeling studies .
while hypoxic volume measurements are available 
 for other systems , such as chesapeake bay , estimates for the gulf have only been available for a subset of 
 years , using an unpublished methodology .
the primary goal of this study is to improve our knowledge 
 of the 
 gulf s hypoxic extent over the 27-year study period ( 19852011 ) 
 by systematically estimating midsummer hypoxic area and volume .
instrument 
 bias adjustments  that account for the use of different sampling 
 instruments , capable of being lowered to within different proximities 
 of the sea floor , in different years .
our approach also uses trends between do , hypoxic fraction , depth , 
 and spatial position to account for consistent , large - scale spatial 
 patterns in do , thereby improving the explanatory power of the model .
our extent estimates are developed using a monte 
 carlo - type simulation approach , rather than traditional interpolation , 
 allowing for uncertainty quantification ( i.e. , confidence intervals ) .
these confidence intervals reflect the spatial stochasticity of 
 the system , uncertainties in trends among variables , and uncertainties 
 in instrument bias adjustments .
hypoxic volume is estimated using 
 a novel two - step approach where do is simulated first , and hypoxic 
 fraction ( i.e. , the fraction of the water column that is hypoxic ) 
 is simulated second using do as a trend variable .
we use 
 do data from lumcon midsummer sampling cruises conducted between 1985 
 and 2011 .
data for 19982008 are retrieved from the national 
 ocean data center , while data for other years are obtained directly 
 from lumcon .
sampling locations are geo - referenced 
 using the universal transverse mercator ( utm ) zone 15 projection , 
 and bathymetry is determined from a 3-arc - second digital elevation 
 model obtained from the national oceanic and atmospheric administration 
 ( noaa ) . during cruises ,
do is sampled 
 using one or two types of instruments ( figure 1 ) : a rosette - mounted do probe and a hand - held do probe , with the 
 latter capable of being lowered closer to the sea floor .
( hereafter , we refer to these instruments as the  rosette 
 sampler  and  handheld sampler  . )
while instrument 
 technology has changed over time , all instruments were calibrated 
 against winkler titrations , and postcruise corrections were made as 
 necessary . for sampling events where both instruments were used ,
( here , a  sampling event  refers to the depth profile 
 of data collected at a specific latitude / longitude and time . )
number of locations 
 sampled during the annual midsummer shelfwide 
 cruises using hand - held and rosette instruments .
we extract bottom water dissolved oxygen ( bwdo ) and minimum 
 dissolved 
 oxygen ( mindo ) concentrations from the do profiles .
bwdo and mindo 
 concentrations greater than 7 mg l ( less than 
 1% of total samples ) , are treated as 7 mg l , because 
 higher concentrations are outliers representing supersaturation conditions 
 that are not of interest in this study . for sampling events
where 
 hypoxia was present , we also extract the fraction of the water column 
 that is within the hypoxic bottom layer ( hereafter , bottom water hypoxic 
 fraction , or bwhf ) and , where other hypoxic layers are present , the 
 total hypoxic fraction ( thf ) of all hypoxic layers . bottom and upper 
 layers of hypoxia are observed at 49% and 12% of sampling events , 
 respectively .
here , we focus on the models for bwdo and bwhf because 
 bottom layer results are most comparable to previous studies . the models for mindo and thf are formulated in the same way and 
 yield comparable results , but for brevity are discussed only in supporting information ( si ) section s7 . in this study , adjustments are made to address biases that arise 
 from the use of different sampling instruments . for cruises when only 
 the rosette sampler was used ( figure 1 )
, one 
 would expect the hypoxic extent to be underestimated because the rosette 
 sampler does not reach as close to the sea floor as the hand - held 
 sampler .
we quantify this bias using data from sampling events where 
 both instruments were used together . for these cases ,
observed bwdo 
 and bwhf are calculated for both the synthesized profile and rosette - only 
 profile .
probabilistic relationships are then developed between the 
 synthesized results and the rosette - only results ( section s1 of the si ) . when performing 
 simulations ( described below ) , we adjust the rosette - only observations 
 by sampling from these relationships . a larger bias adjustment is 
 required for the first 38 sampling events in 1991 because the ship s 
 fathometer was not properly calibrated , causing the rosette sampler 
 to be lowered 1.5 m less than it would have been otherwise ( n.n .
rabalais , 
 cruise records ) . for these observations , the probabilistic relationships 
 between synthesized and rosette - only results are again determined 
 based on sampling events where both instruments were functioning properly 
 ( as described above ) but with the bottom 1.5 m of the rosette - only 
 profiles removed .
geostatistical modeling is performed on a 
 5  5 km grid of estimation points ( figure 2 ) , covering 
 342.5837.5 km utm easting and 3122.53292.5 km utm 
 northing , reflecting the general extent of sampling .
the estimation 
 grid is limited to depths between 3 and 80 m , which are typical of 
 the louisiana  texas shelf region where hypoxia occurs .
geostatistical 
 methods provide an effective means to model data that exhibit spatial 
 correlation .
the efficacy of these methods 
 has been demonstrated in previous environmental analyses of rainfall , 
 snow depth , soil phosphorus , and water quality indicators .
many of these 
 studies have focused on mapping through spatial interpolation , and 
 comparison studies have demonstrated that geostatistical interpolation 
 outperforms simpler interpolation methods , such as inverse distance 
 weighting . in a geostatistical model , 
 the response variable , z ,
is represented as a combination 
 of deterministic and stochastic components:1 this formulation is similar to linear regression , but includes 
 an additional term ,  , representing spatially correlated 
 stochasticity , in addition to the more commonly modeled uncorrelated 
 stochasticity , . as in linear regression , x is the portion of z that can be expressed 
 as a deterministic function of categorical and/or trend variables 
 ( x ) and their corresponding regression coefficients (  ) .
if the x term includes only cruise - specific 
 categorical variables ( vectors of zeros and ones that bin the data 
 by cruise ) , then the model can essentially be used to perform ordinary 
 kriging ( ok ) .
if x also includes trend variables , then 
 it can be used to perform universal kriging ( uk ) . in this study ,
the stochastic components (  and  ) are found to be well represented by the 
 commonly 
 used exponential covariance function with a nugget effect , defined as follows:2where q is 
 the covariance between two observations ( zi and zj ) separated by distance ( hi , j ) ;  and 
  are parameters 
 representing the variances of the stochasticity that is spatially 
 correlated and is not spatially correlated , respectively ; and r is a range parameter ( 3r is approximately 
 the distance at which observations are no longer spatially correlated ) . 

we allow for anisotropy in the covariance model by scaling hi , j using parameter  , 
 which represents the ratio of east  west to north  south 
 correlation ranges .
covariance function parameters are estimated using 
 restricted maximum likelihood and deterministic component 
 parameters (  ) are estimated using generalized 
 least - squares , as outlined in a previous study by obenour et al . that used geostatistical methods to explore biophysical 
 drivers of do depletion .
geostatistical models are developed 
 for both bwdo and bwhf to determine 
 bottom layer hypoxic area and volume , respectively .
the model for 
 bwdo uses all observations , while the model for bwhf uses only observations 
 from locations where bwdo is below the hypoxic threshold .
both models 
 use the uk formulation to take advantage of potential relationships 
 between the response variables and available trend variables .
represent relationships 
 between response and trend variables using regression coefficients 
 that are the same for all cruises .
for the bwdo model , potential constant 
 trends include linear and quadratic trends with depth ( depth and depth ) , easting ( easting and easting ) , and northing 
 ( northing and northing ) . on the basis of an examination 
 of model residuals
, it was noted that the overall trend between depth 
 and bwdo does not continue for depths greater than 40 m , and so depths 
 of greater than 40 m are treated as 40 m in the final bwdo model development . 

this approach is justified because model residuals are evenly distributed 
 across all depths when depths greater than 40 m are modified as described 
 here ( see figure s3 of the si ) . for bwhf ,
potential constant trends include 
 linear and quadratic trends with easting , northing , and bwdo .
second ,  cruise - specific trends  ( with cruise - specific 
 regression coefficients for easting ) represent relationships that 
 are specific to individual cruises . allowing these trends was motivated 
 by previous studies that indicate the east 
west 
 distribution of hypoxia is influenced by alongshore current velocity , 
 which can vary interannually in response to prevailing winds . when 
 selected
, these cruise - specific trends modify the constant trend . to prevent overparameterization of the model
, we use only trend 
 variables selected through a geostatistical adaptation of the bayesian 
 information criterion ( bic ) , wherein models with different subsets of trend variables are compared 
 based on their bic score , and the model with the lowest bic score 
 is optimal in terms of its explanatory power relative to its complexity . 
 because of the large number of variables considered in this study , 
 an initial optimal model is first selected among those formed from 
 constant trend variables only .
cruise - specific trend variables are 
 then individually tested relative to that model , and only cruise - specific 
 trends that improve the bic score are included in the final model .
the deterministic portion of the model also includes categorical 
 variables that bin data by cruise , reflecting the fact that the mean 
 bwdo varies from cruise to cruise . in uk
, the categorical variables 
 essentially allow for different  intercepts  ( as in 
 linear regression ) , that shift the trend up or down to best fit the 
 observations from each cruise . in ok
the covariance 
 model , selected deterministic variables , and categorical 
 variables are used to determine a set of geostatistical weights ,  , that are applied to observations when 
 performing geostatistical interpolation and simulation . for each cruise ,
y , a unique set of weights , y , are determined by solving 
 a system of linear equations:3where qoo is an n  n covariance matrix for the n observation 
 locations 
 ( all cruises ) , with elements determined from eq 2 . because we did not assume correlation among stochasticity from 
 different cruises , intercruise covariances are assigned a value of 
 zero .
similarly , qoe , y is an n  m covariance 
 matrix of n observation locations and m estimation locations , and the rows of qoe , y that correspond to observations 
 from cruises other than cruise y are assigned a value 
 of zero .
the matrix xo is n  p and includes the p deterministic variables ( trend and categorical ) for the observation 
 locations , and the matrix xe ( m  p ) includes the same 
 variables for the estimation locations .
the trend variables are normalized 
 to a mean of zero and variance of one .
in eq 3 , y is an n  m matrix of cruise - specific weights that reflect both the 
 spatial correlation structure and the deterministic trends .
also , gy is a p  m matrix of lagrange multipliers that can 
 be used with y to determine location - specific estimation uncertainties . using 
 the weights , along with the observations zo ( an n  1 vector )
, one can 
 develop estimates of the response variable across the estimation grid:4where ze , y0 is an m  
 1 vector of interpolated bwdo or bwhf values for cruise y. although we use a uk model formulation , hypoxic area and 
 volume 
 are not determined by kriging ( i.e. , spatial interpolation ) , but instead 
 by developing  spatially consistent monte 
 carlo simulations  which are often 
 referred to as conditional realizations ( crs ) .
while both kriging 
 and crs provide equivalent information for individual estimation locations , 
 crs are necessary to estimate spatially aggregated quantities ( e.g. , 
 area and volume ) probabilistically .
a cr is performed by first creating 
 an  unconditional realization  and then  conditioning  
 it to the observed data and deterministic trend .
unconditional realizations ( eq 5 ) 
 include simulated values at the estimation locations ( ze , yu ) as well as at the observation locations ( zo , yu):5note that the vector zo , yu includes simulated values corresponding to 
 the observations from 
 all cruises , but only observations from cruise y have 
 their stochasticity correlated with that of the estimation locations . 
 here , qee is the m  m covariance matrix between estimation 
 locations , and u is an ( m + n )  1 vector of random independent samples from the 
 standard normal distribution .
the operator c ( ) returns 
 the triangular matrix resulting from cholesky decomposition of the 
 subject matrix .
the unconditional realizations are then conditioned 
 to the observed 
 data and deterministic trends through:6here , zo is the n  1 vector 
 of the observed 
 values , and ze , yc is the resulting cruise - specific 
 cr , an m  1 vector of values corresponding 
 to the estimation locations .
then , bwhf is simulated 
 over those locations where the simulated bwdo is below the hypoxic 
 threshold . at hypoxic estimation locations ,
the simulated bwdo values 
 are used as a trend variable in the bwhf model .
simulated values are 
 limited to within realistic ranges ( determined to be 07 mg 
 l for bwdo and 0.0010.8 for bwhf , based 
 on a review of the observed data ) . without these constraints ,
the 
 estimated hypoxic volume would be 6% smaller , on average , as a result 
 of ( unrealistic ) negative hypoxic thicknesses .
the two - step 
 cr process is repeated 1000 times for each cruise , 
 resulting in 1000 realizations of both bwdo and bwhf . for each realization 
 of bwdo ,
the hypoxic area is calculated as the number of estimation 
 locations simulated to be hypoxic multiplied by the grid cell area 
 ( 25 km ) ; and for each realization of bwhf , the hypoxic 
 volume is calculated as the vector of simulated bwhf values multiplied 
 by the corresponding vector of water column depths , all multiplied 
 by the grid cell area . from this ensemble of results ,
we determine 
 the mean and 95% confidence intervals for the hypoxic area and volume 
 for each cruise .
the bwdo 
 and bwhf models include several parameters that characterize 
 the deterministic and stochastic model components ( eq 1 ) .
regression coefficients for the bic - selected trend variables , 
 explaining a portion of the spatial variability in observed bwdo and 
 bwhf , are provided in table 1 .
the standard 
 errors of these coefficients are low ( i.e. , p < 
 0.05 ) , suggesting these trends are statistically significant .
the 
 coefficients for the cruise - specific categorical variables , accounting 
 for year - to - year variability in the responses , are included in section s3 of the si .
overall , the deterministic model components explain 28% and 32% 
 of the total ( spatial plus interannual ) variability in bwdo and bwhf , 
 respectively , while the stochastic components of the models explain 
 the remainder of the spatial variability .
if trend variables are omitted 
 ( i.e. , the ok formulation ) , then the deterministic components only 
 account for interannual variability , and only 12% and 11% of the total 
 variability in bwdo and bwhf is explained , respectively . a check of 
 the linearity assumption implicit in these deterministic relationships 
 is included in section s2 of the si .
c.s.e 
 = cruise specific trend 
 for easting , n.s.=not selected , n.a.=not available .
the trend between bwdo and easting is quadratic with a minimum at 
 794 km easting ( utm coordinate , figure 2 ) , 
 between the mississippi and atchafalaya river outfalls , which is reasonable 
 given that these rivers provide the freshwater flows and nutrients 
 that are important for hypoxia formation . in 1998 , the only year for 
 which a cruise - specific east 
west trend with bwdo is selected , 
 the minimum is shifted to 1111 km easting , outside of the study area , 
 indicating that bwdo concentrations decrease monotonically with easting 
 within the study area .
the unique spatial distribution in 1998 has 
 been noted previously , and is generally attributed to unusually persistent 
 eastward currents . a similar bwdo pattern 
 has been noted for 2009 , but it was not 
 sufficiently strong to result in a cruise - specific trend in this analysis . 

the trend between bwdo and northing is linear , and suggests that bwdo 
 concentrations are higher to the south .
the trend between bwdo and 
 bathymetry ( depth ) is quadratic with a minimum at 22 m. the overall 
 deterministic trends in bwdo ( figure s7 of the si ) compare well with spatial 
 patterns in hypoxic frequency , as mapped by walker and rabalais . because bwhf is modeled as a function 
 of bwdo , it effectively inherits 
 the trends from the bwdo model .
the trend with bwdo is linear and 
 suggests that bwhf is larger where bwdo concentrations are lower . 

the model for bwhf also includes a quadratic trend with easting , suggesting 
 a maximum at 701 km , about 30 km east of the atchafalaya outfall location 
 ( although the overall trend in bwhf is also affected by the trends 
 in bwdo ) .
. for bwdo ,  and  , commonly referred to as the nugget and partial 
 sill , were 0.39 and 2.33 ( mg
since  is smaller 
 than  , the majority 
 of the stochasticity in the model is spatially correlated .
the approximate 
 range of spatial correlation ( 3r , per eq 2 ) is 94 km in the east  west direction and 
 53 km in the north  south direction .
the greater correlation 
 distance in the east  west direction was expected due to the 
 dominant east  west current pattern . for bwhf ,
the spatial correlation of the stochasticity is somewhat 
 weaker with  and  similar 
 in magnitude , having values of 0.009 and 0.011 ( unit - less ) , respectively . 

also , the correlation range was 63 km in all directions ( anisotropy 
 was negligible ) .
overall , spatially correlated stochasticity accounts 
 for the greatest portion of the variability in both bwdo and bwhf . 

, the spatial 
 correlation of the stochasticity is consistent with the effects of 
 varying coastal current patterns , influencing the distribution of 
 hypoxia over the spatial scales described above . in general , the correlation 
 ranges are considerably longer than the typical distances between 
 sampling locations ( figure 2 ) , especially on 
 the eastern shelf , suggesting that the sampling network is adequate 
 for resolving the spatially correlated stochasticity of the system . using cr , we determined the mean and 95% confidence intervals for 
 hypoxic area and volume of each cruise ( figure 3 ) .
the largest estimated hypoxic area was for 1996 , but it was not 
 significantly different from 1991 , 1993 , 1995 , 1997 , 1999 , 2001 , 2002 , 
 2007 , or 2008 ( p > 0.05 ) , given the uncertainty 
 in 
 the area estimates .
the largest hypoxic volume estimate was for 2008 , 
 but it was not significantly different from 1993 , 1996 , 1999 , 2004 , 
 or 2007 ( p > 0.05 ) .
bottom 
 layer hypoxic extent estimates with 95% confidence intervals 
 by year ; estimates prior to making adjustments for instrument bias 
 as triangles ; previous lumcon area estimates as open squares ; revised 
 lumcon area estimates as solid squares .
two sets of hypoxic area estimates , determined by lumcon , 
 are also 
 included in figure 3 for comparison .
the first 
 set are the original lumcon estimates , determined by hand contouring , 
 as described above .
the second set are revised estimates 
 developed as part of this study ; also by hand contouring , but using 
 the updated bwdo values for this study ( which are superior because 
 they reflect postcruise do calibrations ) . in most cases , the original 
 lumcon estimates ( figure 3 , open squares )
overlap 
 the revised estimates ( solid squares ) , and the majority of the lumcon 
 estimates fall within the 95% confidence intervals determined by this 
 study .
however , for the first third of the study period ( 19851993 ) , 
 the new estimates are consistently higher than the lumcon estimates .
lumcon did not make an estimate for 1989 , because only a portion 
 of the shelf was sampled that year .
our method does allow for extent 
 estimates for 1989 , and the wide confidence intervals ( figure 3 ) reflect the relatively large uncertainties . even 
 though data for 1989 were limited , the overall spatial trends ( based 
 on the data from all years ) help to constrain the variability in the 
 crs to within realistic ranges , across the study area .
we note that 
 the 1989 estimate presented here is consistent with previous estimates 
 developed from nutrient loading models .
there 
 was considerable interannual variability in the spatial distribution 
 of hypoxia and in the thickness of the hypoxic bottom layer ( figure 4 ) .
years with similar hypoxic areas may have very 
 different average hypoxic thicknesses and thus volumes .
for example , 
 2002 and 2008 have similar hypoxic areas , but they are very different 
 in terms of hypoxic thickness , such that 2008 has approximately twice 
 the hypoxic volume of 2002 .
the results indicate that the average 
 bottom layer hypoxic thickness for the 27-year study period was 3.9 
 m , with the thickest average layers of approximately 6.2 and 6.3 m 
 occurring in 2008 and 2009 , respectively . over the 27-year study period
interestingly , hypoxic area and thickness are also somewhat 
 positively correlated over most years ( r = 0.37 , not including 1998 and 2009 , which were subject to unusually 
 strong eastward currents ) , such that volume 
 increases exponentially relative to area . tabulated results and additional 
 maps showing estimated bwdo , estimated hypoxic thickness , probability 
 of hypoxia , and example crs for the entire study period are included 
 in sections s4 and s5 of the si .
example maps of estimated bottom layer hypoxic 
 thickness ( median 
 values from crs ) , 20012008 ; observation locations shown as 
 white dots .
the geostatistical modeling results can be used to assess temporal 
 trends in hypoxic zone size for the 27-year study period .
hypoxic 
 volume increased by an average of 2.3% per year as a linear trend 
 with time ( 1.4 km yr , percent increase 
 determined by dividing by the 27-year mean hypoxic extent ) , but this 
 trend was not significant ( p = 0.12 ) .
hypoxic area 
 increased at a lesser rate of 0.9% per year ( 140 km yr ) and was also not significant ( p = 0.42 ) .
however , the relatively large increase in volume relative 
 to area reflects a significant increasing trend in hypoxic layer thickness 
 ( 1.8% per year , 0.069 m yr , p = 0.05 ) .
note that trend significance is affected by the uncertainty 
 in the geostatistical estimates ; trend significance is determined 
 using a monte carlo approach where temporal trend coefficients are 
 developed for each of the 1000 sets of crs , and an overall probability 
 distribution is developed by sampling from the uncertainty in these 
 trend coefficients ( note p - values are based on a 
 two - sided test ) . without accounting for uncertainty ,
the volume , area , 
 and thickness trends would have p - values of 0.09 , 
 0.38 , and 0.01 , respectively .
the new hypoxic area estimates 
 can be compared to the previous 
 hypoxic area estimates developed by lumcon . for 
 the 27-year study period ,
the previous area estimates increased at 
 a highly significant rate of 2.6% per year ( 360 km yr , p = 0.01 ) , substantially greater 
 than the 0.9% rate for the new estimates .
this dissimilarity is due 
 primarily to differences in estimates for the earlier years of the 
 study period ( figure 3 ) . for 19851987 , 

the new estimates are consistently higher than the previous estimates 
 primarily because the geostatistical methodology accounts for the 
 possibility of hypoxia occurring outside the envelopes of the cruises , 
 which were relatively small in these years . for 1990 , 1991 , and 1993 , 

the new estimates are also higher because of the instrument bias adjustments 
 developed in this study . for 19851993 ( except 1989 ,
for which 
 there is no previous estimate ) , the mean geostatistically determined 
 hypoxic area is 39% ( 3650 km ) greater than that of the 
 previous estimates . for 19942011 ,
the new area estimates 
 are in general agreement 
 with the previous area estimates . over this time period ,
the two sets 
 of estimates are highly correlated ( r= 0.88 ) , the means of the two data sets are of negligible difference , 
 and only three of the previous estimates ( 1996 , 2003 , and 2010 ) fall 
 outside the 95% confidence intervals of the new estimates ( figure 3 ) .
this suggests that when the biasing issues noted 
 above are avoided , there is approximate agreement between the geostatistical 
 and hand - contouring estimates .
the new extent estimates have 
 implications for our understanding 
 of how hypoxia is changing over long temporal scales .
multiple studies 
 using nutrient loads to predict the previous hypoxic area estimates 
 have suggested the gulf is becoming increasingly susceptible to hypoxia , 
 based on increasing hypoxic area relative to nutrient loading .
( nutrient loading increased greatly in the 1970s but remained relatively 
 stable throughout the study period . )
increases 
 in hypoxic susceptibility have also been suggested in other studies 
 of the gulf and other coastal systems .
the new hypoxic area 
 estimates , however , exhibit relatively little increase over time , 
 especially when compared to the previous estimates , potentially suggesting 
 less system change during the study period than previously thought . 

but , the new hypoxic thickness and volume estimates do increase to 
 a greater degree , potentially suggesting a more vertically oriented 
 increase in hypoxic extent .
future studies could focus on recalibrating 
 existing nutrient - loading models to these new hypoxic extent estimates 
 to develop a refined understanding of how gulf hypoxia may be changing 
 over time .
uncertainties in the new hypoxic extent estimates , 
 represented 
 by the 95% confidence intervals ( figure 3 ) , 
 reflect the limited spatial scope and resolution of shelfwide cruise 
 sampling .
uncertainties are generally greatest in the earlier years 
 when cruises were smaller and did not always use instruments that 
 reached the sea floor . from 1985 to 1993 , the mean relative standard 
 error for hypoxic area was 23% , but it decreased to 11% for 19942011 . 
 uncertainties also appear to be larger in years with relatively severe 
 hypoxia ( i.e. , lower average bwdo and larger hypoxic area ) , likely 
 because more of the estimation grid is subject to the possibility 
 of hypoxia ( in years with higher average bwdo , most of the estimation 
 grid is determined to be well above the hypoxic threshold , such that 
 simulated values rarely fall below the hypoxic threshold ) . in the 
 future
, the modeling framework presented here could be used to evaluate 
 different sampling designs based on how well they constrain estimate 
 uncertainty .
the modeling approach could also be used to compare hypoxic 
 extent results and associated uncertainties using hypoxic threshold 
 choices other than 2 mg l. regardless of 
 the precision of the estimates , it is important to 
 remember that they represent conditions at only certain points in 
 time ( i.e. , during the shelfwide cruises ) , and that hypoxic extent 
 can vary substantially throughout the summer due to changes in organic 
 matter production , wind - driven mixing events , and fluctuating current 
 patterns . as a result , 
 these estimates do not necessarily reflect hypoxic conditions over 
 the entire summer ( although back - to - back cruises have demonstrated 
 fairly consistent hypoxic areas under stable weather conditions ) .
more detailed monitoring and biophysical modeling 
 are needed to better understand the short- and long - term dynamics 
 of hypoxia formation , and to mechanistically interpret the temporal 
 variability in the extent estimates presented here .
the primary 
 feature of the geostatistical approach , when compared 
 to more traditional , interpolation - based approaches , is the use of 
 simulations ( i.e. , crs ) .
the most obvious 
 benefit is the ability to quantify the uncertainty in extent estimates . 

a second benefit is that crs provide more realistic extent estimates 
 than can be derived from kriging ( or largely equivalent methods such 
 as gauss - markov smoothing ) alone . in 
 this study ,
hypoxic area estimates derived from kriged maps ( uk interpolation , figure s14 of the si ) were substantially lower than both our cr estimates and the lumcon 
 estimates .
this is because the kriged maps ( figure 
 s9 of the si ) tend to characterize 
 large portions of the estimation grid as slightly above the hypoxic 
 threshold .
however , due to the stochasticity of the system , these 
 locations still have some probability of being hypoxic .
the cr approach , 
 which samples from the uncertainty in the system , accounts for this 
 possibility , and thus results in larger hypoxic area estimates .
this 
 is consistent with chiles and delfiner , who argue that cr , rather than kriging , is the more appropriate 
 approach for determining spatially aggregated quantities .
a third 
 advantage of the cr approach is that it provides a framework for performing 
 instrument bias adjustments probabilistically , such that adjustment 
 uncertainty is propagated to the hypoxic extent estimates .
finally , 
 cr was fundamental to probabilistically determining hypoxic volume 
 ( in addition to area ) .
the two - step cr approach , developed here , can 
 be compared to other volume estimation methods , such as multilayer 
 kriging and three - dimensional cr , which have been applied in chesapeake bay .
a 
 benefit of our approach is that it allows uncertainty quantification 
 of hypoxic volume within a relatively simple ( two - dimensional ) geostatistical 
 framework .
this study also demonstrates the benefits of including 
 trend variables 
 within the geostatistical model ( i.e. , the uk formulation ) because 
 the deterministic trends help reduce model uncertainty and result 
 in more realistic extent estimates ( section s6 of the si ) .
the advantages of uk have 
 been demonstrated previously for lake erie and chesapeake bay hypoxia , but are perhaps even more salient for an open system such as the 
 gulf shelf . when trends in do are not modeled , it has been necessary 
 to limit the estimation grid around the bounds of the sampling cruise , such that the size of the cruise can potentially bias the inferred 
 hypoxic area . by including trend variables that explain the large - scale 
 spatial patterns in bwdo and bwhf ,
it is possible to develop realistic 
 crs of do across the entire study area , so the same estimation grid 
 can be used for all cruises .
finally , the results confirm that 
 the hypoxic area on the louisiana - texas 
 shelf greatly exceeds the hypoxia task force goal of 5000 km as a five - year running average .
the most recent 
 five - year period from our study , 20072011 , has a mean hypoxic 
 area of 16  600 km with a 95% confidence range of 
 15  10018  000 km .
clearly , additional 
 management measures are required if the hypoxic extent is to be reduced 
 to comply with the task force goal .
 robust estimates of hypoxic extent 
 ( both area and volume ) are important 
 for assessing the impacts of low dissolved oxygen on aquatic ecosystems 
 at large spatial scales . 
 such estimates are also important for calibrating 
 models linking hypoxia to causal factors , such as nutrient loading 
 and stratification , and for informing management decisions . in this 
 study 
 , we develop a rigorous geostatistical modeling framework to 
 estimate the hypoxic extent in the northern gulf of mexico from data 
 collected during midsummer , quasi - synoptic monitoring cruises ( 19852011 ) . 
 instead of a traditional interpolation - based approach , 
 we use a simulation - based 
 approach that yields more robust extent estimates and quantified uncertainty . 
 
 the modeling framework also makes use of covariate information ( i.e. , 
 trend variables such as depth and spatial position ) , to reduce estimation 
 uncertainty . 
 furthermore , adjustments are made to account for observational 
 bias resulting from the use of different sampling instruments in different 
 years . 
 our results suggest an increasing trend in hypoxic layer thickness 
 ( p = 0.05 ) from 1985 to 2011 , but less than significant 
 increases in volume ( p = 0.12 ) and area ( p = 0.42 ) . the uncertainties in the extent estimates vary 
 with sampling network coverage and instrument type , and generally 
 decrease over the study period . 
Introduction
Methods
Results
Discussion