Hydrodynamic and Biochemical Impacts on the Development of 1 Hypoxia in the Louisiana–Texas Shelf Part II: Statistical Modeling 2 and Hypoxia Prediction

. In this study, a novel ensemble regression model was developed for hypoxic area (HA) forecast in the Louisiana– 10 Texas (LaTex) Shelf. The ensemble model combines a zero-inflated Poisson generalized linear model (GLM) and a quasi- 11 Poisson generalized additive model (GAM) and considers predictors with hydrodynamic and biochemical features. Both 12 models were trained and calibrated using the daily hindcast (2007–2020) by a three-dimensional coupled hydrodynamic– 13 biogeochemical model embedded in the Reginal Ocean Modeling System (ROMS). A promising HA forecast is provided by 14 the ensemble model with a low RMSE (3,204 km 2 ), a high R 2 (0.8005), and a precise performance in capturing hypoxic area 15 peaks in the summers. To test its robustness, the model was further applied to a global forecast model and produces HA 16 prediction from 2019 to 2020 with the adjusted predictors from the HYbrid Coordinate Ocean Model (HYCOM). Predicted 17 HA shows a high agreement with the ROMS hindcast time series (RMSE=4,571 km 2 , R 2 =0.8178). Our model can also predict 18 the magnitude and onsets of summer HA peaks in both 2019 and 2020 with high accuracy. To the best of our knowledge, this 19 ensemble model is by far the first one providing fast and accurate daily HA predictions for the LaTex Shelf while considering 20 both hydrodynamic and biochemical effects. This study demonstrates that it is feasible to perform regional ocean HA prediction 21 using global ocean forecast.


Introduction 23
The Louisiana-Texas (LaTex) Shelf has become a center of hypoxia (bottom dissolved oxygen, DO<2 mg L-1) study since where is the rate of surface heat input, is the volume expansion coefficient, is water specific heat capacity, is 100 coefficient of wind mixing, 1 is drag coefficient, 1 is humid air density near the sea surface, and is the wind speed near 101 sea surface. The first term on the right-hand side of Eq. (2) represents the rate of change of water stratification due to surface 102 heating, while the second term is the rate of working by wind stress contributing negatively to water stratification. Therefore, 103 the heat-induced change of PEA is proportional to the product of heat input and water depth, which is, 104 The total net heat flux, a sum of net shortwave and net longwave radiation flux, is derived from the National Centers for 108 Environmental Prediction Climate Forecast System (NCEP) Reanalysis (CFSR) 6-hourly products (Saha et al., 2010;2011) in 109 this study. The term (Qh) is added to the candidate list of predictors and is denoted as PEA 4567 (heat-induced PEA changes) 110 for simplification. 111 112 Daily variability of term ( 1 1 2 ) is dominated by that of 2 , since the 1 fluctuates much less than the 2 in a daily 113 scale ( Figure A1). We obtained the 1 according to (Picard et al., 2008) : 114 where represents the absolute air pressure, & (=28.96546 g mol -1 ) is the molar mass of dry air, = (=18.01528 g mol -1 ) is 118 the molar mass of water vapor, indicates compressibility, (=8.314472 J mol -1 K -1 ) is the molar gas constant, is 119 thermodynamic temperature, = is the mole fraction of water vapor. We assumed that air parcels at the sea surface are ideal 120 gases ( = 1) and are always saturated with water vapor. Thus, = is a function of absolute air pressure ( ) and saturation 121 vapor pressure of water ( >1, ) and can be calculated as follows: 122 The 1 is then estimated using sea surface air pressure and air temperature 2 meters above the sea surface provided by NCEP 135 CFSR 6-hourly products. Correlation of daily 1 2 and 2 (provided by NCEP CFSR 6-hourly products) is significantly 136 high as 0.9989 (p<0.001, Figure A1) emphasizing the importance of term 2 in controlling the daily variability of wind-137 induced PEA changes over the shelf. We, thus, approximated the relationship as: 138 Information System (NWIS). Due to lateral transports and vertical settling of particulate organic matter, a leading period should 152 be introduced to the time series of riverine nutrient loads. The optimal length of leading days is obtained by examining the 153 highest linear correlation of regionally averaged ROMS-hindcast SOC and SOCalt following Eq. (9) and is calculated as 19 154 days ( Figure A2a). The exponential term in Eq. (9) estimates the temperature-dependent decomposition rate of organic matter. 155 A significant correlation coefficient between daily SOCalt and ROMS-hindcast SOC is found as 0.8157 (p<0.001, Figure A2). 156 157 = Mississippi River inorganic nitrogen loads (led by 19 days) • e F.FHI2< / , where J indicates bottom water temperature (in °C ). Along with SOCalt, the temperature-dependent decomposition rate 160 F.FHI2•< / is also considered as a candidate predictor in statistical models and is denoted as DCPTemp for simplification. 161

