A numerical model study of the main factors contributing to hypoxia and its interannual and short-term variability in the East China Sea

A three-dimensional physical-biological model of the marginal seas of China was used to analyze interannual and intra-seasonal variations in hypoxic conditions and identify the main processes controlling their generation off the Changjiang (or Yangtze River) estuary. The model was compared against available observations and reproduces the observed temporal and spatial variability of physical and biological properties including bottom oxygen. Interannual variations of hypoxic extent in the simulation are partly explained by variations in river discharge but not nutrient load. As riverine inputs of freshwater and nutrients are consistently high, promoting large productivity and subsequent oxygen consumption in the region affected by the river plume, wind forcing is important in modulating interannual and shortterm variability. Wind direction is relevant because it determines the spatial extent and distribution of the freshwater plume, which is strongly affected by either upwelling or downwelling conditions. High-wind events can lead to partial reoxygenation of bottom waters and, when occurring in succession throughout the hypoxic season, can effectively suppress the development of hypoxic conditions, thus influencing interannual variability. A model-derived oxygen budget is presented and suggests that sediment oxygen consumption is the dominant oxygen sink below the pycnocline and that advection of oxygen in the bottom waters acts as an oxygen sink in spring but becomes a source during hypoxic conditions in summer, especially in the southern part of the hypoxic region, which is influenced by open-ocean intrusions.

able observations and reproduces the observed temporal and spatial variability of physical and biological properties including bottom oxygen. Interannual variations of hypoxic extent in the simulation are partly explained by variations in river discharge but not nutrient load. As riverine inputs 10 of freshwater and nutrients are consistently high, promoting large productivity and subsequent oxygen consumption in the region affected by the river plume, wind forcing is important in modulating interannual and short-term variability. Wind direction is relevant because it determines the spa-15 tial extent and distribution of the freshwater plume which is strongly affected by either upwelling or downwelling conditions. High-wind events can lead to partial reoxygenation of bottom waters and, when occurring in succession throughout the hypoxic season, can effectively suppress the devel-20 opment of hypoxic conditions thus influencing interannual variability. A model-derived oxygen budget is presented and suggests that sediment oxygen consumption is the dominant oxygen sink below the pycnocline and that advection of oxygen in the bottom waters acts as an oxygen sink in spring but 25 becomes a source during hypoxic conditions in summer, especially in the southern part of the hypoxic region, which is influenced by open-ocean intrusions.

