Coupling and Decoupling of High Biomass Phytoplankton Production and Hypoxia in a Highly Dynamic Coastal System: The Changjiang (Yangtze River) Estuary

The global increase in coastal hypoxia over the past decades has resulted from a considerable rise in anthropogenically-derived nutrient loading. The spatial relationship between surface phytoplankton production and subsurface hypoxic zones is often explained by considering the oceanographic conditions associated with basin size, shape, or bathymetry, but this is not the case where nutrient-enriched estuarine waters merge into complex coastal circulation systems. We investigated the physical and biogeochemical processes that create high-biomass phytoplankton blooms and hypoxia off the Changjiang (Yangtze) Estuary in the East China Sea (ECS). Extensive in situ datasets were linked with a coupled Regional Ocean Modelling Systems (ROMS) and Carbon, Silicate and Nitrogen Ecosystem (CoSiNE) model to explain the connect and disconnect of phytoplankton production and hypoxia. Diatoms were the major contributor of carbon export, and even though phytoplankton concentrations generally were three times greater above the hypoxic zones, high-biomass distributions during the summer-fall period did not closely align with that of the hypoxic zones. A major cause for this decoupling was the non-uniform offshore advection and detachment of segments of water underlying the Changjiang river plume, which carried organic-rich subsurface water northeast of the river mouth. The remineralization of this dissolved organic matter during transit created offshore patches of hypoxia spatially and temporally independent of the nearshore high biomass phytoplankton blooms. The absence of high phytoplankton biomass offshore, and the 1–8 weeks time-lag between the inshore diatom production and offshore hypoxia, made it otherwise difficult to mechanistically explain the in situ observations. The findings here highlight the value of developing integrated physical and biogeochemical models to aid in forecasting the dynamics of coastal hypoxia, under both contemporary and future coastal ocean conditions.

The global increase in coastal hypoxia over the past decades has resulted from a considerable rise in anthropogenically-derived nutrient loading. The spatial relationship between surface phytoplankton production and subsurface hypoxic zones often can be explained by considering the oceanographic conditions associated with basin size, shape, or bathymetry, but that is not the case where nutrient-enriched estuarine waters merge into complex coastal circulation systems. We investigated the physical and biogeochemical processes that create high-biomass phytoplankton production and hypoxia off the Changjiang (Yangtze River) Estuary in the East China Sea (ECS). Extensive in situ datasets were linked with a coupled Regional Ocean Modeling Systems (ROMS) and carbon, silicate, and nitrogen ecosystem (CoSiNE) model to explain the temporary decoupling of phytoplankton production and hypoxia. The CoSiNE model contains two functional groupings of phytoplankton-diatoms and "other"-and the model results show that diatoms were the major contributors of carbon export and subsurface hypoxia. Both observations and simulations show that, although surface phytoplankton concentrations generally were much higher above hypoxic zones, highbiomass distributions during the summer-fall period did not closely align with that of the bottom hypoxic zones. Model results show that this decoupling was largely due to non-uniform offshore advection and detachment of subsurface segments of water underlying the Changjiang River plume. The near-bottom water carried organic-rich matter northeast and east of the major hypoxic region. The remineralization of this particle organic matter during transit created offshore patches of hypoxia spatially and temporally separated from the nearshore high-biomass phytoplankton production. The absence of high phytoplankton biomass offshore, and the 1-8 weeks' time lag between

