Factors Controlling Hypoxia Occurrence in Estuaries, Chester River, Chesapeake Bay

: The Chester River, a tributary of Chesapeake Bay, provides critical habitats for numerous living species and oyster aquaculture, but faces increasing anthropogenic stresses due to excessive nutrient loading and hypoxia occurrence. An application of the Integrated Compartment Water Quality Model (ICM), coupled with the Finite-Volume Community Ocean Model (FVCOM), was carried out to study the controlling mechanisms and interannual variability in hypoxia occurrence from 2002 to 2011. Our study shows that hypoxia occurs mostly in the main stem in July, followed by August and June. On an interannual scale, 2005 had the highest hypoxia occurrence with an accumulative hypoxia volume of about 10 km 3 -days, whereas 2008 had the lowest occurrence with an accumulative hypoxia volume of about 1 km 3 -days. Nutrient loading is the predominant factor in determining the intensity and interannual variability in hypoxia in the Chester River estuary, followed by stratiﬁcation and saltwater intrusion. Phosphorus has been found to be more e ﬃ cient in controlling hypoxia occurrence than nitrogen due to their di ﬀ erent limiting extent. On a local scale, the Chester River estuary is characterized by several meanders, and at certain curvatures helical circulation is formed due to centrifugal forces, leading to better reaeration and dissolved oxygen (DO) supply to the deeper layers. Our study provides valuable information for nutrient management and restoration e ﬀ orts in the Chester River.


Introduction
Hypoxia in oceanography is conventionally defined as a condition of low dissolved oxygen (<2 mg L −1 ), where living organisms are under stress [1,2]. This is commonly referred to as a Dead Zone [3]. Stratification, which reduces vertical water exchange and dissolved oxygen (DO) supply to deeper layers, is a physical condition for hypoxia development, but nutrient loadings resulted from anthropogenic activities lead to an ever-increasing trend in eutrophication and hypoxia occurrence in the coastal oceans and estuaries [4][5][6][7]. The North Sea, the Baltic Sea, and the Adriatic Sea are examples among other larger bodies of water that are subject to severe seasonal and episodic hypoxia events [2,8,9]. Elevated nitrogen concentrations are present in 28% of the U.S. stream length and 40% with elevated phosphorus concentration [10]. The Northern Gulf of Mexico, the Chesapeake Bay, the Long Island Sound, the Narragansett Bay, and the Puget Sound are examples of large bodies of water that regularly experience serious hypoxia events and ecosystem deterioration in the US. The Chesapeake Bay, located on the east coast along the Middle Atlantic Bight, is the largest estuary in the US. One major characteristic of the Chesapeake Bay is its large drainage watershed as compared to the Bay surface, 166.5 × 10 3 vs. 11.4 × 10 3 km 2 [11]. The population in the surroundings has more than doubled since 1950 to about 18 million, as have the nutrient loadings to the Bay [12,13]. Even though hypoxia has been observed since the early 1900s [14,15], hypoxic water volume has tripled from 1950 to 2000 [16], with oyster landing having declined by up to 6 folds [17,18]. the diffusivity in both horizontal and vertical). Zooplankton were not explicitly simulated and the grazing on phytoplankton was parameterized as a closure term by applying the grazing coefficient (α p ) to the phytoplankton biomass. This term included all losses, such as filter feeder clearance and aggregation to larger particulates, so that the quadratic function was used. In a bottom-up control system where the predator abundance covaries with the prey, the predation loss of the prey takes a mass-dependent quadratic function of the prey [39,40]. Temperature, light, and nutrient control on phytoplankton growth rate (µ) were parameterized based on the Jassby and Platt function [41]: where K T is an exponential coefficient for temperature control on phytoplankton growth, K I a constant in the Jassby-Platt light function, K N and K P are the half-saturation constants for nitrogen (N) and phosphorous (P), T O is the optimal temperature. The source and sink terms of dissolved organic carbon (DOC) without the transport terms are where a DOC is the DOC fraction of phytoplankton metabolism and predation losses, α LPOC and α LDOC are dissolution coefficients of labile and refractory particulate organic carbon (LPOC and RPOC), K OC and K DN are the half-saturation coefficients for DOC remineralization and denitrification, and α DOC and α DOC are DO-based and nitrate-based remineralization coefficients of DOC. The last term in Equation (3) determines the denitrification rate in the model that occurs at low DO and high nitrate concentration. The remineralization of particulate organic carbon (LPOC and RPOC) is through the dissolution (or hydrolysis) to DOC followed by the remineralization of the latter, so that the source terms are the metabolism and predation of phytoplankton, and sink terms are the dissolution and sinking to the bottom. Dissolved and particulate organic nitrogen and phosphorus are simulated in a similar way to carbon, with different rates of geochemical processes. As mentioned earlier, inorganic nitrogen cycles involved nitrification and denitrification. The ammonium equation is written as where a NC is the N:C ratio in phytoplankton, α NH , α DON , and α NT are the coefficients for metabolism leading to NH4, DON mineralization, and nitrification rate, p NH is the ammonium uptake preference by phytoplankton, K ONT and K NNT are the half-saturation of DO and NH4 for nitrification, and k NT is the exponential temperature coefficient for nitrification. The last term, nitrification, depends on DO and NH4 concentration and temperature (NH4 represents a variable in the model and not the chemical form of ammonium). Dissolved oxygen (DO) constitutes the core water quality variable in the ICM. DO is determined by photosynthesis production, consumption by phytoplankton respiration, DOC remineralization, nitrification, COD (Chemical Oxygen Demand) oxidation, and reaeration at the sea surface: where a OC is the oxygen to carbon ratio in phytoplankton and organic matter, p NH is NH4 sustained primary production (ammonium preference), a CD fraction of metabolism leading to DOC, a ON oxygen Water 2020, 12, 1961 4 of 17 to nitrogen ratio in nitrification, a DOC DOC remineralization rate, a COD COD oxidation rate, K ODC and K COD half-saturation of DO for DOC remineralization and COD oxidation, a air reaeration rate at the sea surface, DOs DO saturation and DZs the thickness of the surface layer. The reaeration rate was linked to wind speed (V w ) to the 1.5th power and surface water temperature (T) and salinity (S) due to their influence on the viscosity [39]: a air = a 0 V w 1.5 (0.54 + 0.0233T − 0.002S). DO saturation concentration was computed using the Garcia and Gordon formulation [42]. One mole of carbon fixation released 1.3 mole of oxygen based on nitrate assimilation, but only 1 mole of oxygen based on ammonium (constants in the first term of Equation (5)) [43]. Parameter definition and values are given in Table 1. ICM uses the Ditoro and Fitzpatrick model as the diagenesis module, which simulates the sinking flux of organic matter to the sediment, remineralization, burial, and fluxes of nutrients and COD from the sediment to the water column [44]. Two sediment layers (aerobic and anaerobic) were simulated and organic substances were divided into three categories based on the rate of remineralization [45]: G1, G2, and G3. Basically, G1 is labile (half-life about 20 days), G2 is refractory (half-life about 1 year), and G3 is almost inert (half-life about 30 years). Linear approaches were implemented in terms of decay from organic particulate substances to dissolve inorganic components and the decay rates were linked to temperature. Nutrient fluxes at the sediment-water interface were computed with an exchange coefficient applied to the difference in concentration between the overlying water and the sediment porewater. Fluxes of methane and sulfite were expressed together in oxygen equivalent as chemical oxygen demand (COD).
The simulation domain covers from the fall line to the Chester River mouth (Figure 1). Grid resolution varies from 1 km at the river mouth to 100 m near the coastline, with 13,824 cells in total, 8351 nodes, and 10 vertical sigma layers. The simulation time step was set at 5 s. The model was first spun up for 5 years using the forcing and currents of 2002. Initial condition was based on the