Data pre-processes 178
We applied the spatially averaged daily ROMS-derived predictors over the LaTex Shelf, then applied the min-max 179 normalization (Eq. (10)) to the one-dimensional time series. Predictive models can be beneficial from the min-max 180 normalization when applying to a new dataset since the method guarantees that the normalized predictors from different 181 datasets range from 0 to 1 as the minimum and maximum values are prescribed. Note that the response is not normalized. 182  method is applied to each model yielding 10 root-mean-square errors (RMSEs) and 1 corresponding mean. The candidate 206 predictors of PEA and SOCalt are forced into each subset. Thus, the number of fitted models with a subset size of k is 207 2 ≤ ≤ 6 (the total number of candidate predictors is 6). The optimal subset of this size is 208 found as the one with the lowest mean CV RMSE among these models. The best subset is then obtained by comparing mean 209 CV RMSEs of the optimal subsets of different sizes.

Regular GLMs and zero-inflated GLMs 222
The response variable can be treated as count data. Regular Poisson (function glm in R package "stats" version 3.6.2), quasi-223 Poisson (function glm in R package "stats" version 3.6.2), and negative binomial (function glm.nb in R package "MASS"

Model interpretation for GLMzip3 284
We applied the complete ROMS training set to the model construction of GLMzip3 and found the coefficients for PEA, 285 SOCalt, and DCPTemp (Table 2) are all significantly positive (p<0.001) in the count model, while coefficients for these 286 predictors are significantly negative (p<0.001) in the zero-excess model. The count model simulates the HA while the zero-287 excess model estimates the probability of HA being zero. Higher PEA is consistent with stronger water stratification, while 288 higher SOCalt and DCPTemp are both corresponding to higher sediment oxygen consumption. Therefore, there is no surprise 289 that higher PEA, SOCalt, and DCOTemp are related to greater HA and higher hypoxia occurrence or lower probability of HA 290 being zero. Results indicate that the GLMzip3 essentially builds up reasonable relationships between the response and 291 predictors variables with a high agreement with physical and biochemical mechanisms. 292  sensitive to the predictors in the low-value ranges but becomes nearly stable in the medium-and high-value ranges of 315 predictors. It implies that bottom hypoxia develops rapidly in early summer when water stratification and sediment oxygen 316 demand start to increase. The bottom hypoxic water further extends with a much lower expansion speed as the stratification 317 and SOC further intensify. Nevertheless, the smooth function of PEA is slightly greater also with a more acute slope than those 318 found for SOCalt and DCPTemp in the medium-and high-value regimes of the predictors. It indicates that the HA variability is 319 more related to the hydrodynamic changes in the shelf than the biochemical effects. The result is consistent with the findings 320 by previous studies of the shelf hypoxia (Yu et al., 2015;Mattern et al., 2013) emphasizing that the physical impacts are 321 stronger than the biological impacts on HA estimates. A short conclusion is made that the GAMqsp3 model provides reasonable 322 interpretations on the hypoxic area mechanisms. 323

328
The prediction performance of GAMqsp3 is estimated using the Bagging ensemble method (Figure 4b). The RMSE and R 2 329 between the Bagging mean and ROMS-hindcast HA is 3,134 km 2 and 0.8093, respectively. They are 13 % lower and 8 % 330 higher than the corresponding statistics found for the GLMzip3, respectively, suggesting that GAMqsp3 outcompetes 331 GLMzip3 in terms of overall performance. However, GAMqsp3 tends to produce underestimated predictions at HA peaks 332 (like peaks around the 310th and 920th samples) some of which are overestimated by the GLMzip3. Therefore, instead of 333 determining the best model out of the two, ensemble HA predictions blending efforts of both GLMzip3 and GAMqsp3 are 334 carried out with an expectation to improve model performance in the peak forecast. We assumed that the contributions of 335 GLMzip3 and GAMqsp3 are equally weighted since there is no clue showing the apparent superiority of either model in HA 336 peak predictions. We thus averaged the predicted HA by GLMzip3 and GAMqsp3 and calculated the 95 % PIs given the 337 Bagging results of these models ( Figure 4c). As expected, the overall performance of the ensemble forecast is somewhere 338 between the performance of GLMzip3 and GAMqsp3 with an RMSE of 3,204 km 2 and an R 2 of 0.8005. However, some HA 339 peak events (like peaks around the 310th and 920th samples) which are overestimated by GLMzip3 but are underestimated by 340

Conclusion 389
In this study, an ensemble HA forecast model for the LaTex Shelf is developed using the state-of-the-art statistic programming 390 language R. The model is trained using numeric simulations from 1 January 2007 to 26 August 2020 generated by a coupled 391 hydrodynamic-biogeochemical model. Before splitting data into a training set and a test set, we applied regional average over 392 the LaTex Shelf and min-max normalization to the hindcast data. The ensemble model is then migrated to the GOFS 3.1 products based on HYCOM,s which provides eight-day forecast of 404 global hydrodynamics. The ensemble model is trained using the ROMS training set and then is used for the HA prediction 405 covering the period from January 1st, 2019 to August 26th, 2020. The prediction is robust when compared to the ROMS 406 simulations (2019-2020), with a low overall RMSE (4,571 km 2 ) and a high R 2 (0.8178). The model can also accurately predict 407 the magnitude and onset of summer HA peaks in 2019 and 2020, respectively. To our best knowledge, this ensemble model is 408 the first model providing efficient yet accurate daily HA forecast for the LaTex Shelf while considering both hydrodynamic 409 and biochemical effects. This model is also the first model successfully applying global hydrodynamic forecast in regional HA 410 predictions. 411