INTRODUCTION
The flux of anthropogenically derived nutrients into estuaries and coastal oceans has been increasing worldwide over the past few decades (Anderson et al., 2002;Zhou et al., 2008;Conley et al., 2009;Liu et al., 2015). A fundamental effect of this increasing nutrient availability has been more frequent, intense, and widely distributed phytoplankton blooms (Smith, 2003;Anderson et al., 2008;Heisler et al., 2008) that often differ in species composition from earlier times (Glibert et al., 2001;Quay et al., 2013;Jiang et al., 2014). When these high-biomass blooms become nutrientlimited, the algae die and sink below the photic zone, where microbial respiration consumes dissolved oxygen (Officer et al., 1984;Cloern, 2001;Rabalais et al., 2002;Carrick et al., 2005). If bottom water replenishment rates are slow, the elevated decay rates below the pycnocline lead to hypoxia or even anoxia in severe cases (Anderson et al., 2002;Conley et al., 2002;Kasai et al., 2007).
There are numerous examples of increased hypoxia in nearshore waters over the past several decades, including the Baltic Sea (Conley et al., 2002;Carstensen et al., 2014;Neumann et al., 2017), the Chesapeake Bay (Newcombe and Horne, 1938;Hagy et al., 2004;Kemp et al., 2005), the Gulf of Mexico (Rabalais et al., 2002;Turner et al., 2008), as well as many other coastal waters (Legovi and Petricioli, 1991;Li et al., 2002;Dai et al., 2006;Conley et al., 2007;Nakayama et al., 2010;Tishchenko et al., 2011;Zhai et al., 2012;Ram et al., 2014;Zhang et al., 2015;Krogh et al., 2018), and inland waters (Zaitsev, 1992;Zhou et al., 2013). In most of these examples, the long-term increasing trend in hypoxia was shown to link with increased anthropogenic nutrient loading and the development of largescale or frequent phytoplankton blooms. It is believed that the vast majority of anthropogenically enhanced inputs of nitrogen (N) and phosphorus (P) to the Changjiang (Yangtze River) watershed result from the tremendous use of fertilizers, leading to increased N:P and N:silicate (Si) ratios in estuarine waters (Glibert et al., 2006). Altered nutrient ratios ultimately can drive shifts in the dominant phytoplankton speciation away from Si-dependent diatoms that commonly support high fisheries productivities to less important primary producers or even toxic species (Quay et al., 2013;Jiang et al., 2015), although only after the Si concentrations decrease to growth rate-limiting levels. While the overall long-term trends in nutrient concentration and composition show enormous regional variations, strategies that control the input of both N and P are strongly suggested and found to be effective in lessening the negative impacts of cultural eutrophication (Conley et al., 2009;Paerl, 2009). Waters beneath the highly productive surface layer become hypoxic when bacterial respiration associated with the decay of the sinking biomass reduces the dissolved oxygen concentration to <2-3 mg/L, and even anoxic when continued oxygen consumption results in dissolved oxygen <0.5 mg/L (Officer et al., 1984;Chan et al., 2008;Pitcher and Probyn, 2011). Most multicellular organisms cannot survive under hypoxic conditions, and in some areas mass mortalities of fish and benthos typically accompany the onset of hypoxia. However, hypoxia does not necessarily develop in regions with very high rates of natural or anthropogenically fueled primary production as long as the replenishment rates of the deep water are sufficiently strong to supply oxygen in excess of the demand for the decay of organic matter (Ishikawa et al., 2004;Fennel and Testa, 2018). These factors vary over time and often interact with estuary dynamics. The interactive complexity of these processes generates low dissolved oxygen conditions having different spatiotemporal features in coastal regions (Rabouille et al., 2008).
In many cases, the spatial relationship linking surface phytoplankton production and restricted subsurface water exchange can be reasonably predicted by considering the nutrient input rates, regional geomorphology, and oceanographic conditions (Justić et al., 1993;Scavia et al., 2003;Turner et al., 2006). More challenging are those estuarine systems that merge into open coastal or large embayment waters with dynamic circulation patterns. Here, a difference in surface and subsurface advection can lead to a temporal and spatial mismatch or decoupling of surface high-biomass regions and hypoxia zones in coastal waters.
The Changjiang Estuary (CJE) is subject to very high nutrient loading and high N:P and N:Si ratios (Tian et al., 1993;Zhang, 1996;Liu et al., 2003;Glibert et al., 2006;Jiang et al., 2014). Dissolved inorganic nitrogen (DIN) loading from the Changjiang has increased from 20 to 120 µM over the five decades between 1960 and 2010 (Li et al., 2007;Siswanto et al., 2008). In addition to this estuarine inflow, the coastal flow outside the CJE also brings oceanic nutrients to the estuarine region through a complex shelf-circulation system (Figure 1). Of particular importance is the subsurface injection of N and P from the strong western boundary current (the Kuroshio) (Liu et al., 1992;Chen, 1996;Zhang et al., 2007). The Kuroshio water mass that reaches the CJE originated from strong upwelling north of Taiwan (Su and Pan, 1987;Su, 1998;Ding et al., 2019) and contains nitrate concentrations that have increased by 25% over the same time frame (Guo et al., 2012). Once on shelf, the Kuroshio water joins the Taiwan Warm Current (TWC) (Su and Pan, 1987;Su, 1998;Zhu et al., 2004;Jan et al., 2006; FIGURE 1 | The circulation pattern (arrows) in the coastal and shelf waters adjacent to the Changjiang Estuary. The expansion of the Changjiang Diluted Water (CDW) off the estuary depends on interplays among the Yellow Sea Coastal Current (YSCC) from the north, the Zhejiang Coastal Current (ZCC) from the south, the inshore and offshore Taiwan Warm Current (TWC) originated from the strong upwelling north of Taiwan and other forcing. The distributions of harmful algal blooms (HABs) during 1979(Liu et al., 2013 are shown by red squares and those for 2010-2016 by pink circles (Bulletin of China Marine Environment; Supplementary Table S2). The composite distribution of the hypoxic zones in this region since the early 1990s is depicted with green shading (Supplementary Table S1). Black dots represent the station locations for the three repeated surveys referred to in this study. Wang et al., 2013;Wang and Oey, 2016;Xuan et al., 2017; Figure 1). In addition to the TWC, the shelf circulation systems in this area also include the CJE outflow (Zhu and Shen, 1997;Zhou et al., 2009Zhou et al., , 2015Wu et al., 2011), the Zhejiang Coastal Current (ZCC) (Yang et al., 2013), and the Yellow Sea Coastal Current (YSCC) (Zhu et al., 1998; Figure 1). Over the past five decades, these combined inflows have generated a fivefold increase (from 25 to 100) in the N:P and a 20-fold increase (from 0.2 to 5) in N:Si ratios off the CJE (Dai et al., 2010), making it one of the most severely impacted eutrophic regions in the coastal seas of China (Ning et al., 2004;Dai et al., 2013;Lu et al., 2014).
Beyond leading to higher phytoplankton standing stock (Zhou et al., 2008), this coastal eutrophication has resulted in a shift in the phytoplankton community composition in the East China Sea (ECS) (Zhou et al., 2001(Zhou et al., , 2008Yu and Liu, 2016). In particular, there has been an increase in harmful algal blooms (HABs), including the toxic Karenia mikimotoi and Alexandrium catenella as well as the ecosystem-disruptive species Prorocentrum donghaiense and Noctiluca scintillans (Figure 1; Liu et al., 2013;Dai et al., 2014;Lu et al., 2014;Yu and Liu, 2016). Along with changes in plankton assemblages has come the annual development of large-scale hypoxia in the region (Li et al., 2002;Chen et al., 2007;Wei et al., 2007;Wang, 2009;Zhou et al., 2010;Ning et al., 2011;Zhu et al., 2011Zhu et al., , 2016Liu et al., 2012;Wang et al., 2012Wang et al., , 2017Zhu et al., 2017;Luo et al., 2018;Supplementary Figure S1 and Supplementary Table S1). However, the quantitative linkages among these large-scale changes in the frequency and species composition of phytoplankton blooms, dissolved oxygen depletion, and the distribution and expansion of hypoxic areas are poorly understood, mainly due to the limitations of shipbased, discrete, and intermittent sampling of this highly dynamic coastal environment (Figure 1).
Like other anthropogenically affected coastal regions (Silva et al., 2016), marine environmental issues off the CJE demand an enhanced mechanistic understanding of the physicalbiogeochemical processes to enable better forecasting of hypoxic conditions (Zhang et al., 2018). In this work, we use a coupled physical-biogeochemical model to explore the relationship between the spatial-temporal variability of highbiomass phytoplankton production and offshore hypoxia and to reveal the causal linkages among physical processes, nutrient cycling, phytoplankton production, and hypoxia at both seasonal and single-event timescales in the ECS. Although our understanding of this complex system is not yet complete, the findings here help to validate this coupled model approach and its predictive capacity. We hope that this approach can be adopted to support the design of ecosystem restoration strategies to mitigate hypoxia off the CJE.

Field and Satellite Data
Three repeated multidisciplinary field surveys with coarsely spaced sampling stations were carried out off the CJE (122.0-125.0 • E, 27.5-33.5 • N) during June 2-11 (Jun cruise), August 19-30 (Aug cruise), and October 3-13 (Oct cruise) in 2006 (Figure 1, dots). Detailed information about these cruises are given in Zhou et al. (2010) and Zhu et al. (2011). A single intensive survey with high-resolution sampling stations into the Changjiang and in the adjacent coastal region was conducted during July 13-August 23, 2006 Li et al., 2011). Additional cruise observational data were included from May (May 10-June 3)  and September (Zou et al., 2008;Wang et al., 2012) cruises in 2006 to facilitate our understanding of the seasonal development of hypoxia.
Vertical profiles of temperature and salinity were measured using a Sea-Bird conductivity, temperature, and depth (CTD). Discrete water samples were collected using a CTD rosette at near surface, 10 m, within the pycnocline, and near bottom. Chlorophyll a (Chl a) was measured by extraction into 90% acetone using the acidification method (Holm-Hansen et al., 1965). Remotely sensed Chl a data were used to resolve the high spatial-temporal variability of phytoplankton biomass. Daily level 3 surface Chl a data from the GlobColour 1 (hereafter, SAT Chl a), i.e., SeaWiFS, MERIS, and MODIS, were extracted for the study region and showed good agreement with in situ observations (Supplementary Figure S2A). A composite of these SAT Chl a datasets using the weighted averaging algorithm among different satellite sensors was used from the GlobColour. Chl a was overestimated by the satellite where in situ Chl a was <10 µg/L, but was underestimated where in situ Chl a was 1 http://globcolour.info >10 µg/L (Supplementary Figure S2A). In spite of that, SAT Chl a showed overall good agreement (r = 0.68) with in situ Chl a data and the root mean squared difference (RMSD) of 0.39 µg/L (Supplementary Figure S2A).
Vertical profiles of dissolved oxygen were measured with a high-accuracy dissolved oxygen sensor (Sea-Bird SBE43), mounted on the CTD rosette. Dissolved oxygen was also measured from the sampled water by the titration method (Bryan et al., 1976) at discrete depths (surface, two or three depths spanning the pycnocline to improve the sampling resolution according to the pycnocline thickness, and near the seabed) to verify sensor accuracy (Supplementary Figure S2B). The sensor-measured dissolved oxygen data were validated with the titrated dissolved oxygen data and showed a high correlation coefficient of 0.97 with a low RMSD of 0.53 mg/L. Then, the highvertical-resolution sensor-measured dissolved oxygen data were calibrated based on a linear regression y = 0.9063x + 0.7265. Here, y represents the calibrated sensor-measured dissolved oxygen (hereafter, DO) and x represents the titrated dissolved oxygen (Supplementary Figure S2B).

Physical-Biogeochemical Model
The physical processes were simulated by a customized Regional Ocean Modeling Systems (ROMS) (Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008). The model covers the domain of 117.5-132.0 • E and 23.5-41.0 • N, with a curvilinear grid and an approximately homogeneous 4-km horizontal resolution and 30 vertical levels. The model was spun up for 3 years starting from a climatological January mean that was derived from multiyear observations between 1958 and 1987 (Chen, 1992) and driven by Comprehensive Ocean-Atmosphere Data Set (COADS) climatological monthly forcing (Da Silva et al., 1994). The last-year model output was used as the initialization for the realistic run of 2006, for which the surface heat and freshwater fluxes were obtained from the 3-hourly European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim (Dee et al., 2011). The wind stress was calculated from wind vector data as described by Large and Pond (1981), and the wind vector data were from the Blended Sea Winds (Zhang et al., 2006;Peng et al., 2013), which provide daily 0.25 • gridded ocean surface vector wind based on multiple satellite observations. The same wind stress was also used for analyzing the wind mixing. Boundary conditions such as temperature, salinity, velocity, and sea surface elevation were extracted from Global Hybrid Coordinate Ocean Model (HYCOM) climatology averaged from 1979 to 2003 and HYCOM realistic run for 2006 (Kelly et al., 2007), respectively. Both the climatological run and the realistic run include the Changjiang runoff (Table 1). See Zhou et al. (2015) for more details on the circulation model setup and validation. The realistic run of 2006 was driven by the realistic forcing repeatedly for two cycles, and the last cycle output was analyzed in this article.
Observational evidence indicates that diatom blooms are the major source of organic matter export leading to hypoxia off the CJE (Wang et al., 2017). The 13-component carbon, silicate, and nitrogen ecosystem (CoSiNE) model was thus selected, which includes: silicate, nitrate, phosphate, ammonium, picophytoplankton, diatoms, microzooplankton, mesozooplankton, small/suspended particles, large/sinking particles, oxygen, total CO 2 , and total alkalinity. The detailed descriptions of the CoSiNE model are given in Chai et al. (2002) and Xiu and Chai (2014). For more details on the CoSiNE model setup and validation for the ECS and the CJE, see Zhou F. et al. (2017). All biological components were coupled online with physical processes, such as advection and diffusion. The simulated Chl a is referred to as modeled Chl a (MOD Chl a) to distinguish it from satellite-determined Chl a (SAT Chl a) and water column measurements (in situ Chl a).

Modeled Chl a and Currents
The spatial patterns of the bimonthly MOD Chl a (Figures 2A-D) were in reasonable agreement with the composite bimonthly SAT Chl a over the May through October study period, as was the magnitude of production in July and August. Both MOD Chl a and SAT Chl a showed a strong concentration gradient between the nearshore and offshore regions. The MOD Chl a showed high production mostly occurring off the CJE and adjacent nearshore regions. The simulations showed that phytoplankton production in the area had significant seasonality and peaked in July-August, which was similar to those of SAT Chl a. There were some mismatches between MOD Chl a and SAT Chl a. The spring phytoplankton production off the Zhejiang coast occurred later in the simulations than that of the satellite data ( Figure 2E). The high-concentration Chl a tongue extending from the river mouth to beyond the 20-m isoline was underestimated in the simulations (Figures 2F,G). Meanwhile, the modeled Chl a levels were overestimated compared to SAT Chl a values during May-October, particularly in the nearshore regions such as the Subei Shoal and the Qiantangjiang Estuary (Hangzhou Bay) (Figures 2F-H).
The model-simulated current field suggested a prevailing monsoon-driven northeastward current at all layers from the surface to near bottom, namely, the TWC as suggested by Beardsley et al. (1985), Su and Pan (1987), and Su (1998), showing two branches in the southern ECS (Figure 3) similar to the current structure presented by Yuan et al. (1987), Ichikawa and Beardsley (2002), and Xuan et al. (2017). The nearshore branch was along the Zhejiang coast and separated from the offshore branch south of 28 • N. The nearshore branch again bifurcated into two parts south of 30.5 • N, with one part moving northward and pushing the Changjiang Diluted Water (CDW) to north of the river mouth until 33.5 • N, and there the two parts turned together eastward and southeastward. The offshore part turned eastward at 30.5 • N and finally joined the Cheju Warm Current, as suggested by Lie et al. (2000). The northeastward and later southeastward of the CDW tongue north of the river mouth was a well-recognized feature in previous in situ observations (Mao et al., 1963;Beardsley et al., 1985). The simulated current fields at different depths showed some similar patterns, but also significant differences off the CJE. For instance, the bottom layer suggested a notable northwestward onshore flow right outside of the river mouth, while the surface layer showed a bifurcation that extended northeastward and southeastward, respectively. The bottom layer also produced a band of convergence close to the 50-m isobath northeast of the river mouth, and at the convergent band place the CDW at the surface turned eastward and then southeastward. In addition, the 30-m and bottom layers seemed to have more trend of eastward movement in the offshore of the CJE. A high concentration of SAT Chl a continued to occur in June at the surface over an in situ observed hypoxic patch of bottom waters along the Zhejiang coast south of the CJE. The SAT Chl a concentration along the Zhejiang coast (32 µg/L) was similar to that in May, as were the high in situ Chl a levels over the submerged river valley (see Figure 1) off the CJE (40 µg/L) ( Figure 4B). The peak in situ Chl a measured during the cruise (10.5 µg/L) occurred on June 10, which was approximately 25% of the peak value estimated by the satellite on June 21. Episodic high Chl a levels were observed both in situ and from satellite primarily off the Zhejiang coast during both May and June with less variation than in the other regions, consistent with previous observations (Zhou et al., 2003;Zhao et al., 2004;Liu et al., 2013).
In contrast to June, high SAT Chl a (i.e., surface waters) during July coincided with low DO concentrations in bottom waters very close to the river mouth ( Figure 4C). The high in situ Chl a was also observed by the cruise in July (not shown). These cruise data, benefitting from the greater sample station density in the CJE, clearly show the co-occurrence of high phytoplankton biomass and low DO, mainly between the 20-and 50-m depth east of the river mouth. Smaller patches of hypoxia were measured not only in this region but also to the east (126 • E, 31 • N and 124 • E, 31.5 • N) in late July and early August. The highest in situ Chl a concentration measured in July was 21.3 µg/L, in close agreement with the SAT Chl a values (20 µg/L). The Changjiang and Qiantangjiang River estuaries had relatively low in situ Chl a near the river mouth in July (5 µg/L) compared to that occurring in May.
The high SAT Chl a concentration pattern (peak at 15 µg/L) in August was similar to that observed in July; in contrast, the low DO zone expanded significantly in both meridional and zonal directions ( Figure 4D). DO concentrations were <3 mg/L across a spatial zone of 77,100 km 2 , or about 10 times of the hypoxic area observed in July, and hypoxia was also more intense, with the minimum DO deceasing from 2.0 mg/L in July to 1.0 mg/L in August. Although the hypoxic region expanded and covered most of the study area, except close to the 100-m isobath, peak in situ Chl a concentrations decreased to 16.1 µg/L and were slightly lower than those in July. SAT Chl a showed a similar decreasing trend during August.
September brought a drastic decrease in SAT Chl a concentrations, and presumed ventilation of subsurface waters led to some relaxation in hypoxic intensity and contraction of the hypoxic zone. The minimum DO increased to 2.0 mg/L (Zou et al., 2008), and hypoxia north of the river mouth disappeared ; Figure 4E). SAT Chl a in the study region showed a reduction by ∼30% relative to August, and the spatial extent of high biomass (Chl a > 10 µg/L) was less than 5% of that measured in August.
Moderately high SAT Chl a remained at the surface in the area north of the CJE during October, slightly lower than that in August but much higher than that in September ( Figure 4F). The high Chl a patch along the Zhejiang coast that disappeared in September appeared again in October at more moderate levels. The peak in situ Chl a during the cruise (27.7 µg/L) was twice that of the SAT Chl a (14.0 µg/L) averaged over the larger area. The spatial extent of the hypoxic zone increased substantially, from 16,400 km 2 in September to 26,500 km 2 in October, and extended from the Zhejiang coast to the northern tip of the submerged river valley ( Figure 4F).
As noted, the spatial distribution of high Chl a biomass did not match that of hypoxia during the summer-fall stratified period, illustrating a complex interaction between biological and physical processes. This variation is illustrated by the presence of minor hypoxic patches in regions with low Chl a (Figures 4C,D). In addition, the highest Chl a at the surface on the 31 • N transect occurred along the outer edge of the upper front immediately outside of the river mouth, while the minimum DO occurred in the shore side of the front (Figure 5). Indeed, hypoxic bottom waters were overlain by low levels of Chl a biomass more often than that of high-biomass waters, according to the comparison of both SAT Chl a and in situ Chl a at the surface with DO measurements at the bottom layer (Figure 6). Given that advection would have contributed to the variable distribution of the high-biomass phytoplankton patches, the annual spatial dynamics of surface phytoplankton biomass was compared to the

