Modeling Summer Hypoxia Spatial Distribution and Fish Habitat Volume in Artiﬁcial Estuarine Waterway

: This study analyzes the dissolved oxygen (DO) depletion or hypoxia formation affecting the ecological vulnerability of Gyeongin-Ara Waterway (GAW), an artiﬁcial estuarine waterway. The physical, chemical, and biochemical factors affecting the summer hypoxia dynamics and distribution are simulated and the habitat volumes of major ﬁsh species are calculated. CE-QUAL-W2, a two-dimensional hydrodynamic and water quality model, is applied for the simulation. Comparison with observation reveals that the salinity stratiﬁcation, vertical DO gradient, and summer hypoxia characteristics are realistically reproduced by the model. Comprehensive analysis of the spatial distributions of the residence time, salinity, and DO concentration reveal that the residence time is longest at the bottom of a freshwater inﬂow zone. Accordingly, residence time is identiﬁed as the physical factor having the greatest inﬂuence on hypoxia. It is also clear that a hypoxic water mass diffuses towards the entire waterway during neap tides and summer, when the seawater inﬂow decreases. Based on the modeling results, the DO depletion drivers are identiﬁed and the hypoxic zone formation and distribution are sufﬁciently explained. Finally, ﬁsh habitat volumes are calculated. In particular, the survival habitat volume of Mugil cephalus is found to decrease by 32–34% as a result of hypoxia from July to August. The model employed in this study could be utilized to establish an operational plan for the waterway, which would increase ﬁsh habitat volumes.


