Temporally resolved coastal hypoxia forecasting and uncertainty assessment via Bayesian mechanistic modeling

Abstract. Low bottom water dissolved oxygen conditions (hypoxia) occur almost every summer in the northern Gulf of Mexico due to a combination of nutrient loadings and water column stratification. Several statistical and mechanistic models have been used to forecast the midsummer hypoxic area, based on spring nitrogen loading from major rivers. However, sub-seasonal forecasts are needed to fully characterize the dynamics of hypoxia over the summer season, which is important for informing fisheries and ecosystem management. Here, we present an approach to forecasting hypoxic conditions at a daily resolution through Bayesian mechanistic modeling that allows for rigorous uncertainty quantification. Within this framework, we develop and test different representations and projections of hydrometeorological model inputs. We find that May precipitation over the Mississippi River basin is a key predictor of summer discharge and loading that substantially improves forecast performance. Accounting for spring wind conditions also improves
forecast performance, though to a lesser extent. The proposed approach generates forecasts for two different sections of the Louisiana–Texas shelf (east and west), and it explains about 50 % of the variability in the total hypoxic area when tested against historical observations (1985–2016). Results also show how forecast uncertainties build over the summer season, with longer lead times from the nominal forecast release date of 1 June, due to increasing stochasticity in riverine and meteorological inputs. Consequently, the portion of overall forecast variance associated with uncertainties in data inputs increases from 26 % to 41 % from June–July to August–September, respectively. Overall, the study demonstrates a unique approach to assessing and reducing uncertainties in temporally resolved hypoxia forecasting.



Contents
), is based on a steady-state solution to mass balance differential equations presented in Obenour et al. (2015).Using time-varying inputs, DMO20 predicts daily BWDO (Ob) concentration (mg/L) from 1 June to 30 September of each year: where ka is the reaeration rate (m/d), Dw is the water column oxygen demand (WCOD), and Ds is the sediment oxygen demand (SOD) at Of, a reference oxygen concentration set to 3 mg/L.
The net WCOD (g/m 2 /d) of each lower compartment is represented as: Here, J is the downward carbon flux (g/m 2 /d), γ is the mass ratio of oxygen demand to organic carbon set to 3.5, λ is the ratio of organic carbon to nitrogen set to 5.68, A is the area of the shelf section (Gm 2 ), ω is an oxygen demand adjustment factor (accounting for photosynthesis, off-shelf losses, etc.), and  is the effective settling velocity (m/d), which incorporates both the production and sinking of organic matter.The variables Q and L represent the near-term flows (Gm 3 /d) and N loads (Gg/d) entering the surface-layer model compartments, respectively.Subscripts r, u, and g denote the origin of these fluxes: Mississippi and Atchafalaya Rivers, upstream (i.e., eastern) shelf section, and the greater Gulf of Mexico.Qg is approximated as a dilution factor (3.2, derived from surface salinity data) multiplied by mean Mississippi River discharge (1.6 Gm 3 /d).
The reaeration rate for each section is determined as a function of wind stress (representing shear-induced turbulence) and freshwater flow (representing buoyancy): where U is the 14-day weighted mean wind speed for the shelf section (m/s), Qs is the river discharge entering the section (Gm 3 /d), and β0 and β1 are calibration parameters.
Partitioning of riverine inputs is computed through: where FW is the fraction of abovementioned flows and loads transported westward over the shelf, ve is the mean eastward wind velocity (m/s), and βe is a calibration parameter.The 0.5 indicates that, in absence of wind, inputs from both rivers would equally partition westward and eastward.
SOD is represented as: where L (Mg/mo) is the combined nutrient loading from the Mississippi and Atchafalaya Rivers, averaged November-March.
We normalize these pre-spring loads relative to their long-term average ( � ) for the study period.SOD is temperature dependent and is based on the Arrhenius model with θ = 1.07.Rates are corrected when temperatures deviate from  � , the summertime average.Here, T is the monthly mean temperature.
A quadratic polynomial function g is used to transform modelled bottom water dissolved oxygen (mg/L) into hypoxic area (km 2 ).For each section of the shelf, west and east, g is: For Eq. S6 the R 2 = 0.98 and the residual standard deviation is 706 km 2 , while for Eq.S7 the R 2 = 0.99 and the residual standard deviation is 216 km 2 .

S2 Bias adjustment for model predictions in June
Preliminary analysis indicated that hindcasted BWDO was somewhat lower than observations for the west section of the shelf in June.This bias remained after conversion of BWDO to HA.A linear regression with zero intercept and a sequence of numbers 29 to 0 representing period from June 1 to June 30 as the predictor was used to estimate a bias adjustment, defined as difference between average observed and hindcasted BWDO divided by hindcasted BWDO.The resulting regression indicates a gradual decline in bias from the beginning of June (Fig. S2.1).The bias adjustment increases the R 2 of relationship between observed and hindcasted BWDO from -0.15 to 0.36 (Fig. S2.2), and between observed and hindcasted HA from -0.14 to 0.45 (Fig. S2.3).

Figure S2. 1 :
Figure S2.1:Bias adjustment factor vs day number (before 1 July) for the west section with red line showing the regression fit with slope mean and standard error of 0.007 and 0.001, respectively.The adjusted R 2 of this regression is 0.20.

Figure S2. 2 :
Figure S2.2:Month by month comparison of observed with hindcasted (black) and bias-adjusted hindcasted (red) averaged BWDO in the west and east sections.Diagonal lines represent perfect prediction.

Figure S2. 3 :
Figure S2.3:Month by month comparison of observed with hindcasted (black) and bias-adjusted hindcasted (red) averaged HA in the west and east sections.Diagonal lines represent perfect prediction.

Figure S4. 1 :
Figure S4.1:Daily hindcasted HA and observed HA versus pseudo-forecasted HA for the west and east sections.Diagonal line represents perfect prediction.

Figure S4. 2 :
Figure S4.2:Averaged daily variance of total HA due to different sources of uncertainty.In this case, the "residual error" variance includes transformation and bias uncertainties, in addition to the DMO20 residuals.Note that the relative magnitudes of the variance components are somewhat different from the magnitudes of the IQR components (e.g., Fig. 5) because variance has squared units.

Figure S5. 1 :
Figure S5.1:Daily pseudo-forecasts of HA for the west and east sections in 1985 (top) and 1986 (bottom), including 95% IQR of the predictive distribution, distinguishing between i) parameter, ii) hydrometeorological inputs, iii) mechanistic model error, and iv) regressions related to transformation of BWDO to HA and bias adjustment uncertainties (shades of gray from lightest to darkest).Yellow dashed line is hindcasted estimate, black dashed line is the 32-year average hindcast, orange points and error bars represent the mean and associated 95% confidence interval of the (geostatistically estimated) hypoxia observations.