Modeled Hypoxic Zone
As a first step, the low-level DO distribution during the year was extracted from the model to identify the most probable regions of oxygen depletion ( Figure 7A). The modeled low DO concentration generally occurred below 20-m depth. To the north of the CJE, the model results suggested that oxygen depletion was distributed on both the western and eastern sides of the 30-m isobath, with the exception of an isolated low-DO patch further offshore at the 50-m isobath. To the south of the river mouth, the simulated oxygen depletion occurred along the coast between the 20-and 50-m isobaths and at the submerged river valley deeper than 50 m ( Figure 7A).
The maximum hypoxic zone (modeled MHZ) during the year, defined as the assemblage of all grid points where hypoxia was detected in the model-see Zhou F. et al. (2017) for a detailed discussion-was extracted to show the largest extent of hypoxic zone. The modeled MHZ existed as a narrow elongated band from the southern Jiangsu coast to the northern Zhejiang coast ( Figure 7B). Three separate patches could be identified from the modeled MHZ, likely associated with the diversion of the CDW. Running north to south, the first hypoxic water patch, Patch A, lay approximately 300 km northeast of the river mouth near the 50-m isobath ( Figure 7B). It had a long-duration core, reflecting the fact that the CDW passed through this region over most of the season. The second most intense patch, Patch B, lay meridionally distributed in the offshore vicinity of the CJE, while the third hypoxic patch, Patch C, occurred to the south along the Zhejiang coast ( Figure 7B).
To better understand how organic matter production in the upper layers is linked to hypoxia, the modeled oxygen depletion under the pycnocline was compared with the time series of SAT Chl a inside and outside the modeled MHZ within the rectangular zone (120.5-128.5 • E, 25.5-34.5 • N) (Supplementary Figure S3A). The annual mean SAT Chl a within this MHZ was almost three times larger than that outside the region, consistent with hypoxia being attributable to the decay of nutrient-driven phytoplankton production in the overlying waters. There was a positive correlation (r = 0.52) between the fluctuations of the modeled Chl a and SAT Chl a within the MHZ. These SAT Chl a data show that phytoplankton biomass increased slowly from January through April and peaked with a bloom during May-June, followed by lower but still substantial chlorophyll levels in July-August. The modeled Chl a concentrations followed the measured SAT Chl a from January through April, did not capture the bloom in May-June, but reached maximum levels in late July that matched the SAT Chl a concentrations (Supplementary Figure S3A). Thereafter, there was generally good quantitative agreement between the SAT Chl a and modeled Chl a, with diminishing concentrations from August to mid-September, followed by a smaller phytoplankton bloom during October, and then decreasing to low winter values by December. Discrete measurements of the in situ Chl a >10 µg/L during the four cruises of this study (June through October) showed the same general pattern (Supplementary Figure S3B). Augmenting these data with the most current findings of the Bulletin of China Marine Environment (for 2006) also shows a vast bloom extent in May-June.