Introduction
In coastal seas, hypoxic conditions (oxygen concentrations 30 lower than 2 mg L −1 or 62.5 mmol m −3 ) are increasingly caused by rising anthropogenic nutrient loads from land (Diaz and Rosenberg, 2008;Rabalais et al., 2010;Fennel and Testa, 2019). Hypoxic conditions are detrimental to coastal ecosystems leading to a decrease in species diversity and ren-35 dering these systems less resilient (Baird et al., 2004;Bishop et al., 2006;Wu, 2002). Hypoxia is especially prevalent in coastal systems influenced by major rivers such as the northern Gulf of Mexico (Bianchi et al., 2010), Chesapeake Bay , and the Changjiang Estuary (CE) in the East 40 China Sea (Li et al., 2002).
The Changjiang is the largest river in China and fifth largest in the world in terms of volume transport, with an annual discharge of 9 × 10 11 m 3 year −1 via its estuary (Liu et al., 2003). The mouth of the CE is at the confluence of the 45 southeastward Yellow Sea Coastal Current and the northward Taiwan Warm Current (Figure 1). Hydrographic properties in the outflow region of the CE are influenced by several different water masses including fresh Changjiang Diluted Water, relatively low-salinity coastal water, more saline wa-50 ter from the Taiwan Warm Current, and high-nutrient, lowoxygen water from the subsurface of the Kuroshio (Wei et al., 2015;Yuan et al., 2008). The interactions of these water masses together with wind forcing and tidal effects lead to a complicated and dynamic environment. 55 Freshwater (FW) discharge by the Changjiang reaches its minimum in winter when the strong northerly monsoon (dry season) prevails and peaks in summer during the weak southerly monsoon (wet season) resulting in a large FW plume adjacent to the estuary. Along with the FW, the 60 Changjiang delivers large quantities of nutrients to the East China Sea (ECS) resulting in eutrophication in the plume region (Li et al., 2014;Wang et al., 2016). Since the 1970s, nutrient load has increased more than twofold with a subsequent increase in primary production (PP) in the outflow 5 region of the estuary . Hypoxia off the CE was first detected in 1959 and, with a spatial extent of up to 15,000 km 2 , is among the largest coastal hypoxic zones in the world (Fennel and Testa, 2019). Although no conclusive trend in oxygen minima has been observed (Wang, 2009;10 Zhu et al., 2011), hypoxic conditions are suspected to have expanded and intensified in recent decades (Li et al., 2011;Ning et al., 2011) due to the increasing nutrient loads from the Changjiang .
It is generally accepted that water-column stratification 15 and the decomposition of organic matter are the two essential factors for hypoxia generation, and this is also the case for the shelf region off the CE (Chen et al., 2007;Li et al., 2002;Wei et al., 2007). High solar radiation and FW input in summer contribute to strong vertical stratification which 20 is further enhanced by near-bottom advection of waters with high salinities (> 34) and low temperatures (<19 • C) by the Taiwan Warm Current. The resulting strong stratification inhibits vertical oxygen supply (Li et al., 2002;Wang, 2009;Wei et al., 2007). At the same time, high organic matter sup- 25 ply fuels microbial oxygen consumption in the subsurface (Li et al., 2002;Wang, 2009;Wei et al., 2007;Zhu et al., 2011). It has also been suggested that the Taiwan Warm Current brings additional nutrients contributing to organic matter production (Ning et al., 2011) and that the low oxygen 30 concentrations (∼5 mg L −1 ) of the Taiwan Warm Current precondition the region to hypoxia (Ning et al., 2011;Wang, 2009). While observational analyses suggest that hypoxia off the CE results from the interaction of various physical and bio-35 geochemical processes, quantifying the relative importance of these processes and revealing the dynamic mechanisms underlying hypoxia development and variability require numerical modeling (Peña et al., 2010). Numerical modeling studies have proven useful for many other coastal hypoxic 40 regions such as the Black Sea northwestern shelf (Capet et al., 2013), Chesapeake Bay Scully, 2013), and the northern Gulf of Mexico (Fennel et al., 2013;Laurent and Fennel, 2014). Models have also been used to study the hypoxic region of the CE. Chen et al. (2015a) used a 3D circulation model with a highly simplified oxygen consumption parameterization (a constant consumption rate) to investigate the effects of physical processes, i.e. FW discharge, and wind speed and direction, on the dissipation of hypoxia. Chen et al. (2015b) examined the tidal modulation of hy-50 poxia. The model domain in these two previous studies was relatively limited encompassing only the CE, Hangzhou Bay and the adjacent coastal ocean but did not cover the whole area affected by hypoxia (Wang, 2009;Zhu et al., 2011). Zheng et al. (2016) employed a nitrogen cycle model cou- 55 pled with a 3D hydrodynamic model to examine the role of river discharge, wind speed and direction on hypoxia, and also emphasized the physical controls. These previous modeling studies focused on the response of hypoxia to physical factors only and did not address seasonal evolution and inter-60 annual variations of hypoxia or the influence of variability in biological rates.
More recently, Zhou et al. (2017) analyzed the seasonal evolution of hypoxia and the importance of the Taiwan Warm Current and Kuroshio intrusions as a nutrient source using 65 an advanced coupled hydrodynamic-biological model. However, the baseline of their model does not include sediment oxygen consumption (SOC), which is thought to be a major oxygen sink in the hypoxic region off the CE (Zhang et al., 2017) and other river-dominated hypoxic regions including 70 the northern Gulf of Mexico (Fennel et al. 2013, Yu et al. 2015a. Zhou et al. (2017) acknowledged the importance of SOC based on results from a sensitivity experiment but did not quantify its role in hypoxia generation.
Here we introduce a new 3D physical-biological model 75 implementation for the ECS that explicitly includes nitrogen and phosphorus cycling and SOC. The model is a new regional implementation for the ECS of an existing physicalbiogeochemical model framework that has been extensively used and validated for the northern Gulf of Mexico (Fennel 80 et al., 2011(Fennel 80 et al., , 2013Laurent et al., 2012;Laurent and Fennel, 2014;Yu et al., 2015b;Fennel and Laurent, 2018). The hypoxic zones in the northern Gulf of Mexico and off the CE have similar features including the dominant influence of a major river (Changjiang and Mississippi), a seasonal recur-85 rence every summer, a typical maximum size of about 15,000 km 2 , documented P-limitation following the major annual discharge in spring and a significant contribution of SOC to oxygen sinks in the hypoxic zone (Fennel and Testa 2019).
We performed and assessed a 6-year simulation of the 90 ECS and use the model results here to identify the main factors driving hypoxia variability on interannual and short-term (days to seasons) timescales in the simulation. More specifically, we investigate the role of interannual variations in riverine inputs of nutrients and FW versus short-term vari-95 ations in coastal circulation and mixing. We also present an oxygen budget to quantify the relative importance of SOC and the influence of lateral advection of oxygen in the model. The companion study by Grosse et al. (2020) used the same model to quantify the importance of intrusions of nutrient-100 rich oceanic water from the Kuroshio for hypoxia development off the CE. The physical model is based on the Regional Ocean Mod-105 eling System (ROMS; Haidvogel et al., 2008) and was im-plemented for the ECS by Bian et al. (2013a). The model domain extends from 116 • E to 134 • E and 20 • N to 42 • N (Figure 1), covering the Bohai Sea, the Yellow Sea, the ECS, part of the Japan Sea and the adjacent northwest Pacific, with a horizontal resolution of 1/12 • (about 10 km) and 5 30 vertical layers with enhanced resolution near the surface and bottom. The model uses the recursive Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) for the advection of tracers (Smolarkiewicz and Margolin, 1998), a third-order upstream advection of momentum, and the Generic Length Scale (GLS) turbulence closure scheme (Umlauf and Burchard, 2003) for vertical mixing.
The model is initialized with climatological temperature and salinity from the World Ocean Atlas 2013 (WOA13) (Locarnini et al., 2013;Zweng et al., 2013), and is forced 15 by 6-hourly wind stress, and heat and FW fluxes from the ECMWF ERA-Interim dataset (Dee et al., 2011). Open boundary conditions for temperature and salinity are also from WOA13, and horizontal velocities and sea surface elevation at the boundaries are specified from the monthly 20 SODA data set (Carton and Giese, 2008). In addition, eight tidal constituents (M2, S2, N2, K2, K1, O1, P1 and Q1) are imposed based on tidal elevations and currents extracted from the global inverse tide model TPXO7.2 of Oregon State University (OSU, Egbert and Erofeeva, 2002). At the 25 open boundaries, Chapman and Flather conditions are used for free surface and barotropic velocities, respectively, and the radiation condition for the baroclinic velocity. Eleven rivers are included in the model. FW discharge from the Changjiang uses daily observations from the Datong Hydro-30 logical Station. For the other rivers, we prescribed monthly or annual climatologies (Liu et al., 2009;Tong et al., 2015;Zhang, 1996).

Biological model
The biological component is based on the pelagic nitro-35 gen cycle model of Fennel et al. (2006Fennel et al. ( , 2011Fennel et al. ( , 2013 and was extended to include phosphate (Laurent et al., 2012;Laurent and Fennel, 2014) and riverine dissolved organic matter (Yu et al., 2015b). The model includes two forms of dissolved inorganic nitrogen (DIN), nitrate (NO3) and 40 ammonium (NH4), phosphate (PO4), phytoplankton (Phy), chlorophyll (Chl), zooplankton (Zoo), two pools of detritus, suspended and slow-sinking small detritus (SDet) and fastsinking large detritus (LDet), and riverine dissolved organic matter (RDOM). Here, riverine dissolved and particulate organic nitrogen enter the pools of RDOM and SDet, respectively. The remineralization rate of RDOM is an order of magnitude lower than that of SDet to account for the more refractory nature of the riverine dissolved organic matter (Yu et al., 2015b).

50
At the sediment-water interface, SOC is parameterized assuming "instantaneous remineralization," i.e. all organic matter reaching the sediment is remineralized instantaneously and oxygen is consumed due to nitrification and aerobic remineralization at the same time. In the "instantaneous reminer-55 alization," all phosphorus is returned to the water column as PO4 while a constant fraction of fixed nitrogen is lost due to denitrification. All biogeochemical model parameters are given in Table S1 in the Supplement. A more detailed model descriptions can be found in the Supplement to Laurent et al. 60 (2017). Light is vertically attenuated by chlorophyll, detritus and seawater itself. In addition, to account for the effects of colored dissolved organic matter (CDOM) and suspended sediments, which show relatively high values near the coast and in the river plume (Bian et al., 2013b;Chen et al., 2014), 65 a light-attenuation term dependent on water depth and salinity is introduced which yields higher attenuation in shallow areas and in the FW plume. Initial and boundary conditions for NO3, PO4 and oxygen are prescribed using WOA13 (Garcia et al., 2013a,b). 70 A small positive value is used for the other variables. NO3 is nudged towards climatology in the northwest Pacific at depth > 200 m. Monthly nutrient loads of NO3 and PO4 from the Changjiang are from the Global-NEWs Model  but were adjusted by multiplicative factors of 1.20 and 1.66, respectively, to ensure a match between simulated and observed nutrient concentrations in the CE (see July and Aug 2012 in Figure 2). Nutrient loads in other rivers are based 5 on other published climatologies (Liu et al., 2009;Tong et al., 2015;Zhang, 1996). Due to a lack of data on organic matter loads, river load concentrations of SDet and LDet and RDOM were assumed conservatively at 0.5, 0.2 and 15 mmol N m −3 , respectively.

30
The simulated current systems in the ECS and YS show typical seasonal variations as follows (see also Figure S6). In winter, currents mainly flow southward on the Yellow Sea and ECS shelves driven by the northerly wind. In contrast, the ECS Coastal Current and the Korean Coastal Current 35 flow northward in summer. The Kuroshio is stronger in summer than in winter. The model captures the seasonal pattern of the current system and resolves currents in the ECS and Yellow Sea (also see Grosse et al. 2020).
Simulated monthly averaged (2008-2013) surface chloro-40 phyll concentrations in May, August and November are compared with satellite-derived fields and agree well with spatial correlation coefficients of 0.77, 0.94 and 0.64, respectively ( Figure   Together, these comparisons show that the model is able to reproduce important aspects of the physical-biogeochemical dynamics in the study region.

Simulated oxygen dynamics
First, we describe the timing and distribution of simulated bottom-water oxygen off the CE to set the stage for our investigation into the drivers underlying hypoxia variability. The model simulates annually recurring hypoxic condi- 125 mmol m −3 (or 4 mg L −1 ), i.e. twice the hypoxic threshold. It is known from observations that there are two centers of recurring hypoxic conditions: the northern core is located just to the east of the CE and Hangzhou Bay and the southern core to the southeast of Hangzhou Bay. The model is 30 consistent with these observations and simulates two distinct core regions of low-oxygen conditions centered at 31 • N and 29.3 • N. The northern core region is larger than the southern core region (9,050 km 2 for a threshold of 80 days per year of <4 mg L −1 compared to 5,230 km 2 ). We will refer to the 35 region defined by a threshold of 40 days of <4 mg L −1 of per year (solid black line in Figure 1 and 4d) as the "typical low-oxygen zone" for the remainder of the manuscript and demarcate the northern and southern sections by 30.1 • N latitude (dashed line in Figures 1 and 4d). 40 There are marked differences in the phenology of simulated hypoxic extent (Figure 4c). Among the four years with largest hypoxic areas, hypoxia establishes relatively late (mid-August) and lasts long (into November) in 2008 and 2009. In contrast in 2012, hypoxic conditions establish ear-45 lier (June), are most pronounced in August and are eroded by mid-October. In 2010, the year with the largest peak extent, hypoxia establishes already at the beginning of June and is maintained until the end of October, leading to the largest time-integrated hypoxia by far among the 6 years (Figure 50 4b). In all years there are times when hypoxic extent decreases rapidly.
In the following sections, we explore the drivers underlying these interannual and short-term variations, specifically the contribution of year-to-year variations in nutrient loads 55 and FW inputs from the Changjiang, and the potential reasons for shorter-term variability in hypoxia by assessing the role of biological processes and physical forcing.

Interannual variations in hypoxia
The first question we address is: Do year-to-year variations 60 in nutrient load and FW input from the Changjiang explain interannual variability in hypoxic conditions? We do this by investigating correlations of time-integrated hypoxic area, average PP, total oxygen consumption (OC) by respiration, SOC, and bottom oxygen in the typical low-oxygen zone 65 ( Figure 5 a-f). We also consider the correlation between the spatial extent of the FW plume, defined as the horizontal ex- There is a significant negative correlation between annual FW input and mean bottom-water oxygen concentration in the low-oxygen zone of -0.86 and a weaker, statistically insignificant positive correlation of 0.69 between annual FW input and integrated hypoxic area (Figure 5a, d). This indicates that variations in FW input at least partly explain variability in hypoxic conditions. Perhaps surprisingly, there is no convincing correlation between annual FW input and annual DIN load (Figure 5h). Although the correlation coefficient is 0.56 when all 6 years are considered, the correlation reverses to -0.17 when the low-flow year 2011 is excluded and neither of these correlations is statistically significant. As expected, there is a strong positive correlation of 0.84 be-15 tween the annual FW input and time-integrated plume area ( Figure 5g). Plume area can thus be interpreted as a proxy of FW input.
In contrast to the positive correlations between FW input and hypoxia, and FW input and bottom oxygen, correlations 20 between the annual DIN load with integrated hypoxic area and mean bottom-water oxygen are much weaker and insignificant (Figure 5b, e). This implies that interannual variations in DIN load do not lead to year-to-year variations in hypoxia. However, the correlations between integrated hypoxic 25 area and mean rates of PP and OC (especially SOC) in August are significant and strong at 0.94 and 0.93 (0.97), respec-tively (Figure 5c, f, i). The high correlation between hypoxic area and OC is primarily driven by SOC. Clearly, biological processes are important drivers of hypoxia and contribute 30 to its interannual variability, but they do not appear to result from variations in DIN load. More relevant are variations in FW load, which explain interannual variations in hypoxia at least partly.
Other factors than riverine inputs of nutrients and FW must 35 be contributing to interannual variations. For example, the years 2010 and 2012 both had very similar FW input and DIN load but differed in severity of hypoxia (Figure 5a, b). Likewise, the years 2009 and 2013 were very similar in terms of FW input and DIN load, but very different in hypoxic ex-40 tent. Next, we investigate the potential reasons for shorterterm variability in hypoxia, i.e. the processes leading to the differences in hypoxia phenology in Figure 4c.

Biological drivers of short-term variability in hypoxia 45
In the previous subsection, we identified biological rates as important drivers of low-oxygen conditions on interannual timescales but unrelated to variations in riverine DIN load.
Here we attempt to elucidate what drives variations in biological rates and low-oxygen conditions on shorter timescales 50 by addressing the following two questions. Do low-oxygen conditions correlate with biological rates on these shorter timescales? If yes, what drives variations in biological rates? For this analysis it seems prudent to distinguish between the northern and southern hypoxic regions for the following reasons. The bathymetry in the northern zone is slightly deeper than in the southern zone (median depth of 28.5 m versus 24.6 m) and several biological rates with di- In the water column, the difference is reversed and WOC larger in the southern than the northern zone (52.9 versus 46.7 mmol O 2 m −2 d −1 ). Because of these different characteristics, we consider the northern and southern zones of the typical low-oxygen re-20 gion separately.
First, we explore whether significant relationships exist between daily biological rates and bottom-water oxygen by determining the correlations of daily averaged rates of PP, OC and SOC with daily mean bottom oxygen concentration 25 (Figure 7 and Table 1).
Indeed, daily PP, OC, and SOC are all significantly and negatively correlated with bottom-water oxygen. This confirms that local production of organic matter and the resulting biological oxygen consumption are important for hy-30 poxia development and that variations in these rates partly explain variations in low-oxygen conditions. However, it is also obvious that variability around the best fit is large (Figure 7). The next question is: What drives variations in the biological rates? Since the annual correlations presented in 35 the previous section indicate that variability in annual FW input partly explains interannual variability in hypoxia, we consider whether FW variability is related to variations in biological rates. Using daily plume extent as a measure of FW  Table 1. presence and comparing it to daily rates of PP, OC, SOC, and bottom oxygen, we find that bottom oxygen and biological rates are significantly correlated with the extent of the FW plume with correlation coefficients ranging from 0.43 to 0.62 (Table 1). In other words, variability in the extent of the FW 5 plume explains roughly half of the variability in biological rates. Mechanistically, the presence of a large FW plume not only affects hypoxia by increasing vertical stratification and thus inhibiting vertical supply of oxygen to the subsurface but also because PP and respiration is larger in the plume.

10
Large FW plumes stimulate more widespread biological production and thus oxygen consumption. Since annual FW input is highly correlated with the extent of the FW plume (see Figure 5g), variability in its extent is partly due to variations in riverine input, but coastal 15 circulation and mixing processes must be playing a role as well. Next, we analyze the impact of the underlying physical drivers.

Physical drivers of short-term variability in hypoxia
We focus our analysis of physical drivers on wind direction and wind strength, and their relation to FW plume location and extent because the latter has already been identified as 5 an explanatory variable for interannual variations in the previous section. Wind direction is relevant because for most of June, July, and August winds blow predominantly from the south, but switch to predominantly northerly winds between the sec-10 ond half of August and the end of September. As a result of the northward, upwelling favorable winds in the early summer, the FW plume is spread offshore and overlaps primarily with the northern zone. After the switch to mostly southward, downwelling-favorable directions, the FW plume 15 moves southward, becomes more contained near the coast, and grows in its southward extent as it is transported by a coastal current. Wind direction has a demonstrable impact on PP and the extent of the FW plume as shown in Figure  8 for the month of September. Especially in the northern re-20 gion, PP and plume extent are notably larger during southerly winds when the FW plume is more spread out, than during northerly winds when the plume is more restricted within the coastal current.
Wind strength is relevant because storm events can erode 25 vertical stratification and thus lead to resupply of oxygen to bottom waters due to vertical mixing. We investigated the effect of wind strength on bottom oxygen, hypoxia, and the extent of the FW plume by first inspecting time series of these variables ( Figure S8). We isolated all events during the 30 months June to September and, in Figure 10, show the corresponding changes in wind stress, mean bottom oxygen in the northern and southern zones, and the extent of the FW plume. We diagnosed these events as follows. First, we identified all days when the wind stress exceeded 0.12 Pa. Then 3 days prior and 3 days after the high-wind days. The periods within these minima are used as analysis period for each wind event. In four instances the wind stress exceeded 40 the threshold within 5 days of a previous wind event. Those subsequent high-wind events were combined into one. We identified the minimum in bottom oxygen (maximum in FW plume area) at the beginning of the event and the maximum in oxygen (minimum in FW area) after the maximum in wind 45 stress was reached. Figure 9a illustrates rapid increases in wind stress typically within 2 to 4 days. The only exceptions are the 4 events where two storms occurred in rapid succession and the combined event lasted longer (up to 8 days) until maxi-50 mum wind stress was reached. The year with the most wind events is 2013 (with 8 in total including one of the combined long-lasting event). The year with the least events is 2010 (2 events) followed by 2009 (3 events). Most of these events resulted in notable increases in mean bottom oxygen, typically 55 by 10 to 30 mmol m −3 , but up to 100 mmol m −3 in 2010 in the southern zone (Figure 9b). In the rare cases where bottom oxygen did not increase or slightly decreased, bottom oxygen was already elevated before the wind event. The wind events strongly affected the extent of the FW plume (Figure 9c) by mixing the FW layer with underlying ocean water. The effects were largest when the FW plume was most expansive.

5
This analysis shows the significant role of storm events in disrupting the generation of low-oxygen conditions and ventilating bottom waters.
In section 3.2.1 above, where we discussed interannual variability, we noted that while the years 2010 and 2012 had 10 very similar FW input and DIN load, 2010 had a much larger hypoxic area. Likewise, the years 2009 and 2013 were very similar in terms of FW input and DIN load, but 2009 had a much larger hypoxic area. It now becomes obvious that the frequency and severity of high-wind events, i.e. variations on 15 short timescales, explains the interannual differences in both cases. Figure 10 shows the wind stress, mean bottom oxygen in the northern and southern zones, and total hypoxic extent and FW plume extent in 2012 and 2010. In 2012, there were 5 20 high-wind events during the months of August, September, and October that all coincided with increases in bottom oxygen, decreases in hypoxic extent when a hypoxic zone was established at the beginning of the event, and decreases in FW plume extent. Inspection of the evolution of bottom oxy- 25 gen is especially instructive. While bottom oxygen concentrations declined during periods with average or low wind, they were essentially reset at a much higher level during each wind event. Whenever the FW plume was extensive at the beginning of a high-wind event, it was drastically reduced dur-30 ing the event. In 2010, bottom oxygen was at similar levels to 2012 at the beginning of August but dropped to low levels throughout August, especially in the northern zone, and remained low with widespread hypoxia until a major wind event in the second half of October ventilated bottom waters. 35 Except for a very short event in the second half of September, there were no high-wind events from August until mid-October in 2010.
The differences in hypoxia in 2009 and 2013 can also be explained by the frequency and intensity of high-wind 40 events. In 2013, there were 8 high-wind events from July to October that led to an almost continuous ventilation of bottom waters while in 2009 there were only 3 such events during the same period ( Figure S8). Low to average winds from mid-August to early October of 2009 coincided with a 45 decline in bottom oxygen and establishment of an expansive hypoxic zone throughout most of September.
These analyses show that wind direction and strength play an important role in determining the location of the hypoxic zone (i.e. northern versus southern region) and the extent and 50 severity of hypoxic conditions.

Oxygen budgets for the northern and southern regions
In order to further investigate the roles of physical and biological processes in regulating hypoxia, oxygen budgets 55 were calculated from daily model output for the period from March to November for the northern and southern hypoxic regions. Considering that hypoxic conditions occur near the bottom, we evaluate an oxygen budget not only for the whole water column but also for its lower portion which typically 60 becomes hypoxic. To account for variations in the thickness of the hypoxic layer, which tends to be thicker in deeper waters (similar to observations by Ning et al., 2011), we include the bottommost 12 layers of our model grid. Because of the model's terrain-following vertical coordinates, the thick-65 ness of these 12 model layers varies with total depth. The terms considered in the budget are air-sea flux, lateral physical advection and diffusion, vertical turbulent diffusion (for the subsurface budget only), PP, WOC (including respiration and nitrification), and SOC. Each term was integrated verti-70 cally over the whole water column and also over the bottommost 12 layers and then averaged for the northern and southern regions for each month (Figure 11). We also report these terms for the months during which oxygen decreases (March to August) in Table S2.

75
For the whole water column (Figure 11a, b), biological processes (PP, WOC, and SOC) greatly exceed physical processes (air-sea exchange and advective transport) in affecting oxygen. PP is always greater than the sum of WOC and SOC in the whole column indicating autotrophy in spring and 80 summer. Advection is negative, acting as an oxygen sink and offsetting 21% of PP on average in the northern and southern regions. Of the two biological oxygen consumption terms (WOC and SOC), WOC accounts for half of total respiration. Negative air-sea flux indicates oxygen outgassing into the at-85 mosphere and is due to photosynthetic oxygen production and decreasing oxygen solubility. However, since hypoxia only occurs in the subsurface, the subsurface budget below is more instructive. When considering only subsurface waters (Figure 11c, d), the influence of PP decreases markedly, ac-90 counting for less than 2% of that in the whole water column. Vertical turbulent diffusion acts as the largest oxygen source in the subsurface layer. SOC is the dominant oxygen sink accounting for 80% of the total biological oxygen consumption. As photosynthetic oxygen production increases gradu-95 ally from spring to summer (Figure 12a, b) WOC and SOC also increase as they are closely associated with photosynthetically produced organic matter. Vertical oxygen diffusion tends to covary with PP, implying an oxygen gradient driven by photosynthetic oxygen production in the upper layer. Lat-100 eral advection of oxygen is negative in March only (early in the hypoxic season) mainly in the southern region but becomes positive later. This suggests that early in the hypoxic season, import of low-oxygen water contributes to hypoxia generation but advection switches to an oxygen source later. 105 Overall, oxygen sources and sink terms are similar in the northern and southern regions.

Discussion
We implemented and validated a state-of-the-art physicalbiological model for the ECS. The implementation is based 5 on a model that was previously developed and extensively used for the northern Gulf of Mexico (Fennel et al. 2011, Laurent et al. 2012, Yu at al.2015b), a region that is similar to the ECS in that it receives large inputs of FW and nutrients from a major river and develops extensive, annually recurring hypoxia (see Table 1 in Fennel and Testa (2019). Our model is more comprehensive than previous models for the ECS. A 6-year simulation was performed and compared to available observations. The model faithfully represents patterns and variability in surface and bottom temperature and salin-15 ity, surface chlorophyll and nitrate distributions, bottom oxygen, and correctly simulates the major current patterns in the region (see Section 3.1 and Supplement). We thus deem the model's skill as sufficient for the analysis of biological and physical drivers of hypoxia generation presented here.

20
The model simulates annually recurring hypoxic conditions but with significant interannual and short-term variability and marked differences in phenology of hypoxic conditions from year to year (Figure 4a, b, c). Interannual variability in hypoxic conditions is much larger than variations in 25 FW input, nutrient load, and bottom oxygen concentrations (Figure 4b) because small variations in oxygen can lead to large changes in hypoxic area when bottom oxygen is near the hypoxic threshold. Interannual variability in hypoxic area is partly explained by variations in annual FW input, con-30 sistent with previous studies (Zheng et al., 2016;Zhou et al., 2017). While the correlation between time-integrated hy- Figure 11. Monthly averaged (2008)(2009)(2010)(2011)(2012)(2013) oxygen budgets for the whole water column and subsurface water from March to November in the northern and southern hypoxic regions. Adv represents lateral advection and lateral diffusion which is comparatively small, while v.mixing represents vertical turbulent diffusion, which is only relevant for the subsurface budget. Thin color bars represent individual years whereas the black bars are the 6-year average. poxic area and FW input is insignificant, there is a strong and significant negative correlation between mean bottom oxygen in August and annual FW input ( Figure 5). Annual FW input is also correlated strongly and significantly with the annually integrated spatial extent of the FW plume, which is a 5 useful metric for extent of the region directly influenced by riverine inputs which induce strong density stratification and high productivity.
Surprisingly, DIN load is not correlated with FW input, hypoxic area, and mean bottom oxygen in August ( Figure   10 5). This is in contrast to the northern Gulf of Mexico where DIN load is highly correlated with both FW input and nutrient load and frequently used as a predictor of hypoxic extent (Scavia et al. 2017;Laurent and Fennel 2019). However, the lack of correlation between hypoxia and DIN load in the 15 ECS should not be interpreted as biological processes being unimportant in hypoxia generation, just that variations in DIN load do not explain year-to-year differences. In fact, hypoxic area and biological rates (i.e. mean August PP, OC, and SOC) are strongly and significantly correlated ( Figure   20 5), emphasizing the dominant role of biological oxygen consumption. The fact that riverine variations in DIN load do not seem to have an effect suggests that riverine nutrient inputs are large enough to saturate the region with nutrients, similar to the northern Gulf of Mexico where small reductions in nu-25 trient load have a relatively small effect (Fennel and Laurent 2018).
Variations in riverine FW input only partly explain interannual variations in hypoxia. For example, the years 2010 and 2012 had similar FW inputs and DIN loads but the hy-30 poxic area was 4 times larger in 2010 than 2012 (Figure 5a). Similarly, 2009 and 2013 had the same FW inputs and nutrient loads but 2009 experienced extensive hypoxia while there was almost none in 2013. In order to elucidate these differences, we investigated biological and physical drivers 35 on shorter time scales.
In the ECS, two distinct zones of low oxygen have been observed (Li et al., 2002;Wei et al., 2007;Zhu et al., 2016Zhu et al., , 2011. The model simulates these two zones, referred to as the northern and southern zones, consistent with observa-40 tions ( Figure 4d) and with generally higher PP and SOC in the northern zone ( Figure 6). Because of these differences we treated the two zones separately in our analysis of intraseasonal drivers.
We found daily biological rates (i.e. PP, OC, SOC) to be 45 significantly correlated with bottom oxygen in both zones, but with relatively large variability around the best linear fit (Figure 7). The biological rates and bottom oxygen are also significantly correlated with the extent of the FW plume (Table 1). Again, these results emphasize the dominant role of 50 biological oxygen consumption, and its relation to riverine inputs, in hypoxia generation but leave a significant fraction of the variability unexplained.
Intra-seasonal variability in hypoxic conditions is significantly related to the extent of the FW plume which is partly 55 explained by variations in riverine FW input but strongly modulated by coastal circulation and mixing. Their influence is elucidated by our analysis of the effects of wind direc-tion and strength on hypoxia. Wind direction has a notable effect on the geographic distribution of hypoxia. Southerly, upwelling-favorable winds lead to a more widespread eastward extension of the FW plume with elevated PP and vertical density stratification (Figure 8). Northerly, downwelling-5 favorable winds create a coastally trapped southward jet that moves FW southward and constrains the plume close to the coast. A similar behavior has been described for the northern Gulf of Mexico (Feng et al., 2014).
Wind strength turned out to be one of the dominant fac-10 tors in hypoxia evolution. We identified high-wind events and showed that whenever bottom oxygen is low, the occurrence of a high-wind event will lead to a partial reoxygenation of bottom waters and decrease hypoxic extent (Figure 9). The impact of high-wind events is also visible in the extent of the 15 FW plume, which is drastically reduced during high winds because FW is mixed. The frequency of high-wind events during summer explains the interannual differences in hypoxic area between 2010 and 2012 ( Figure 10) and 2009 and 2013 ( Figure S8). In 2009 and 2010 there were only few 20 high-wind events during summer while 2012 and 2013 experienced a sequence of storms that led to partial reoxygenation of the water column throughout the summer and thus impeded the development hypoxia. We calculated oxygen budgets for the northern and south-25 ern regions considering the whole water-column and the near-bottom layer only. The subsurface budget is particularly useful in providing insights into when and where lateral advection amplifies or mitigates hypoxia and illustrates that SOC is the dominant oxygen sink in the subsurface. The rel-30 ative importance of WOC and SOC had not previously been quantified for this region due to lack of concurrent WOC and SOC observations and lack of models that realistically account for both processes. The budget for the whole water column is less useful because it is dominated by the oxygen 35 sources, sinks and transport in the surface layer, which does not experience hypoxia and thus is not relevant. The importance of SOC suggested by our model is consistent with recent observational studies in the ECS. SOC on the coastal shelves in the Yellow Sea and ECS has been  Zhang et al. (2017). Simulated SOC in the typical lowoxygen zone falls within the range observed by Zhang et al. (2017) with a mean rate of 20.6±19.2 mmol O 2 m −2 d −1 between April and October. Based on observations, Zhang et al. (2017) already suggested that SOC is a major contributor 50 to hypoxia formation in below-pycnocline waters, which is further corroborated by our model results. It is also consistent with the modelling study of Zhou et al. (2017), who did not include SOC in the baseline version of their model but showed in a sensitivity study that inclusion of SOC simu-55 lates hypoxic extent more realistically. Our results are in line with findings from the northern Gulf of Mexico hypoxic zone where WOC is much larger than SOC below the pycnocline, while SOC is dominant in the bottom 5 m where hypoxia occurs most frequently in summer (Quiñones-Rivera et al., 60 2007;Yu et al., 2015b).
The finding that lateral oxygen transport can act as a net source to subsurface water is also new. On seasonal scales, oxygen advection in the subsurface varies from an oxygen sink in spring to a source in summer, especially in the south-65 ern hypoxic region, implying that the TWC becomes an oxygen source when oxygen is depleted in the hypoxic region. This aspect was neglected in previous studies which only emphasized the role of advection as an oxygen sink promoting hypoxia formation (Ning et al., 2011;Qian et al., 2015). The 70 Taiwan Warm Current originates from the subsurface of the Kuroshio northeast to Taiwan Island, and thus represents an intrusion onto the continental shelves from the open ocean (Guo et al., 2006). In addition to oxygen advection, nutrients are transported supporting PP on the ECS shelves (Zhao and 75 Guo, 2011;Grosse et al., 2020). The intrusion of the Taiwan Warm Current and the Kuroshio accompanied by relatively cold and saline water, and nutrient and oxygen transport, is thought to influence hypoxia development (Li et al., 2002;Wang, 2009;Zhou et al., 2017) but no quantification of the 80 relative importance has occurred until now (see companion paper by Grosse et al., 2020, using the same model).

Conclusions
In this study, a new 3D coupled physical-biological model for the ECS was presented and used to explore the spatial and 85 temporal evolution of hypoxia off the CE and to quantify the major processes controlling interannual and intra-seasonal oxygen dynamics. Validation shows that the model reproduces the observed spatial distribution and temporal evolution of physical and biological variables well. A 6-year sim-90 ulation with realistic forcing produced large interannual and intra-seasonal variability in hypoxic extent despite relatively modest variations in FW input and nutrient loads. The interannual variations are partly explained by variations in FW input but not DIN load. Nevertheless, elevated rates of bi-95 ological oxygen consumption are of paramount importance for hypoxia generation in this region, as shown by the high correlation between hypoxic area, bottom oxygen, and biological rates (PP, OC, SOC) on both annual and shorter time scales.

100
Other important explanatory variables of variability in hypoxia are wind direction and strength. Wind direction affects the magnitude of PP and the spatial extent of the FW plume, because southerly, upwelling favorable winds tend to spread the plume over a large area while northerly, downwelling-105 favorable winds push the plume against the coast and induce a coastal current that contains the FW and moves it downcoast. Wind strength is important because high-wind events lead to a partial reoxygenation whenever bottom oxygen is low and can dramatically decrease the extent of the FW plume. The frequency of high-wind events explains some of the interannual differences in hypoxia, where years with 5 similar FW input, nutrient load, and mean rates of oxygen consumption display very different hypoxic extents because high-wind events lead to partial reoxygenation of bottom waters.
A model-derived oxygen budget shows that SOC is larger 10 than WOC in the subsurface of the hypoxic region. Lateral advection of oxygen in the subsurface switches from an oxygen sink in spring to a source in summer especially in the southern region and is likely associated with open-ocean intrusions onto the coastal shelf supplied by the Taiwan