Introduction
The number of artificial waterways constructed in estuarine areas is increasing significantly worldwide. However, many of those artificial waterways are eutrophic, because of high-concentration nutrients from nearby cities, and are affected by water pollution due to long retention times and poor water circulation [1]. Hypoxia, i.e., dissolved oxygen (DO) depletion, is a representative water quality problem that has been drastically increasing in shallow coastal and estuarine areas, as well as in estuarine waterways worldwide [2][3][4][5]. Hypoxia, which is defined as a DO level of less than 62.5 µmol/L (2.0 mg/L), develops wherever the consumption of oxygen by organisms or chemical processes exceeds the oxygen supply from adjacent water, the atmosphere, and photosynthesizing organisms [2][3][4][5].
Oxygen depletion may have various causes. Specific hydrodynamic conditions (e.g., the limited water exchange and long water retention times observed in enclosed and semi-enclosed water bodies, insufficient oxygen supply due to density stratification, and presence of wind-shielded freshwater lenses in coastal seas near river mouths) may render aquatic systems naturally prone to oxygen depletion. In addition, oxygen depletion is often caused by the oxidation of organic matter and reduced substances as well as by the advective intrusion of oxygen-poor waters from deeper layers or adjacent water bodies. As noted above, the hypoxia threshold is defined as a DO level below 2 mg/L in many studies [2][3][4][5]; however, the actual threshold for hypoxia and its associated effects can vary over time and space, depending on water temperature and salinity [6,7]. Hypoxia also has diverse effects on aquatic ecosystems. Generally, a DO level of less than 5 mg/L is potentially hazardous to a marine ecosystem. At a level exceeding 3 mg/L but below 5 mg/L, bottom fish may begin to leave the area, and the growth of sensitive species will be reduced. For a DO level of less than 2 mg/L, some juvenile fish and crustaceans that cannot leave the area may die. Finally, at a DO level of less 1 mg/L, fish will totally avoid the area or begin to die in large numbers [6][7][8][9].
As DO depletion of an estuary area is the product of an intricate combination of physical, chemical, and biochemical processes, various methods have been adopted to model this phenomenon. For example, Nielsen et al. [10] used a simple model in a study of the shallow artificial estuary known as Ringkøbing Fjord in Denmark, revealing that DO depletion in deep parts of the estuary is caused by the long retention of intermittently inflowing seawater and the rapid consumption of oxygen under eutrophic conditions. Furthermore, when a three-dimensional model was applied to St. Lucie Estuary, a shallow subtropical estuary in the USA [11], excessive nutrients and algae from upstream freshwater were found to be responsible for algal bloom and DO depletion in the estuary. Bottom hypoxia due to stratification and deteriorated circulation was also noted. In addition, Chen et al. [12] assessed the physical factors (river discharge, wind speed, and direction) of hypoxia using a three dimensional circulation model in the Yangtze estuary in China. They showed that the area exhibiting hypoxia was greatly affected by wind speed and direction, but was insensitive to river discharge. Lajaunie-Salla et al. [13] also applied a modeling method to assess the impact of urban wastewater treatment plants (WWTPs) and sewage overflows (SOs) on hypoxia in the Gironde Estuary, France. The modeling results revealed that particulate organic carbon (POC) from urban waste was resuspended during the spring tide, thereby expanding the hypoxia area. Those researchers also found that hypoxic events could be mitigated by blocking the POC from the WWTPs and SOs. Overall, the findings to date can be summarized as follows: DO depletion of an estuarine area is mainly caused by excessive nutrients and organic matter, and the formation of a hypoxia zone is closely related to deteriorated circulation due to stratification and upstream freshwater inflow [10][11][12][13][14][15][16][17][18]. In addition, the spatial distribution of hypoxia in a small estuarine area is affected by river inflow [10,11,13], while that in a large estuarine area is influenced by wind [12,[14][15][16][17][18]. However, although diverse modeling studies have considered hypoxia in estuaries [10][11][12][13][14][15][16][17][18], and the impact of hypoxia on fish has been investigated in various ways [19][20][21][22][23][24], few studies have attempted to assess the spatial distribution of hypoxia in connection with fish habitats.
The present study has the following aims: (1) to analyze the impact of hydraulic factors and water quality on the spatiotemporal formation of summer hypoxia in an artificial waterway using CE-QUAL-W2, a two dimensional hydrodynamic and water quality model; and (2) to calculate the habitat volumes of major fish species in relation to temperature and oxygen to analyze seasonal vulnerability. The study area is Gyeongin-Ara Waterway (GAW) in South Korea, in which stagnant zones are hydraulically formed, salinity stratification occurs, and organic matter and nutrients are abundant. As a result, GAW has many water quality problems. The results of this study are expected to contribute to establishment of reasonable operational measures for an artificial waterway, which could mitigate existing water pollution problems.

Study Area
GAW is South Korea's first navigational waterway artificially constructed to connect the West Sea and Han River. It is located in the northern part of South Korea between the cities of Incheon and Seoul (the South Korean capital). GAW was originally designed to prevent chronic flooding around the Gulpo River region, but was redesigned to serve as both a transport route and a flood prevention facility. The main components of the waterway include an 18-km-long navigation channel (width: 80 m, operating depth: 6.3 m) connecting terminals (T/Ms) at Incheon (area: 2.84 km 2 ) and Gimpo (area: 1.91 km 2 ), and lock gates at each end ( Figure 1a). Additional structures include a diversion culvert and weir where the Gulpo River joins GAW. GAW transports seawater from the West Sea and freshwater from the Han River (through lock gates at the ends of the two terminals), and discharges water through the flood gate of the Incheon T/M. The Gulpo River is normally diverted to the Han River, and flows into the waterway only when the Gyulhyeon Weir is overflowing because of rainfall. The quantity of the saltwater from the West Sea is approximately twice that of the freshwater from the Han River. Therefore, a brackish water zone is artificially formed in the waterway. The tidal cycle range in the West Sea area, which is adjacent to the Incheon T/M, is between sea level +5.30 and −4.45 EL. m, and the waterway operation level is EL. 2.7 m. Accordingly, there are periods during the tidal cycle when water flows both into and out of the waterway. As a result of these large water level fluctuations and changes in flow direction, the waterway exhibits spatially complex hydraulic characteristics.

Study Area
GAW is South Korea's first navigational waterway artificially constructed to connect the West Sea and Han River. It is located in the northern part of South Korea between the cities of Incheon and Seoul (the South Korean capital). GAW was originally designed to prevent chronic flooding around the Gulpo River region, but was redesigned to serve as both a transport route and a flood prevention facility. The main components of the waterway include an 18-km-long navigation channel (width: 80 m, operating depth: 6.3 m) connecting terminals (T/Ms) at Incheon (area: 2.84 km 2 ) and Gimpo (area: 1.91 km 2 ), and lock gates at each end ( Figure 1a). Additional structures include a diversion culvert and weir where the Gulpo River joins GAW. GAW transports seawater from the West Sea and freshwater from the Han River (through lock gates at the ends of the two terminals), and discharges water through the flood gate of the Incheon T/M. The Gulpo River is normally diverted to the Han River, and flows into the waterway only when the Gyulhyeon Weir is overflowing because of rainfall. The quantity of the saltwater from the West Sea is approximately twice that of the freshwater from the Han River. Therefore, a brackish water zone is artificially formed in the waterway. The tidal cycle range in the West Sea area, which is adjacent to the Incheon T/M, is between sea level +5.30 and −4.45 EL. m, and the waterway operation level is EL. 2.7 m. Accordingly, there are periods during the tidal cycle when water flows both into and out of the waterway. As a result of these large water level fluctuations and changes in flow direction, the waterway exhibits spatially complex hydraulic characteristics. Although the basin area of GAW (including the Gulpo River area) does not exceed 134 km 2 , the Han River flowing in from the Gimpo T/M has a very large basin area of 25,953 km 2 . Further, as the Han River flows through Seoul, the country's capital with a population exceeding 10 million, its pollution levels are usually high. The Gulpo River flows into the waterway in the case of flooding only. However, the majority of the inflow is from urban WWTPs. Therefore, the Gulpo River has the poorest water quality among the three inflow sources (the Han River, West Sea, and Gulpo River). Consequently, GAW is eutrophic year-round, which causes water quality problems such as algal Although the basin area of GAW (including the Gulpo River area) does not exceed 134 km 2 , the Han River flowing in from the Gimpo T/M has a very large basin area of 25,953 km 2 . Further, as the Han River flows through Seoul, the country's capital with a population exceeding 10 million, its pollution levels are usually high. The Gulpo River flows into the waterway in the case of flooding only. However, the majority of the inflow is from urban WWTPs. Therefore, the Gulpo River has the poorest water quality among the three inflow sources (the Han River, West Sea, and Gulpo River).
Consequently, GAW is eutrophic year-round, which causes water quality problems such as algal bloom and DO depletion. Furthermore, as fish enter the waterway when the rock gates are opened, the waterway has a very peculiar environment in which freshwater, brackish, and saltwater fish species are observed simultaneously.

Water Quality
The GAW water quality was surveyed on a weekly basis from January to December 2014. The inflow water outside the waterway was monitored once a week at three sites, namely the Han River (outside the Gimpo T/M lock gate), Gulpo River (the upper region of Gyulhyeon Weir), and the West Sea (outside the Incheon T/M lock gate). Similarly to the case of the outside waterway, the water quality of the navigational channel was monitored weekly at four main sites, namely Gimpo T/M (S1), Gyeyang Bridge (S2), Sicheon Bridge (S3), and Incheon T/M (S4) (Figure 1a). Inflow water samples were collected from the surface layer using a Van Dorn water sampler. In the channel, samples were collected at 1, 3, and 5 m depth, corresponding to the surface, middle, and bottom layers, respectively. A water quality measurement instrument (YSI EXO1, Yellow Spring Instrument Inc.) was used to measure the water temperature, DO, salinity, pH, conductivity, and turbidity during the sampling. The chlorophyll-a (Chl-a), phosphate (PO 4 -P), nitrate (NO 3 -N), ammomium (NH 4 -N), total phosphorus (TP), and total nitrogen (TN) levels were analyzed using the standard method [25].

Fish
The fish abundance data for GAW used in this study were based on an environmental impact assessment survey by K-water [26]. The investigation was conducted in each season of 2014 at the four water quality monitoring sites of the waterway and in the Han and Gulpo Rivers. Details of sampling methods are provided in Appendix A.

Model Selection and Description
The CE-QUAL-W2 model is a two-dimensional laterally averaged hydrodynamic and water quality model co-developed by the U.S. Army Corps of Engineers Waterways Experiment Station and Portland State University [27]. Because the model assumes lateral homogeneity, it is best suited for relatively long and narrow waterbodies exhibiting longitudinal and vertical water quality gradients. This model has been applied to various stratified lakes and reservoirs worldwide [28][29][30][31][32]. In particular, Park et al. [33] and Brown et al. [34] have shown successful application of this model to estuaries. The CE-QUAL-W2 version 4.0 model was selected in this study because it is appropriate for GAW, for which the lateral variations in water quality are unimportant but the vertical variations in velocity, salinity, and DO are significant. CE-QUAL-W2 is based on finite difference solution of the laterally averaged equations of fluid motion, including the continuity equation, x-momentum equation, z-momentum equation, free surface equation, equation of state, and advective-diffusion equation. Details of these equations are provided in Appendix B.

Configuration and Application of Model
The geographical model of GAW constructed by CE-QUAL-W2 is illustrated in Figure 1b. The model grid consisted of 173 segments and 52 layers in the water depth direction. The average segment length was 103 m and the layers were finely distinguished at a vertical spacing of 0.25 m. For the inflow/outflow data, the GAW water management system (ARAWMS) constructed and input datasets representing inflow from the West Sea, Han River, and Gulpo River and datasets representing outflow through the West Sea flood gate at hourly intervals. The water quality of the inflow from the three sources was monitored every week, and the survey results were input to establish the water quality boundary conditions. The 12 constituents taken as the boundary conditions of CE-QUAL-W2 were the chemical oxygen demand (COD Mn ), salinity, suspended solid (SS), PO 4 -P, NH 4 -N, NO 3 -N, labile dissolved organic matter (LDOM), refractory dissolved organic matter (RDOM), labile particulate organic matter (LPOM), refractory particulate organic matter (RPOM), algae, and DO. Observational data recorded at the Incheon meteorological station of the Korean Meteorological Administration were used to calculate water temperatures [35]. The hourly temperature, wind direction, velocity, dew-point temperature, and cloudiness data were configured as input files. The simulation was conducted from January to December 2014 so that seasonal variation could be reflected. The minimum and maximum computational time steps were set to 1 and 3600 s, respectively. Within the assigned value, the average time step equaled 61 s in the range of 2-98 s. Table 1 presents the hydraulic and water quality parameters used for calibrating the CE-QUAL-W2 model.

Model Performance Assessment
The model performance was assessed using the root mean square error (RMSE) and absolute mean error (AME). These statistics are frequently used to examine the goodness of fit between modeling results and observed values. They were calculated using the following equations: where n is the number of measured and predicted pairs, O is an observed value, and P is a predicted value. Figure 2 shows the total inflow of GAW in 2014. The West Sea provided the largest annual inflow of 597.8 × 10 6 m 3 , followed by 247.2 and 16.3 × 10 6 m 3 from the Han River and Gulpo River, respectively. As for the relative proportions of annual inflow, the Han River, Gulpo River and West Sea accounted for 28.7%, 1.9%, and 69.4%, respectively. The Gulpo River is only diverted to the GAW navigational channel when Gyulhyeon Weir overflows because of heavy rainfall; this occurred seven times during the study period. The inflow from the West Sea showed a remarkable tidal fluctuation during the spring and neap tides. During the neap tides, inflow did not occur for 2 to 5 consecutive days twice a month. In addition, the inflow decreased from July to August due to flood control. In contrast, the Han River exhibited a relatively constant inflow rate throughout the year.    Table 2 presents the water quality parameters of the inflow water during the study period. As regards the annual average concentration of organic matter, the Gulpo River had a COD Mn of 12.4 mg/L, which far exceeded the 5.8 and 2.8 mg/L values for the Han River and West Sea. As regards the nutrient concentrations, the Gulpo River exhibited a TN of 12.196 mg/L and TP of 2.553 mg/L, indicating a strongly eutrophic state. The standard deviations were also large. Accordingly, because of the influence of these high concentrations of organic matter and nutrients, the DO in the Gulpo River was as low as 5.2 mg/L. The highest salinity level of 25.7 psu was recorded in the West Sea, while low salinity was observed for other rivers. The highest Chl-a level of 21.7 µg/L was recorded in the Han River. The Gulpo River, which had the poorest water quality, featured a lower level of Chl-a than the Han River because of dilution with water discharged from WWTPs. With the exception of the flood season, WWTP discharge accounted for almost 70% of the Gulpo River flow.

Salinity Stratification and DO Depletion
The salinity and DO were monitored in the vertical direction at each monitoring site. Various spatial distributions appeared to depend on the inflow and distance from each source. When there was considerable inflow from the West Sea, the salinity was high in the bottom layer and the strong salinity stratification was created ( Figure 3). The highest salinity was measured at S4, which was adjacent to the West Sea, while the lowest salinity was recorded at S1, which was at the entrance of the Han River. When inflow from the Gulpo River occurred, the salinity decreased dramatically. This phenomenon was observed in May, July, and August. The lowest salinity values of the year were obtained from July to early September, because of the inflow from the Gulpo River and the decrease in the inflow from the West Sea. From the middle of September, the inflow from the West Sea increased, such that the previous salinity and stratification levels were reestablished. other words, the hypoxia occurring in the upper regions from Gimpo T/M spread towards the Incheon T/M area when the current from the West Sea was less influential. The formation time of the hypoxic zone of GAW was similar to other esturine systems (e.g., Chesapeake Bay, Gulf of Mexico, and Tokyo Bay) owing to the increase in temperature and reduction of oxygen saturation in the summer [14][15][16][17][36][37][38]. In particular, there was a similarity in the spatial distribution of the hypoxic zone in the St. Lucie Estuary (USA) and Gironde Estuary (France), which are both affected by river discharge and have narrow width and low water depth. These estuaries were affected by a high concentration of nutrients entering the river, and thus, severe oxygen depletion occurred near the water zone adjacent to the river [11,13].  The spatial and seasonal variations of the DO concentration were more drastic than those of the salinity (Figure 4). DO was more abundant during seasons with lower water temperature (from January to April and from October to December). On the other hand, the DO concentration decreased dramatically from May to September. In May, the DO concentration began to decrease in the bottom layer rather than the surface layer. During summer, the vertical gradient of the DO concentration among the water layers was clear; it began to decrease in October. The most serious DO depletion was recorded at S1. However, greater proximity to the West Sea corresponded to higher DO concentration. At S1, a DO of 3 mg/L or lower persisted for approximately five consecutive months (from May to September) and hypoxia (<2 mg/L), occurred during the month of July. Hypoxia was also observed at S2, which was closest to S1. Lower frequencies of DO depletion were observed at S3 and S4. However, from July to August and in September, during which the inflow from the West Sea decreased, the DO concentration decreased dramatically because of the impact of the S1 hypoxia. In other words, the hypoxia occurring in the upper regions from Gimpo T/M spread towards the Incheon T/M area when the current from the West Sea was less influential. The formation time of the hypoxic zone of GAW was similar to other esturine systems (e.g., Chesapeake Bay, Gulf of Mexico, and Tokyo Bay) owing to the increase in temperature and reduction of oxygen saturation in the summer [14][15][16][17][36][37][38]. In particular, there was a similarity in the spatial distribution of the hypoxic zone in the St. Lucie Estuary (USA) and Gironde Estuary (France), which are both affected by river discharge and have narrow width and low water depth. These estuaries were affected by a high concentration of nutrients entering the river, and thus, severe oxygen depletion occurred near the water zone adjacent to the river [11,13].

Model Validation
The spatiotemporal variations of the salinity stratification, DO depletion, and other water quality parameters were simulated using CE-QUAL-W2. The observed water levels were compared with predicted values to determine whether the model could accurately represent the water level variation in the navigational channel under the tidal influence ( Figure 5). Although the predicted water levels were in relatively good agreement with those observed in practice, the model tended to overestimate the water level of GAW, particularly in the cases of low and neap tides, and RMSE/AME were determined as 0.166/0.130 m, respectively. Since the boundary conditions of the model were implemented at the external/internal flow rather than at the water level, the error of the observed flow affected water level calculations.

Model Validation
The spatiotemporal variations of the salinity stratification, DO depletion, and other water quality parameters were simulated using CE-QUAL-W2. The observed water levels were compared with predicted values to determine whether the model could accurately represent the water level variation in the navigational channel under the tidal influence ( Figure 5). Although the predicted water levels were in relatively good agreement with those observed in practice, the model tended to overestimate the water level of GAW, particularly in the cases of low and neap tides, and RMSE/AME were determined as 0.166/0.130 m, respectively. Since the boundary conditions of the model were implemented at the external/internal flow rather than at the water level, the error of the observed flow affected water level calculations.  Figure 6 illustrates the validation results for water quality constituents in the surface and bottom layers. The RMSEs of temperature, salinity, and DO were relatively low (1.9 °C, 2.7 psu, and 1.9 mg/L, respectively), i.e., the developed model realistically reproduced the water quality levels of the waterway. As for nutrient concentration, the RMSE and AME of TN levels were determined as 1.971 and 1.274 mg/L, respectively, while the corresponding values for TP were obtained as 0.130 and 0.085 mg/L, respectively. High concentrations were mainly underestimated by the model, but the overall variation was reproduced relatively well. Chl-a simulation results yielded an RMSE of 33.2 μg/L and an AME of 15.2 μg/L. Overestimation occurred at S1 and S2, while underestimation was detected at S3 and S4, which indicated an increase of the error range. The low accuracy of Chl-a level prediction was ascribed to model limitations. In particular, the model cannot distinguish between marine and freshwater algae species during Chl-a simulation. In GAW, freshwater and marine algae are affected by high and low salinity, respectively, but the used model does not consider salinity as a limiting factor for algal growth. These limitations affected the accuracy of the Chl-a prediction as well as the prediction accuracy of the DO concentration in the surface layer. Figures S1 and S2, as supplementary documents, provide time-series comparisons of observed and predicted values at the surface and in the bottom layer.
The simulation results for the vertical salinity and DO concentration are discussed with a focus on summer, when the DO depletion was serious (Figure 7). Validation results reveal that the vertical distributions of salinity and DO concentration were relatively well reproduced by the model. In particular, an RMSE of 0.9-3.7 psu and an AME of 0.9-3.7 psu were obtained for salinity modeling, i.e., the corresponding error range was small. The model also suitably simulated salinity stratification caused by density differences and the spatial distribution of concentration (differences in concentrations at S1, S2, S3, and S4 on the same day). For DO modeling, error ranges were determined as 0.2-3.2 mg/L (RMSE) and 0.1-2.8 mg/L (AME). Thus, the model was appropriately validated. In particular, summer hypoxia at S1 (Figure 7e) and DO depletion due to the spread of hypoxia to S4 (Figure 7h) were relatively well reproduced.  Figure 6 illustrates the validation results for water quality constituents in the surface and bottom layers. The RMSEs of temperature, salinity, and DO were relatively low (1.9 • C, 2.7 psu, and 1.9 mg/L, respectively), i.e., the developed model realistically reproduced the water quality levels of the waterway. As for nutrient concentration, the RMSE and AME of TN levels were determined as 1.971 and 1.274 mg/L, respectively, while the corresponding values for TP were obtained as 0.130 and 0.085 mg/L, respectively. High concentrations were mainly underestimated by the model, but the overall variation was reproduced relatively well. Chl-a simulation results yielded an RMSE of 33.2 µg/L and an AME of 15.2 µg/L. Overestimation occurred at S1 and S2, while underestimation was detected at S3 and S4, which indicated an increase of the error range. The low accuracy of Chl-a level prediction was ascribed to model limitations. In particular, the model cannot distinguish between marine and freshwater algae species during Chl-a simulation. In GAW, freshwater and marine algae are affected by high and low salinity, respectively, but the used model does not consider salinity as a limiting factor for algal growth. These limitations affected the accuracy of the Chl-a prediction as well as the prediction accuracy of the DO concentration in the surface layer. Figures S1 and S2, as supplementary documents, provide time-series comparisons of observed and predicted values at the surface and in the bottom layer.
The simulation results for the vertical salinity and DO concentration are discussed with a focus on summer, when the DO depletion was serious (Figure 7). Validation results reveal that the vertical distributions of salinity and DO concentration were relatively well reproduced by the model. In particular, an RMSE of 0.9-3.7 psu and an AME of 0.9-3.7 psu were obtained for salinity modeling, i.e., the corresponding error range was small. The model also suitably simulated salinity stratification caused by density differences and the spatial distribution of concentration (differences in concentrations at S1, S2, S3, and S4 on the same day). For DO modeling, error ranges were determined as 0.2-3.2 mg/L (RMSE) and 0.1-2.8 mg/L (AME). Thus, the model was appropriately validated.
In particular, summer hypoxia at S1 (Figure 7e) and DO depletion due to the spread of hypoxia to S4 (Figure 7h) were relatively well reproduced.

Physical Factors Affecting Hypoxia Formation
To analyze the physical factors affecting the hypoxia formation, the residence time and salinity distributions during the spring and neap tides were examined. In this study, the residence time was defined as the period during which water coming from the boundary remains in the waterbody and was calculated by setting the zero-order decay rate to −1 day −1 and zeroing out all other generic constituent kinetic parameters [27]. In each season, the time difference between a spring tide and a neap tide was between 5 and 9 days. However, there was a stark difference in the spatial distribution of the residence time (Figure 8). During the spring tide, when the inflow from the West Sea peaked, the residence time was shortest at S4 and longest at S1 (Figure 8a,c,e,g). In particular, the bottom layer of S1 exhibited the most stagnant characteristics. This waterbody with the longest residence time was mixed with the inflow from the Han River and spread to the upper layers of S2 and S3. On the other hand, during the neap tide, the West Sea provided no saltwater and the current from the Han River became stronger. Accordingly, the stagnant water mass in the bottom layer of S1 spread to the upper layer of S4 through S2 and S3 (Figure 8b,d,f,h). However, during summer, the bottom layer of S1 continuously exhibited a long residence time. This behavior seems to have occurred because the velocity of the stagnant water mass decreased as the discharge was reduced in accordance with the decreased inflow from the West Sea. Analysis revealed that this residence time spatial distribution created a physical environment in which the strongest hypoxia occurred at S1. Therefore, this simulation result indicates that the physical environment of the waterway undergoes a very drastic change in both time and space, and this feature greatly affects the water quality and ecosystem.

Physical Factors Affecting Hypoxia Formation
To analyze the physical factors affecting the hypoxia formation, the residence time and salinity distributions during the spring and neap tides were examined. In this study, the residence time was defined as the period during which water coming from the boundary remains in the waterbody and was calculated by setting the zero-order decay rate to −1 day −1 and zeroing out all other generic constituent kinetic parameters [27]. In each season, the time difference between a spring tide and a neap tide was between 5 and 9 days. However, there was a stark difference in the spatial distribution of the residence time ( Figure 8). During the spring tide, when the inflow from the West Sea peaked, the residence time was shortest at S4 and longest at S1 (Figure 8a,c,e,g). In particular, the bottom layer of S1 exhibited the most stagnant characteristics. This waterbody with the longest residence time was mixed with the inflow from the Han River and spread to the upper layers of S2 and S3. On the other hand, during the neap tide, the West Sea provided no saltwater and the current from the Han River became stronger. Accordingly, the stagnant water mass in the bottom layer of S1 spread to the upper layer of S4 through S2 and S3 (Figure 8b,d,f,h). However, during summer, the bottom layer of S1 continuously exhibited a long residence time. This behavior seems to have occurred because the velocity of the stagnant water mass decreased as the discharge was reduced in accordance with the decreased inflow from the West Sea. Analysis revealed that this residence time spatial distribution created a physical environment in which the strongest hypoxia occurred at S1. Therefore, this simulation result indicates that the physical environment of the waterway undergoes a very drastic change in both time and space, and this feature greatly affects the water quality and ecosystem. The variation in the salinity distribution was more strongly affected by the change of inflow during the spring and neap tides than by seasonal change (Figure 9). During the spring tide, the inflow from the West Sea increased and, thus, the salinity increased dramatically. For this reason, a saltwater wedge typical for an estuary occurred in the navigational channel. During the neap tide, saltwater tended to concentrate in the bottom layer, exhibiting slow movement. In general, the spring and neap tides in spring, autumn, and winter yielded similar salinity spatial distributions. On the other hand, during summer, as the inflow from the West Sea decreased, the salinity of the upper layer decreased to 10 psu or below, thereby causing freshwater dominance. Figure 10 presents the simulation result for the DO spatial distribution. It seems that the DO concentration distribution was affected by not only the above-mentioned physical factors such as residence time and salinity stratification, but also other complex factors such as water temperature and organic matter decomposition. The model accurately simulated the frequent DO depletion in the bottom layer of S1 and also revealed the spread of the hypoxia generated at S1 towards S4 during the neap tide. In particular, as shown in Figure 10c,d, the hypoxic waterbody in the bottom layer migrated to the middle layer, and the so-called metalimnetic oxygen minima were predicted. This The variation in the salinity distribution was more strongly affected by the change of inflow during the spring and neap tides than by seasonal change (Figure 9). During the spring tide, the inflow from the West Sea increased and, thus, the salinity increased dramatically. For this reason, a saltwater wedge typical for an estuary occurred in the navigational channel. During the neap tide, saltwater tended to concentrate in the bottom layer, exhibiting slow movement. In general, the spring and neap tides in spring, autumn, and winter yielded similar salinity spatial distributions. On the other hand, during summer, as the inflow from the West Sea decreased, the salinity of the upper layer decreased to 10 psu or below, thereby causing freshwater dominance. Figure 10 presents the simulation result for the DO spatial distribution. It seems that the DO concentration distribution was affected by not only the above-mentioned physical factors such as residence time and salinity stratification, but also other complex factors such as water temperature and organic matter decomposition. The model accurately simulated the frequent DO depletion in the bottom layer of S1 and also revealed the spread of the hypoxia generated at S1 towards S4 during the neap tide. In particular, as shown in Figure 10c,d, the hypoxic waterbody in the bottom layer migrated to the middle layer, and the so-called metalimnetic oxygen minima were predicted. This behavior was ascribed to the fact that the operation of the Gimpo terminal only allows freshwater inflow into the Han River while blocking discharge into the Han River, which results in the migration of the hypoxic waterbody mixed with fresh water to the middle layer. As a result, stagnant zones are generated in the waterbody, which yield hypoxic zones exhibiting DO depletion. This behavior is well illustrated by the time-series comparison of residence time and DO values at the bottom of each site ( Figure 11). Temporal changes in DO concentration and residence time clearly showed opposite trends in the bottom layer. Notably, these trends were significant at S1, S3, and S4 but relatively insignificant at S2, which was ascribed to the fact that at S2, the influence of the inflow of the Han River and the West Sea was offset by two-way flow. Consequently, the simulation results showed that GAW has the anthropogenic characteristics of estuarine water and, thus, a unique aquatic environment in which the physical conditions and water quality change dramatically during spring and neap tides. The inflow is artificially controlled, but is also affected by tides. In addition, there are two or three inflow sources, but discharge occurs in one direction only.
Studies on hypoxia formation in the Chesapeake Bay, the Gulf of Mexico, and the Tokyo Bay have examined the influence of wind on hypoxia in various ways [14][15][16][17][36][37][38][39][40]. The intensity and direction of the wind were evaluated as the main factors influencing the formation and duration of the hypoxic layer. In this study, however, the effect of wind on the stratification and depletion of DO in GAW could not be evaluated, and further study is needed. Further, many studies on the hypoxia formation in Lake Erie have quantitatively evaluated the effect of nutrient loading. Scavia et al. [39] and Rucinski et al. [40] suggested that the reduction of phosphorus may reduce the hypoxia area of the Lake Erie. Considering these prior studies, the reduction of nutrients from the Han River and the Gulpo River may have the potential to reduce hypoxia in GAW, and therefore, additional research is needed. behavior was ascribed to the fact that the operation of the Gimpo terminal only allows freshwater inflow into the Han River while blocking discharge into the Han River, which results in the migration of the hypoxic waterbody mixed with fresh water to the middle layer. As a result, stagnant zones are generated in the waterbody, which yield hypoxic zones exhibiting DO depletion. This behavior is well illustrated by the time-series comparison of residence time and DO values at the bottom of each site ( Figure 11). Temporal changes in DO concentration and residence time clearly showed opposite trends in the bottom layer. Notably, these trends were significant at S1, S3, and S4 but relatively insignificant at S2, which was ascribed to the fact that at S2, the influence of the inflow of the Han River and the West Sea was offset by two-way flow. Consequently, the simulation results showed that GAW has the anthropogenic characteristics of estuarine water and, thus, a unique aquatic environment in which the physical conditions and water quality change dramatically during spring and neap tides. The inflow is artificially controlled, but is also affected by tides. In addition, there are two or three inflow sources, but discharge occurs in one direction only. Studies on hypoxia formation in the Chesapeake Bay, the Gulf of Mexico, and the Tokyo Bay have examined the influence of wind on hypoxia in various ways [14][15][16][17][36][37][38][39][40]. The intensity and direction of the wind were evaluated as the main factors influencing the formation and duration of the hypoxic layer. In this study, however, the effect of wind on the stratification and depletion of DO in GAW could not be evaluated, and further study is needed. Further, many studies on the hypoxia formation in Lake Erie have quantitatively evaluated the effect of nutrient loading. Scavia et al. [39] and Rucinski et al. [40] suggested that the reduction of phosphorus may reduce the hypoxia area of the Lake Erie. Considering these prior studies, the reduction of nutrients from the Han River and the Gulpo River may have the potential to reduce hypoxia in GAW, and therefore, additional research is needed.

Assessment of Fish Habitat Volume
The fish numbers surveyed in each season of the study period are presented in Appendix A [26]. In this study, three main fish species were selected, and their habitat environments were investigated in order to determine habitat volumes in relation to temperature and DO (Table 3). As the habitat volumes, the minimum habitat conditions for survival were considered instead of the optimal conditions. Except for the crucian carp (CRC), the other three species were selected because they were frequently observed at every site (S1-S4). In contrast, CRC appeared at S1 only. However, many CRCs live in the Han and Gulpo Rivers and are likely to enter the waterway with freshwater flows. Therefore, CRC was included in the main fish species, since its habitat encompasses regions with broadly variable temperature. Their habitats also have very low DO limit concentrations. Accordingly, these species can endure hypoxia well. The dotted gizzard shad (DGS) also has a wide habitat range, while the flathead mullet (FHM) was observed to have a narrow habitat range and has a low tolerance for DO depletion. Table 3. Habitat temperature and DO criteria for three major fish species in GAW.  Figure 12 presents the fish habitat volumes for each fish species calculated using the model, revealing that CRC could survive in every volume irrespective of season. However, the DGS habitat volume decreased between January and February and in December, when the water temperature decreased dramatically, and between July and August, when strong hypoxia occurred. FHM was most sensitive to environmental change. This species had few habitat volumes for survival between January and March and in December, when the water temperature was low. In late June, its habitat volume was reduced to 77% of the total water volume as a result of hypoxia. Furthermore, the average habitat volumes for survival were restricted to 34% and 32% in July and between August and September, respectively.

Conclusions
This study analyzed the effect of hydraulic and water quality factors on the spatiotemporal formation of summer hypoxia in GAW, an artificial estuarine waterway, by applying the twodimensional CE-QAUL-W2 hydrodynamic and water quality model. In addition, the ecological vulnerability of major fish species to hypoxia was determined, by calculating habitat volumes in relation to temperature and DO.
The salinity and DO of GAW were monitored throughout the year 2014, in both the horizontal and vertical directions. Various spatial distributions were obtained according to inflow and distance Analysis of the fish habitat volumes yielded by the simulation revealed the following: if a fish species sensitive to hypoxia enters the waterway, a poor habitat environment could be maintained for a long period of time. Recall also that the habitat volumes calculated in this study correspond to the minimum habitats for survival, and the optimal habitat volumes are considerably more limited. Besides, as the calculations did not distinguish between life stages, the survival behavior for the larvae and juvenile stages would be more limited.
The modeling results of this study revealed that GAW with its artificial estuarine environment produces hypoxia due to continuous eutrophic inflow and slow physical circulation. Moreover, as the waterway exhibits very complex waterbody behaviors, simulation-based prediction of spatiotemporal production and the spread of hypoxia seems essential. In order to maintain a sound aquatic environment, fish habitats must be assessed continuously. It is expected that the model employed in this study will be effectively utilized to derive operational measures for GAW, which could increase habitat volumes.

Conclusions
This study analyzed the effect of hydraulic and water quality factors on the spatiotemporal formation of summer hypoxia in GAW, an artificial estuarine waterway, by applying the two-dimensional CE-QAUL-W2 hydrodynamic and water quality model. In addition, the ecological vulnerability of major fish species to hypoxia was determined, by calculating habitat volumes in relation to temperature and DO.
The salinity and DO of GAW were monitored throughout the year 2014, in both the horizontal and vertical directions. Various spatial distributions were obtained according to inflow and distance from each inflow source. The salinity level was highest at the Gimpo connecting terminal (S4) and lowest at the Incheon connecting terminal (S1). When the inflow from the West Sea increased, the salinity stratification increased. DO was more abundant when the water temperature was low. However, during summer, the concentration gradient among the water layers became clearer and hypoxia was observed. The longest period of hypoxia was observed at S1. When the inflow from the West Sea decreased, the hypoxia at S1 seemed to be spread towards S4.
The simulation results showed that the CE-QAUL-W2 model realistically reproduced the water quality characteristics of GAW. Both the salinity and DO were relatively well simulated with error rates of RMSE 2.7 psu and 1.9 mg/L, respectively. The simulation results for the vertical salinity and DO also indicated that the salinity stratification, vertical DO gradient, and summer hypoxia were very accurately reproduced.
To analyze the physical factors affecting the hypoxia generation, the residence time, salinity, and DO concentration distributions during the spring and neap tides were comprehensively analyzed for each season. The bottom layer of S1 exhibited the longest residence time, which seemed to be the greatest physical factor influencing hypoxia. This waterbody with limited circulation spread towards S4 during summer and the neap tide, when the inflow from the West Sea decreased. These simulation results provided a basis for identifying the drivers of the observed DO depletion and enabled satisfactory explanation of the generation and distribution of hypoxic zones. The simulation results also clarified that GAW has a unique aquatic environment, with physical conditions and water quality levels that change dramatically according to the tidal flow (inflow and outflow) regulation.
Finally, the fish habitat volumes for survival were calculated. The flathead mullet was found to be most sensitive to environmental change. In summer, the habitat volume of this species was limited to 32-34% of the total volume as a result of hypoxia. Note that, if life stages were considered, the survival habitat volume of this species would be further reduced; therefore, hypoxia countermeasures are required. In order to maintain a sound aquatic environment, fish habitats must be continuously assessed in further studies. The model employed in this study can be effectively utilized to derive operational measures for GAW, with the potential to increase habitat volumes.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4441/10/11/1695/s1, Figure S1: Comparisons of predicted and observed water quality constituents in the surface layer, Figure S2: Comparisons of predicted and observed water quality constituents in the bottom layer.
Turbulent kinetic energy (TKE) model where U and W are the laterally averaged velocity components (m/s) in the x and z directions, respectively; B is the width of the water body (m); t is time (s); g is the gravitational acceleration (m/s 2 ); α is an angle and tan α is defined as the channel slope; ρ is density (kg/m 3 ); P is pressure (N/m 2 ); τ xx and τ xz are the turbulent shear stresses acting in the x direction on the x and z faces of the control volume (N/m 2 ), respectively; A x and A z are the longitudinal and turbulent eddy viscosities (m 2 /s); η is the free water surface location (m); q is the lateral boundary inflow or outflow per cell volume (1/s); B η is the time and spatially varying surface width (m); h is the total depth (m); f (T w , φ TDS , φ SS ) is the density function dependent on water temperature, total dissolved solids or salinity, and suspended solids; φ is the laterally averaged constituent concentration (g/m 3 ); D x and D z are the temperature and constituent dispersion coefficients in the x and z directions, respectively; q φ is the lateral boundary inflow or outflow mass flow rate of the constituent (g/m 3 s); S φ is the source/sink rate term for the constituent concentrations (g/m 3 s); C µ is an empirical constant; k is the turbulent kinetic energy (in the TKE model); ε is the turbulent energy dissipation rate; P 2 is the turbulent energy production from boundary friction; G is a buoyancy term; P k , P ε are production terms from boundary friction; C is a Chezy friction factor; n is a Manning friction factor; R h is the hydraulic radius; N is the Brunt-Vaisala frequency; σ is the turbulent Prandtl number; and C µ and C ε are constants in the TKE model. Typical values of the empirical constants in the above model are σ k = 1.0, σ ε = 1.3, σ t = 1.0, C µ = 0.09, C 1ε = 1.44, and C 2ε = 1.92.