Modeled Physical-Biogeochemical Processes During Phytoplankton Blooms and Hypoxia
The time series of wind-induced mixing, stratification, nutrients, and phytoplankton succession were investigated in the model domain, using averaged values over the entire modeled MHZ, to investigate the linkages of physical processes, nutrient cycling, phytoplankton blooms, and high Chl a with hypoxia (Figure 8).
The study area is dominated mostly by the monsoon wind patterns, with strong northerly wind occurring during the fall and winter, substantially weaker southerly winds during the summer, and more variable wind directions during the late spring and September-October transitional period ( Figure 8A). The averaged wind speed usually was less than 7 m/s from late spring to summer, with the exception of sporadic tropical storm or typhoon events (depicted by the numbered circles in Figure 8), where wind speeds averaged up to ∼14 m/s within the study area.
The stratification intensity was approximately represented by the surface to bottom difference of density, temperature, and salinity. The weaker winds during the summer generally were not able to disrupt the density stratification of the water column, with only some slight reduction, e.g., in mid-August. The significant and rapid loss of stratification occurred only in early September, associated with the strong Tropical Storm Shanshan and subsequent strong northerly wind (Figures 8A,B). Both temperature and salinity difference contributed significantly to the stratification; however, their temporal patterns varied in complex ways, with salinity difference decreasing when the temperature difference increased during mid-June to early August ( Figure 8B). The reduced contribution of salinity to stratification after mid-June was related to the reduction in riverine input ( Table 1). The salinity difference above and below the pycnocline peaked in mid-June instead of August, in associated with the maximum Changjiang runoff in 2006 (Table 1), which was unusual compared to the climatological mean discharge. The increased runoff during late May to June coincided with an increase in the mean nitrate concentrations in surface waters (Figure 8C), but not in phosphate or silicate due to the phytoplankton uptake of these nutrients during the spring blooms. That is, a greater replenishment of N than P or Si in runoff waters led to more rapid depletion of the latter two as a consequence of phytoplankton growth.
The enhanced vertical mixing that accompanied the strong wind events in July, August, and September increased the modeled nutrient flux into surface waters (Figures 8C,D). The simulated diatom growth increased sharply immediately following these wind-induced mixing events, but reached maximum biomass ∼1-2 weeks after the wind events Figure 8D). These two largest diatom bloom events were followed by two large hypoxia events ∼2 weeks later, consistent with the model findings ( Figure 8E).
The model also shows the succession of diatoms and nondiatoms, where increases in non-diatom biomass followed decreases in diatom biomass ( Figure 8D). The onset of nondiatom blooms occurred rapidly on termination of the diatom blooms during mid-April and mid-May, mid-August, mid-September, late October, and in December. The model-derived time lag between the diatom and non-diatom blooms is within a reasonable range of the field observations (data not shown). Diatom concentrations decreased to a very low level after October, and no blooms occurred during mid-November to late January.
A Hovmöller diagram (which displays time on one axis and distance on the other) along a more or less meridional transect was extracted from the simulation to examine the relationship between high-biomass phytoplankton production and oxygen depletion below the pycnocline (Figure 9). The transect cut through the middle of the long side of the modeled MHZ, with its origin set at the southern end. For simplicity, the vertically integrated diatom concentrations are shown, which illustrate a series of diatom blooms between March and November. Intensive blooms occurred on the along-shore transect at the end of June ( Figure 9A). By mid-July, the vertically integrated phytoplankton production reached 2 mol C/m 2 , assuming a C:N ratio of 9.88 (Harrison et al., 1977), and while oxygen began to be depleted in bottom waters at this time, the area and intensity of oxygen depletion did not enlarge substantially until the end of August/early September. The expansion of hypoxia was enhanced by the second and third high diatom biomass events in late August. The time lag was 1-18 days, but varied at different distances along the transect (Figure 9). The lag phase generally decreased south of the river mouth and increased north of the river mouth, except for the region where no significant diatom blooms occurred. Further correlation analysis showed that there was reasonable spatial agreement between diatom abundances and DO at the bottom layer in the model. The inverse correlation between vertically integrated diatom abundance and DO in bottom water is −0.60 for the whole transect (Figure 9D), while the coefficient is −0.58 between diatoms at the surface layer and DO at the near-bottom layer (Supplementary Figure S4). This inverse correlation was always greater than −0.5 along the transect where both hypoxia and diatoms showed significant concentrations, i.e., the distance from 50 to 150 km, and the strongest correlation was −0.73 at 15 km from the southern end of the transect ( Figure 9D). The phase lag analysis showed that the oxygen depletion mostly occurred a week to 2 months later after the diatom blooms in this transect, where the correlation is >0.5 ( Figure 9E). However, this was not always the case. The relationship between oxygen depletion and diatom abundance was significantly correlated at the distance around 50 km along the meridional transect shown in Figure 9C, but showed no remarkable time lag.
The model also shows the succession of diatoms and nondiatoms, where increases in non-diatom biomass followed decreases in diatom biomass ( Figure 8D). The onset of nondiatom blooms occurred rapidly during mid-April and mid-May, respectively, at the termination of the diatom blooms. The non-diatom blooms occurred throughout most of the record, e.g., February, in mid-May, early August, and in late October. Decreases in the integrated concentrations of both diatoms and non-diatoms over the stratified period occurred after the strong northerly wind event, e.g., from August to October. Diatom blooms peaked approximately 1-2 weeks after the wind-mixing events, while non-diatom blooms generally appeared 1-2 weeks after the diatom blooms. The phase lag between these two blooms is within the reasonable range of field observations (Lu et al., 2000). A similar succession of microalgae blooms was also revealed in other modeling studies of this region (Zhou Z. et al., 2017). The simulated diatom concentrations decreased to very low levels after October, and no blooms occurred during mid-November to late January. Note that the two phytoplankton functional groups differed in their sinking rates, with the diatom detritus sinking velocity being ∼25 m each day (i.e., reaching the seabed within ∼1 day compared to the lower rate of 15 m each day for the non-diatoms).