Forcing Data
ICM uses the FVCOM computed currents and turbulence diffusivity fields to drive the water quality simulation [39]. In addition to the physical forcing, ICM needs shortwave radiation for photosynthesis and wind speed to determine the DO reaeration at the sea surface. Wind speed data were from the NOAA Thomas Point Buoy (38°53'56" N 76°26'9" W) located in the northern part of the Chesapeake Bay close to the Chester River estuary, and shortwave radiation data were downloaded from the North America Regional Reanalysis (NARR) ftp site (ftp://ftp.cdc.noaa.gov/Datasets/ NARR). A factor of 0.47 was applied to shortwave radiation to convert to Photosynthetically Active Radiation (PAR), including reflection at the sea surface [46]. NARR provided data with 3 h intervals and a linear interpolation was conducted during the simulation.
Fall line data of nutrient loading was generated using the regulatory watershed model HSPF calibrated with data collected at the USGS River Input Monitoring (RIM) stations on the Chesapeake Bay watershed. The annual total nitrogen loading to the Chester River domain ranged from 0.8 to 2.5  10 6 kg with an average of 1.5  10 6 kg, and that of phosphorus ranged from 30 to 178  10 3 kg with an average of 117  10 3 kg. Open boundary conditions at the river month were based on the CH3D-ICM simulation. The CH3D-ICM simulations were calibrated over the entire Chesapeake Bay from 2001 to 2011 with extensive monitoring data from the Chesapeake Bay Program. State variables at the Chester River open boundary nodes were recorded hourly. As CH3D-ICM and FVCOM-ICM shared the same nodes at the Chester River mouth, the CH3D-ICM outputs were directly used for the FVCOM-ICM simulation without horizontal interpolation. However, CH3D-ICM uses a z-coordinate in the vertical, and a linear interpolation was performed to bring the boundary data from the CH3D-ICM z-layers to the FVCOM-ICM sigma layers. As the FVCOM-ICM uses a smaller time step (5 s) as compared to the CH3D-ICM output (hourly), a linear interpolation in time was also performed during the simulation.

Observation Data and Statistical Analysis
The Chesapeake Bay Program has maintained a monitoring program over the entire Bay since 1984. Two of the long-term monitoring stations are in the Chester River domain: ET4.2 in the mesohaline region (salinity 10-18 psu), and ET4.1 in the tidal fresh area (salinity <0.5 psu; Figure 1).

Forcing Data
ICM uses the FVCOM computed currents and turbulence diffusivity fields to drive the water quality simulation [39]. In addition to the physical forcing, ICM needs shortwave radiation for photosynthesis and wind speed to determine the DO reaeration at the sea surface. Wind speed data were from the NOAA Thomas Point Buoy (38 • 53 56 N 76 • 26 9 W) located in the northern part of the Chesapeake Bay close to the Chester River estuary, and shortwave radiation data were downloaded from the North America Regional Reanalysis (NARR) ftp site (ftp://ftp.cdc.noaa.gov/Datasets/ NARR). A factor of 0.47 was applied to shortwave radiation to convert to Photosynthetically Active Radiation (PAR), including reflection at the sea surface [46]. NARR provided data with 3 h intervals and a linear interpolation was conducted during the simulation.
Fall line data of nutrient loading was generated using the regulatory watershed model HSPF calibrated with data collected at the USGS River Input Monitoring (RIM) stations on the Chesapeake Bay watershed. The annual total nitrogen loading to the Chester River domain ranged from 0.8 to 2.5 × 10 6 kg with an average of 1.5 × 10 6 kg, and that of phosphorus ranged from 30 to 178 × 10 3 kg with an average of 117 × 10 3 kg. Open boundary conditions at the river month were based on the CH3D-ICM simulation. The CH3D-ICM simulations were calibrated over the entire Chesapeake Bay from 2001 to 2011 with extensive monitoring data from the Chesapeake Bay Program. State variables at the Chester River open boundary nodes were recorded hourly. As CH3D-ICM and FVCOM-ICM shared the same nodes at the Chester River mouth, the CH3D-ICM outputs were directly used for the FVCOM-ICM simulation without horizontal interpolation. However, CH3D-ICM uses a z-coordinate in the vertical, and a linear interpolation was performed to bring the boundary data from the CH3D-ICM z-layers to the FVCOM-ICM sigma layers. As the FVCOM-ICM uses a smaller time step (5 s) as compared to the CH3D-ICM output (hourly), a linear interpolation in time was also performed during the simulation.

Observation Data and Statistical Analysis
The Chesapeake Bay Program has maintained a monitoring program over the entire Bay since 1984. Two of the long-term monitoring stations are in the Chester River domain: ET4.2 in the mesohaline region (salinity 10-18 psu), and ET4.1 in the tidal fresh area (salinity <0.5 psu; Figure 1). Data of Water 2020, 12,1961 6 of 17 DO and chlorophyll from monthly (in winter) and biweekly (summer) cruises are available for the entire simulation period from 2002 to 2011 ( Figure 2). The Department of Natural Resource (DNR) of Maryland also undertakes monitoring activities in the Maryland portion of the Bay, including the Chester River and the Corsica River in the simulation domain. Some of the stations are occupied in long terms, such as XGG8251 located on the southern bank of the Chester River estuary (Figure 1), but most of the DNR monitoring stations are rotational, occupied during certain years. There were 23 stations sampled during the simulation period in the DNR monitoring programs and data availability and duration are depicted in Figure 2. Data of DO and chlorophyll from monthly (in winter) and biweekly (summer) cruises are available for the entire simulation period from 2002 to 2011 ( Figure 2). The Department of Natural Resource (DNR) of Maryland also undertakes monitoring activities in the Maryland portion of the Bay, including the Chester River and the Corsica River in the simulation domain. Some of the stations are occupied in long terms, such as XGG8251 located on the southern bank of the Chester River estuary ( Figure 1), but most of the DNR monitoring stations are rotational, occupied during certain years. There were 23 stations sampled during the simulation period in the DNR monitoring programs and data availability and duration are depicted in Figure 2. In addition to cruise-based monitoring programs, DNR conducted Continuous Monitoring programs (CMON) at nine of the monitoring stations (solid line in Figure 2). These are high-frequency (every 15 min) data based on electronic sensor measurements, and tele-communicated to land-based laboratories. The sensors were carefully calibrated on a monthly basis with bottle sample analysis. All these data were used for model calibration and validation. Time-series plots, Taylor Diagrams, and Target Diagrams were employed for simulation-observation comparison, and Principal Component Analysis (PCA) and Generalized Additive Model (GAM) were used for model results analysis. These methods were described in [39], where readers are referred to for more information. Briefly, Taylor diagrams depict the comparison between simulations and observations in terms of the correlation coefficient, the standard deviation of both the simulation and the observation, and the centered root mean squared error (CRMSE) on the same diagram [47,48]. The angle from zero indicates the correlation between simulation and observation, and the distance from the origin is the normalized standard deviation (std) of the simulated values (simulation std divided by observation std), with the result that the overall distance between the observation and the simulation on the Taylor diagram is the CRMSE. Note that the correlation coefficient (R) is used in the Taylor diagram whereas the coefficient of determination (R 2 ) is often referred to in statistical analyses. In the Target diagram, the bias (the difference between the simulation mean and the observation mean) is scaled on the y axis and the centered mean squared error (unbias) on the x axis, and as a result the total bias (the root mean squared error, RMSE) is the radial distance from the origin to the data point [48,49]. As in the Taylor diagram, the bias and the CRMSE are normalized to the data standard deviation. The combined Taylor and Target diagram can thus give a comprehensive assessment of the In addition to cruise-based monitoring programs, DNR conducted Continuous Monitoring programs (CMON) at nine of the monitoring stations (solid line in Figure 2). These are high-frequency (every 15 min) data based on electronic sensor measurements, and tele-communicated to land-based laboratories. The sensors were carefully calibrated on a monthly basis with bottle sample analysis. All these data were used for model calibration and validation. Time-series plots, Taylor Diagrams, and Target Diagrams were employed for simulation-observation comparison, and Principal Component Analysis (PCA) and Generalized Additive Model (GAM) were used for model results analysis. These methods were described in [39], where readers are referred to for more information. Briefly, Taylor diagrams depict the comparison between simulations and observations in terms of the correlation coefficient, the standard deviation of both the simulation and the observation, and the centered root mean squared error (CRMSE) on the same diagram [47,48]. The angle from zero indicates the correlation between simulation and observation, and the distance from the origin is the normalized standard deviation (std) of the simulated values (simulation std divided by observation std), with the result that the overall distance between the observation and the simulation on the Taylor diagram is the CRMSE. Note that the correlation coefficient (R) is used in the Taylor diagram whereas the coefficient of determination (R 2 ) is often referred to in statistical analyses. In the Target diagram, the bias (the difference between the simulation mean and the observation mean) is scaled on the y axis and the centered mean squared error (unbias) on the x axis, and as a result the total bias (the root mean squared error, RMSE) is the radial distance from the origin to the data point [48,49]. As in the Taylor diagram, the bias and the CRMSE are normalized to the data standard deviation. The combined Taylor and Target diagram can thus give a comprehensive assessment of the simulation by providing the correlation coefficient, standard deviation, the centered mean squared error, the bias, and the total RMSE (see "Results" section for graphic examples). The correlation coefficient essentially compares the timing of phytoplankton blooms and the seasonality of biogeochemical events, the bias compares the magnitudes and CRMSE compares the variability between simulation and observation. Harding et al. [50] found that GAM is suitable to extract long-term trend, interannual and seasonal variations from time series data, which was adopted in this paper. Principal component analysis (PCA) is a type of linear transformation through which principal dimensions are found in the variable space that explain the covariance among the variables [51].  simulation by providing the correlation coefficient, standard deviation, the centered mean squared error, the bias, and the total RMSE (see "Results" section for graphic examples). The correlation coefficient essentially compares the timing of phytoplankton blooms and the seasonality of biogeochemical events, the bias compares the magnitudes and CRMSE compares the variability between simulation and observation. Harding et al. [50] found that GAM is suitable to extract longterm trend, interannual and seasonal variations from time series data, which was adopted in this paper. Principal component analysis (PCA) is a type of linear transformation through which principal dimensions are found in the variable space that explain the covariance among the variables [51].

Comparison between Simulation and Observation
Data of dissolved oxygen and chlorophyll collected at 24 stations were used to assess the robustness of the simulation. Given the large quantity of data, only the most representative stations are plotted in time-series (   Given the shallow depth at Station ET4.1, only the surface layer is presented (Figure 3a,b). The time series of DO was dominated by seasonal cycle both in the observation and simulation, with high values in winter and low values in summer and transitional in spring and fall (Figure 3a). In addition to the seasonal cycles, the model revealed significant high-frequency variations in DO at a diurnal scale. In some years, the model tended to overestimate DO concentration at both the high and low ends. A particularly low value of DO (ca. 3 mg L −1 ) was observed in May 2008, which was not produced in the model simulation. The minimum values of DO in 2003 and 2004 were also lower in the observation than in the simulation.
In the chlorophyll time series data at Station ET4.1, (Figure 3b), two peaks of concentration were observed and simulated in most of the years: one in April and another in October. This indicates two phytoplankton blooms in the region: spring bloom in April and fall bloom in October. The spring phytoplankton bloom had a higher amplitude with chlorophyll concentration up to 150 to 200 µg L −1 , whereas the fall bloom had a relatively lower amplitude with chlorophyll concentration close to 100 µg L −1 . The lowest values were reached during the winter season between two successive years. High-frequency variations superposed on the seasonal cycle were also reproduced in the chlorophyll concentration as well as in the DO concentration. In most of the cases, the model simulation compared reasonably well with the observation. In 2009, however, the observation did not show pronounced blooms, whereas the model predicted two blooms as in other years. Additionally, the spring phytoplankton bloom was missing in the observation in 2010 and 2011 but was generated in the simulation.
DO at Station ET4.2 was also dominated by large seasonal variations, with high values in winter and spring, and low values in summer and fall (Figure 3c). The highest values were mostly found in spring during the phytoplankton bloom. At this station, the model simulation of DO compared well with the observation in both the timing and amplitude of the seasonal cycles. The model predicted two phytoplankton blooms at this station as well (Figure 3d), but with an amplitude much smaller than that at Station ET4.1. The spring bloom reached about 70 µg L −1 in chlorophyll concentration and about 50 µg L −1 for the fall bloom, which were less than half of the amplitudes at Station ET4.1. DO simulation in the bottom layer also compared well with the observation (Figure 3e). The amplitude of seasonal variation in the bottom DO was even higher than that in the surface layer, ranging from 0 to 15 mg L −1 .
Taylor Diagrams and Target Diagrams were constructed on DO simulation at all the observation stations ( Figure 4). For surface DO, most of the stations had a correlation coefficient between simulation and observation >0.6 and the centered mean squared error <1 standard deviation of the observation (Figure 4a). Stations 0348 and ET4.1 had a correlation coefficient lower than 0.6 and Station 0077 had a centered root mean squared error higher than 1. These stations are all located in the tidal fresh zone where DO concentration is relatively high (>2 mg L −1 , the critical value for hypoxia). As show in Figure 3a, the model tended to overestimate DO concentration at both the high and low ends in the DO concentration range, which may explain the lower statistical numbers. On the Target diagram (Figure 4b), most of the stations had a total root mean squared error (RMSE) lower than or close to 1 standard deviation of the observation. Again, Stations 0077 and 0348 were outliers with a total RMSE higher or close to 2. At these two stations, the model tended to overestimate the mean of DO concentration. Similar results were obtained for the bottom layer (Figure 4c,d). Most of the stations had a correlation coefficient >0.6 and a centered RMSE less than 1. Stations 0348 and ET4.1 were among the stations with a lower correlation coefficient (<0.6) and Station 0077 had a relatively larger centered RMSE. On the Target Diagram, most of the stations had a total RMSE <1, with some other stations being between 1 and 2. Station 6496 was an outlier with a total RMSE higher than 2. Both the mean and the variability were underestimated by the model at this station. This station is in a secondary tributary on the northern bank of the Chester River and it is unlikely that it can significantly contribute to the total hypoxia volume in the entire Chester River estuary. Water 2019, 11, x FOR PEER REVIEW 9 of 17

Seasonality and Interannual Variations in Hypoxia Volume
Hypoxia is defined as a condition of low dissolved oxygen (<2 mg L −1 ) in a marine environment [1,2,9,52]. Instantaneous and monthly average of hypoxia was given in total volume in the Chester River estuary, while interannual variations were expressed as hypoxia volume-days, i.e., the sum of daily hypoxia volume over the entire year. The timing of hypoxia occurrence can differ from year to year so that changes in daily and monthly hypoxia volume can be caused by shifting in time of the hypoxia occurrence. Hypoxia volume-days integrate all the hypoxia occurrences during the entire year so that it is less subject to changes in the seasonality of hypoxia events. The year 2005 had the highest hypoxia volume-days (black bar in Figure 5

Seasonality and Interannual Variations in Hypoxia Volume
Hypoxia is defined as a condition of low dissolved oxygen (<2 mg L −1 ) in a marine environment [1,2,9,52]. Instantaneous and monthly average of hypoxia was given in total volume in the Chester River estuary, while interannual variations were expressed as hypoxia volume-days, i.e., the sum of daily hypoxia volume over the entire year. The timing of hypoxia occurrence can differ from year to year so that changes in daily and monthly hypoxia volume can be caused by shifting in time of the hypoxia occurrence. Hypoxia volume-days integrate all the hypoxia occurrences during the entire year so that it is less subject to changes in the seasonality of hypoxia events. The year 2005 had the highest hypoxia volume-days (black bar in Figure 5 Daily hypoxia volume is depicted in Figure 6  Daily hypoxia volume is depicted in Figure 6 for   The spatial distribution of bottom DO in June 21 is depicted as an example of low hypoxia day and in July 21 as an example of high hypoxia day in 2005 (Figure 8). In June 21, only a small area in the deep channel adjacent to the Corsica River had bottom DO <2 mg L −1 . In July 21, low DO <2 mg L −1 became widespread, from the lower estuary to the upper reach of the estuary. The core of hypoxia distribution is along the main deep channel. What is interesting is that hypoxia distribution is  Daily hypoxia volume is depicted in Figure 6    The spatial distribution of bottom DO in June 21 is depicted as an example of low hypoxia day and in July 21 as an example of high hypoxia day in 2005 (Figure 8). In June 21, only a small area in the deep channel adjacent to the Corsica River had bottom DO <2 mg L −1 . In July 21, low DO <2 mg L −1 became widespread, from the lower estuary to the upper reach of the estuary. The core of hypoxia distribution is along the main deep channel. What is interesting is that hypoxia distribution is  Daily hypoxia volume is depicted in Figure 6    The spatial distribution of bottom DO in June 21 is depicted as an example of low hypoxia day and in July 21 as an example of high hypoxia day in 2005 (Figure 8). In June 21, only a small area in the deep channel adjacent to the Corsica River had bottom DO <2 mg L −1 . In July 21, low DO <2 mg L −1 became widespread, from the lower estuary to the upper reach of the estuary. The core of hypoxia distribution is along the main deep channel. What is interesting is that hypoxia distribution is The spatial distribution of bottom DO in June 21 is depicted as an example of low hypoxia day and in July 21 as an example of high hypoxia day in 2005 (Figure 8). In June 21, only a small area in the deep channel adjacent to the Corsica River had bottom DO <2 mg L −1 . In July 21, low DO <2 mg L −1 became widespread, from the lower estuary to the upper reach of the estuary. The core of hypoxia distribution is along the main deep channel. What is interesting is that hypoxia distribution is discontinued at several curvatures of the main channel. The bathymetry did not show similar discontinuity in a similar manner at the same location. It is most likely linked to other mechanisms that have the potential to alter DO concentration in the bottom layer. discontinued at several curvatures of the main channel. The bathymetry did not show similar discontinuity in a similar manner at the same location. It is most likely linked to other mechanisms that have the potential to alter DO concentration in the bottom layer.

GAM Decomposition of Time-Series Hypoxia Volume and PCA Causality Analysis
General additive models were fitted to the daily hypoxia volume time series (Figure 9). No general trend was found. The hypoxia time series is dominated by seasonal variations, with seasonal model accounting for 53% of the total variance. Most of the hypoxia occurred within a relatively short period of time in June and July. Adding interannual variations resulted in a fitted function that explains 59% of the total variance. The prediction of the interannual function (red curve, Figure 9

GAM Decomposition of Time-Series Hypoxia Volume and PCA Causality Analysis
General additive models were fitted to the daily hypoxia volume time series (Figure 9). No general trend was found. The hypoxia time series is dominated by seasonal variations, with seasonal model accounting for 53% of the total variance. Most of the hypoxia occurred within a relatively short period of time in June and July. Adding interannual variations resulted in a fitted function that explains 59% of the total variance. The prediction of the interannual function (red curve, Figure 9 discontinued at several curvatures of the main channel. The bathymetry did not show similar discontinuity in a similar manner at the same location. It is most likely linked to other mechanisms that have the potential to alter DO concentration in the bottom layer.

GAM Decomposition of Time-Series Hypoxia Volume and PCA Causality Analysis
General additive models were fitted to the daily hypoxia volume time series (Figure 9). No general trend was found. The hypoxia time series is dominated by seasonal variations, with seasonal model accounting for 53% of the total variance. Most of the hypoxia occurred within a relatively short period of time in June and July. Adding interannual variations resulted in a fitted function that explains 59% of the total variance. The prediction of the interannual function (red curve, Figure 9   . General Additive Model fitting to daily hypoxia volume. Dots: daily data; black line: trend (y = gam(s(time, bs = "cr")) where time is expressed in decimal years) "cr" indicates "cubic regression spline smooth"; blue line: seasonal variability (y = gam(s(DOY, bs = "cr")) where DOY stands as day of the year); red line: Interannual variability (y = gam(s(time, bs = "cr") + s(DOY, bs = "cr") + ti(time, DOY, bs = "cr")) where ti is the tensor product interaction function of time and DOY).
Principal component analysis (PCA) was conducted on the daily data of total hypoxia volume (HPOV), nutrient loads including dissolved inorganic nitrogen (DIN), total nitrogen (TN), phosphate (PO4) and total phosphorus (TP), bottom saltwater intrusion distance (Bottom), difference between bottom and surface saltwater intrusion (Diff ), and stratification using the Brunt-Väisälä frequency (N 2 ). The N 2 was the daily average on the along-channel transect. The first principal component explained 56% of the total variance (x axis in Figure 10) and the second accounted for another 23% of the total variance (y axis in Figure 10). As such, the first two principal components explained 79% of the total variance. All the variables had high loadings on the first principal component and were thus relevant to the hypoxia development. Except the bottom saltwater intrusion distance, all other variables were located on the same side as HPOV on the first principal component, indicating their positive relationships with hypoxia occurrence. On the second principal component, PO4 and TP were located on the same side as HPOV whereas DIN and TN were on the opposite side, indicating that phosphorus loads have more impact on the hypoxia development than nitrogen loads in the Chester River estuary. Both the stratification N 2 and Diff were loaded on the same side as the HPOV on the first principal component, indicating their positive impact on hypoxia development. However, saltwater intrusion had a negative relationship with HPOV. Principal component analysis (PCA) was conducted on the daily data of total hypoxia volume (HPOV), nutrient loads including dissolved inorganic nitrogen (DIN), total nitrogen (TN), phosphate (PO4) and total phosphorus (TP), bottom saltwater intrusion distance (Bottom), difference between bottom and surface saltwater intrusion (Diff), and stratification using the Brunt-Väisälä frequency (N 2 ). The N 2 was the daily average on the along-channel transect. The first principal component explained 56% of the total variance (x axis in Figure 10) and the second accounted for another 23% of the total variance (y axis in Figure 10). As such, the first two principal components explained 79% of the total variance. All the variables had high loadings on the first principal component and were thus relevant to the hypoxia development. Except the bottom saltwater intrusion distance, all other variables were located on the same side as HPOV on the first principal component, indicating their positive relationships with hypoxia occurrence. On the second principal component, PO4 and TP were located on the same side as HPOV whereas DIN and TN were on the opposite side, indicating that phosphorus loads have more impact on the hypoxia development than nitrogen loads in the Chester River estuary. Both the stratification N 2 and Diff were loaded on the same side as the HPOV on the first principal component, indicating their positive impact on hypoxia development. However, saltwater intrusion had a negative relationship with HPOV.

Time-Series Pattern
The seasonal cycle of DO is due to changes in the solubility of DO in response to changes in water temperature and changes in biological and biogeochemical processes. Water temperature ranges from 1 °C in winter to about 30 °C in summer [38]. Based on the DO saturation formulation of Garcia and Gordon, [42] used in the model and assuming a salinity value of 2 psu at Station ET4.1, the DO saturation concentration varies from 14.0 mg L −1 at 1 °C to 7.5 mg L −1 at 30 °C., i.e., the saturation concentration at 1 °C is almost double that at 30 °C. The simulated variation in DO concentration ranges from 5 to 15 mg L −1 , which is larger than the changes in solubility. This extra seasonal variation in DO concentration is linked to seasonal shifts in biological and biogeochemical dynamics. Spring is dominated by phytoplankton blooms, which produce oxygen and lead to supersaturation in DO concentration. Summer and early fall are dominated by remineralization of organic substances, which consume DO and lead to under-saturation in DO concentration. The diurnal variations on top of the seasonal cycles are due to the diurnal cycle between phytoplankton photosynthesis during the daytime and respiration and remineralization during the night. Other

Time-Series Pattern
The seasonal cycle of DO is due to changes in the solubility of DO in response to changes in water temperature and changes in biological and biogeochemical processes. Water temperature ranges from 1 • C in winter to about 30 • C in summer [38]. Based on the DO saturation formulation of Garcia and Gordon, [42] used in the model and assuming a salinity value of 2 psu at Station ET4.1, the DO saturation concentration varies from 14.0 mg L −1 at 1 • C to 7.5 mg L −1 at 30 • C., i.e., the saturation concentration at 1 • C is almost double that at 30 • C. The simulated variation in DO concentration ranges from 5 to 15 mg L −1 , which is larger than the changes in solubility. This extra seasonal variation in DO concentration is linked to seasonal shifts in biological and biogeochemical dynamics. Spring is dominated by phytoplankton blooms, which produce oxygen and lead to supersaturation in DO concentration. Summer and early fall are dominated by remineralization of organic substances, which consume DO and lead to under-saturation in DO concentration. The diurnal variations on top of the seasonal cycles are due to the diurnal cycle between phytoplankton photosynthesis during the daytime and respiration and remineralization during the night. Other mechanisms that can lead to high-frequency variations are heterogeneity in the water property, such as DO concentration and phytoplankton abundance, and horizontal advection. Tide and river flow continuously move water mass through a fixed station such as ET4.1, and differences in water property will lead to variations in observation and simulation at a fixed station.
The amplitude of seasonal variability in bottom DO at ET4.2 is even higher than that in the surface layer, ranging from 0 to 15 mg L −1 . In addition to changes in solubility, stratification is another factor affecting DO concentration in the bottom layer. During the winter season when the water column is quasi-homogeneous, bottom water is mostly mixed with surface water with a high DO concentration. During the summer season when the water column is stratified, bottom water is stagnated with little convection or mixing with the surface layers. Moreover, high temperature and abundant organic substances produced by the spring phytoplankton bloom contribute to the increased remineralization and thus DO depletion in summer, amplifying the seasonal variations in bottom DO. The high similarity between bottom DO observation and simulation provides sound basis for robust model prediction of hypoxia volume.

Major Drivers of Hypoxia Occurrence
PCA causality analysis revealed that nutrient loads are the major driver of hypoxia development in the Chester River estuary. This agrees with previous studies conducted in the entire Chesapeake Bay where nutrient loading is found to be the dominating factor in determining interannual variability in hypoxia occurrence [16,22]. However, phosphorus and nitrogen loads do not have equal impact. Phosphorus loads have a higher effect than nitrogen loads on hypoxia development in the Chester River estuary. This is due to the difference in limitation between the two types of nutrient elements. Usually phosphorus is more limiting in terrestrial environments, whereas nitrogen is the essential limiting factor in the ocean [20]. Estuaries are transitional zones where the relative limiting effect can vary in space and time [53,54]. In the Chester River estuary, phosphorus is more of a limiting factor than nitrogen. Consequently, phosphorus loads have more impact on hypoxia development than nitrogen. This has an important implication in nutrient management and related decision making. Phosphorus load reduction will be more efficient than nitrogen in water quality restoration in the Chester River estuary. On top of nutrient loads, stratification and saltwater intrusion turned out to be significant factors in determining hypoxia occurrence in the Chester River estuary. Both the Brunt-Väisälä frequency (N 2 ) and the difference between bottom and surface saltwater intrusion distances reflect stratification of the water column. The negative relationship between stratification and bottom water saltwater intrusion means that stronger saltwater intrusion leads to less stratification and less hypoxia. When saltwater intrusion extended to the oligohaline zone (salinity 0.5-10 psu) in the middle portion of the Chester River estuary where water depth is relative shallow, bottom friction can cause strong vertical mixing and reduce the overall average of stratification along the main stem and less hypoxia development. Figure 8b shows the spatial distribution of hypoxic water with DO < 2 mg L −1 , and there is a major discontinuity at the curvature of the main stem. Similarly, the stratification metric N 2 is also relatively low at the same location and time (Figure 11). At meandering points, the water momentum and centrifugal force lead to the formation of helical circulation that can strengthen lateral circulation and vertical convection [38,[55][56][57][58]. Indeed, there is stronger lateral circulation on the transect at the curvature where hypoxia distribution is interrupted (Figure 12b). In the surface layer, water moves toward the right-hand shore (the outside bank of the meandering), which is in the opposite direction of the Coriolis force, but in the direction of the centrifugal force. Consequently, the outside bank of the curvature becomes a convergent zone whereas the inside bank a divergent zone, which subsequently generate cross channel circulation with elevated vertical vectors in the circulation field. On the other transect downstream of the meandering where the main channel is relatively straight, water mass mostly moves to the left side (looking up stream), in agreement with the Coriolis force. The major meandering of the main channel at the lower Chester River estuary causes helical circulation, enhances vertical convection and mixing, weakens stratification, and improves DO concentration and hypoxia occurrence.

Conclusions
The model simulation shows that hypoxia is the most severe in July in the Chester River, followed by August and June, whereas there is limited extension in hypoxia occurrence in May and September. The timing of the hypoxia occurrence changes from year to year, but mostly within the three-month window from June to August. In 2010, the total hypoxia volume in June was as high as in July and in 2009, the total hypoxia volume was even slightly higher in August than in July. Most of the hypoxia occurred in the main stem along the deep channels, with limited occurrence in the Corsica River and other sub-tributaries. The Chester River estuary is characterized by several meanders and it has been found that helical circulation is formed at certain curvatures due to the centrifugal force exerted on running water. The helical lateral circulation creates vertical convective vectors, reduces stratification, enhances exchange between the surface and bottom layers, and as such reduces the occurrence of hypoxia at certain locations. Hypoxia was the most severe in 2005 and the lowest in 2004. Statistical analysis shows that nutrient loading is the dominating factor in determining the extension and severity of hypoxia in the Chester River estuary, followed by stratification and saltwater intrusion. It has been found that phosphorus loads have more impact on hypoxia occurrence than nitrogen in the Chester River estuary due to its higher limiting effect on phytoplankton development. Consequently, phosphorus load reduction will be more efficient than nitrogen load reduction in restoration efforts.

Conclusions
The model simulation shows that hypoxia is the most severe in July in the Chester River, followed by August and June, whereas there is limited extension in hypoxia occurrence in May and September. The timing of the hypoxia occurrence changes from year to year, but mostly within the three-month window from June to August. In 2010, the total hypoxia volume in June was as high as in July and in 2009, the total hypoxia volume was even slightly higher in August than in July. Most of the hypoxia occurred in the main stem along the deep channels, with limited occurrence in the Corsica River and other sub-tributaries. The Chester River estuary is characterized by several meanders and it has been found that helical circulation is formed at certain curvatures due to the centrifugal force exerted on running water. The helical lateral circulation creates vertical convective vectors, reduces stratification, enhances exchange between the surface and bottom layers, and as such reduces the occurrence of hypoxia at certain locations. Hypoxia was the most severe in 2005 and the lowest in 2004. Statistical analysis shows that nutrient loading is the dominating factor in determining the extension and severity of hypoxia in the Chester River estuary, followed by stratification and saltwater intrusion. It has been found that phosphorus loads have more impact on hypoxia occurrence than nitrogen in the Chester River estuary due to its higher limiting effect on phytoplankton development. Consequently, phosphorus load reduction will be more efficient than nitrogen load reduction in restoration efforts.

Conclusions
The model simulation shows that hypoxia is the most severe in July in the Chester River, followed by August and June, whereas there is limited extension in hypoxia occurrence in May and September. The timing of the hypoxia occurrence changes from year to year, but mostly within the three-month window from June to August. In 2010, the total hypoxia volume in June was as high as in July and in 2009, the total hypoxia volume was even slightly higher in August than in July. Most of the hypoxia occurred in the main stem along the deep channels, with limited occurrence in the Corsica River and other sub-tributaries. The Chester River estuary is characterized by several meanders and it has been found that helical circulation is formed at certain curvatures due to the centrifugal force exerted on running water. The helical lateral circulation creates vertical convective vectors, reduces stratification, enhances exchange between the surface and bottom layers, and as such reduces the occurrence of hypoxia at certain locations. Hypoxia was the most severe in 2005 and the lowest in 2004. Statistical analysis shows that nutrient loading is the dominating factor in determining the extension and severity of hypoxia in the Chester River estuary, followed by stratification and saltwater intrusion. It has been found that phosphorus loads have more impact on hypoxia occurrence than nitrogen in the Chester River estuary due to its higher limiting effect on phytoplankton development. Consequently, phosphorus load reduction will be more efficient than nitrogen load reduction in restoration efforts.
Funding: This research was funded by the US Environmental Protection Agency (grant number CB-75230480) and the APC was funded by the Integration and Application Network, University of Maryland Center for Environmental Science.