Coupling and Decoupling Between Surface High-Biomass Diatom Blooms and Bottom Hypoxia
Forecasting the timing, position, spatial extent, and duration of hypoxic zones in the coastal and shelf regions surrounding the CJE is compounded by the dynamic linkages among the physical, chemical, and biological processes. The model results suggest that high-biomass phytoplankton production in this region is regulated by seasonal and intra-seasonal fluctuations in both the anthropogenic nutrient-rich Changjiang (Yangtze River) diluted water and the natural nutrient influx carried by several shelf , and hypoxia at the bottom (E). Stratification, nutrients, phytoplankton concentrations, and minimum DO concentration represent the simulated averages within the study area. The simulated phytoplankton biomass concentrations of diatoms and non-diatoms are in nitrogen units (Chai et al., 2002). currents. The model results also revealed that the event-scale variability in the resulting subsurface hypoxia is characterized by episodic disruption of the pycnocline associated with weather events, which enhances nutrient flux to surface waters and can trigger more severe hypoxia in the region. There exists a straightforward coupling between high-biomass diatom blooms and low DO in underlying waters near the core hypoxic region immediately adjacent to the mouth of the estuary; on average, the modeled phytoplankton Chl a biomass was approximately three times higher over the core hypoxic zone than in the surrounding regions (Supplementary Figure S3). However, there was also significant decoupling of production at the surface (or in the water column) and hypoxic zones during 2006, as seen from in situ observations (Figures 5, 6) and simulations (Figures 9A,  10). In general, subsurface hypoxia was not always correlated with Chl a biomass in the overlying waters across the broader region in the shelf sea outside of the CJE. Portions of the high-biomass blooms, or their underlying waters, often became detached from CJE and were advected along the coast and to the outer shelf (Figures 3, 10). Such episodes were poorly captured by ship monitoring surveys, resulting in an underestimation of the spatial extent of hypoxia (Supplementary Figure S3B). Adequate mechanistic-based forecast of the development of these hypoxic events requires quantitative integration of the physical and biogeochemical processes.
Our model findings showed that diatoms accounted for a large portion of the high Chl a biomass, consistent with in situ observations (Wang et al., 2017), and were the major organic matter export that led to DO consumption below the pycnocline. Based on our simulations, the hypoxia tends to occur 1-8 weeks later than the diatom blooms. The findings here demonstrate how combining remote sensing, numerical simulations, and ship survey data can provide a more comprehensive assessment of the development and evolution of high diatom biomass and hypoxic conditions in this region.

Management Implications
Traditionally, the monitoring and research of "red tides" in Chinese coastal waters has been administratively restricted to toxic algae blooms. However, high-biomass productions that lead to hypoxia also are recognized as an important aspect of HABs having substantial ecological impacts to pelagic and benthic systems. Hypoxia off the CJE has become a seasonal phenomenon, and one where the intensity and spatial scale of FIGURE 10 | Schematic diagram of the coupling and decoupling between high-biomass phytoplankton production (pink squares) and subsurface hypoxic patches (green shading). The decoupling is mostly associated with different current structures at the surface (blue arrows) and subsurface off the estuary (red arrows). Note that the current magnitude might not be represented proportionally by the thickness of the arrow.
hypoxia have been expanding. The definition of HABs in Chinese coastal waters should be broadened beyond only toxin-producing organisms to take into account the hypoxia generated by these massive algal blooms in the CJE and other coastal regions. This change would require expansion of the traditional HAB monitoring from not only nearshore regions but also further offshore where hypoxic events also develop.

CONCLUSION
Abundant supplies of nutrients facilitate a continued presence of high-biomass phytoplankton production and their eventual settling in the region and finally lead to multiple patches of hypoxia. Model results support that there is a time lag of a few weeks between bottom hypoxia and surface diatom blooms. Our findings show that, in most cases, hypoxia was spatially linked to primary production in the overlying waters; however, organic detritus in subsurface waters also were transported offshore with different rates and directions compared to the surface waters in this physically dynamic region, producing offshore patches of hypoxic waters. Advection significantly expands the hypoxic zone relative to the highbiomass phytoplankton zone, generating spatially decoupled relationships between bottom hypoxia and phytoplankton production in surface waters.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
FZ, FC, and MW designed the original ideas presented in this manuscript and coordinated integration of the modeling results and observational data. DH conceived the hypoxia projects and the intensive cruise. FZ, JX, XN, CL, and HL collected the field measurements. FZ, XM, CL, and PW analyzed the field data. FZ, MW, HX, and JS wrote the original manuscript draft. FZ, FC, and HX prepared the physical-biological coupled modeling. FZ, QM, and QZ analyzed the results. All authors contributed to refining this manuscript. ACKNOWLEDGMENTS Field data and model outputs used for the figures and tables are available via contacting the first author. We thank Dr. Zuoyi Zhu for providing DO titration data for calibrating sensors used in the three regular surveys and thank Dr. Jianfang Chen and his team for calibrating DO sensors in the intensive surveys. We also thank NOAA for providing gridded wind vector product via ftp://eclipse.ncdc.noaa.gov/pub/seawinds/.