Impact of Submesoscale Vertical Advection on Primary Productivity in the Southern East China Sea

This study aims to investigate the impact of submesoscale vertical advection (SVA) on the primary productivity in the southern East China Sea. The analysis is based on a comparison between two numerical simulations by using a three‐dimensional coupled physical‐biogeochemical model. One simulation directly resolves SVA on a high‐resolution mesh, and the other leaves SVA unresolved on a low‐resolution mesh. The high‐resolution simulation outperforms the low‐resolution simulation in reproducing the observed chlorophyll distribution, particularly in summer. Resolving SVA results in an approximately 40% increase in primary productivity during the summer, though SVA activity is relatively weak in this season than in other seasons. Among multiscale physical processes, SVA, rather than mixing, is found to be the most important vertical nutrient supply pathway from the nutrient‐rich bottom water to the nutrient‐depleted surface water in summer, particularly on the middle and outer shelves. The impact of SVA on the shelf is unique compared to the open ocean in that it efficiently enhances vertical supply of nutrient‐rich subsurface waters to the nutrient‐depleted surface layer. This study highlights the importance of SVA in promoting primary productivity in stratified shelf seas.


Introduction
Shelf seas cover only 10% of the world ocean area, but they contribute approximately 20% of the global marine primary production . The East China Sea (ECS) is located on a wide continental shelf west of the Pacific Ocean ( Figure 1a). When compared with that of other regions, the primary productivity (PP) of the ECS is at the 75th percentile of the 67 continental margin ecosystems compiled by Sherman (2019). The high PP of the ECS supports several well-known fishing grounds in the region (Gong et al., 2003;Sherman, 2019).
The ECS can be bathymetrically divided into three regions, namely, inner shelf, middle shelf, and outer shelf areas, where bottom bathymetry is between 0 and 50, 50 and 100, and 100 and 200 m, respectively. The circulation system in the southern ECS (see Figure 1a) is characterized by several currents separated by sharp fronts (Belkin et al., 2009). The Zhe-Min Coastal Current flows southward and confines the major influence of the eutrophic Changjiang Diluted Water within the inner shelf (Zeng et al., 2012). The shelf sea circulation system is also subject to the strong impact of the Kuroshio, a dynamic western boundary current; the Kuroshio Current (KC) axis is along the 200-m isobath. Dynamic cross-shelf water exchange occurs along the shelf slope, particularly to the northeast of Taiwan (Su, 1998;Su & Pan, 1987;Zhou et al., 2015). The intruded Kuroshio water mixes with Taiwan Strait water to form the Taiwan Warm Current. The cross-isobath water exchange exhibits a three-layer structure; a year-round onshore bottom intrusion takes place with a mean annual onshore flux of 0.83, 0.71, and 0.48 Sv across the 200-, 100-, and 50-m isobaths, respectively . Kuroshio subsurface water can even be found off the Changjiang Estuary as well as in the Zhejiang coastal area (Wang et al., 2018;Xu et al., 2018;Zhou et al., 2015). The upwelled Kuroshio subsurface water is rich in nutrients (Chen & Wang, 1999;Zhang et al., 2007).
Apart from its large-scale hydrographic features, the southern ECS's ephemeral submesoscale physical processes have attracted increasing attention as observational techniques and numerical modeling capabilities continue to be developed (Wu, 2015;Xuan et al., 2019;Yin & Huang, 2016. The horizontal scale of the submesoscale processes in the southern ECS is approximately 1-10 km (Fox-Kemper et al., 2008;, which is 1 order of magnitude smaller than the Rossby radius of deformation (Lévy et al., 2018;McWilliams, 2016;. Submesoscale processes are visible all over the continental shelf and are reflected in a variety of chlorophyll-a (Chl-a) features appearing as filaments or eddies on remote sensing images (e.g., between 50-and 100-m isobaths in Figure 1b). Studies have suggested that submesoscale processes could play an important role in regulating cross-front water exchange (Wu, 2015;Yin & Huang, 2016. In many regions including the southern ECS, submesoscale vertical advection (SVA) has been proposed to be a dominant contributor to vertical exchange by redistributing heat, salinity, and momentum via reversible vertical water motion (Ferrari & Rudnick, 2000;Haine & Marshall, 1998;Spall, 1997;Thomas, 2005;Xuan et al., 2019).
Based on four seasonal cruises covering the southern ECS, Gong et al. (2003) found the highest PP over the shelf occurs in summer (515 mgC/m 2 /day), and it is the lowest in winter (297 mgC/m 2 /day). Based on a combination of survey and remote sensing data, Hao (2010) estimated the annual mean PP in the ECS to be 372 mgC/m 2 /day, with means of 432 mgC/m 2 /day in spring, 540 mgC/m 2 /day in summer, 342 mgC/m 2 / day in autumn, and 163 mgC/m 2 /day in winter. Based on survey data collected from 18 cruises, Liu et al. (2019) concluded that the PP of the ECS in summer is approximately 1,150 ± 1,190 mgC/m 2 /day, which is 3 times the mean rate in winter (410 ± 320 mgC/m 2 /day). Despite the large variation among the estimates, many studies acknowledged the high PP of the ECS in summer. Studies have suggested that on the inner shelf the year-round highly turbid water influences the light penetration (Shi & Wang, 2010) and the PP is limited by light (Hao et al., 2019); the PP on the middle and outer shelves of the southern ECS is limited by light in winter and by nitrogen in summer (Chen et al., 2001;Hung et al., 2000); off the Changjiang Estuary, PP is mostly limited by phosphate (Wang et al., 2015). Therefore, the nitrogen supplies to the euphotic layer determine the summer PP on the middle and outer shelves, and the sources include the Kuroshio subsurface water (Zhang et al., 2007), Changjiang Diluted Water (Zhang et al., 2007), biological nitrogen fixation (Zhang et al., 2012), and atmospheric deposition (Chen & Chen, 2008). Zhang et al. (2012) found that nitrogen fixation supplied 0.01% to 4.6% of nitrogen demand by PP, which is the same order of magnitude as atmospheric deposition (Nakamura et al., 2005). The major sources for the entire ECS are Kuroshio subsurface water and Changjiang Diluted Water; the extent of nutrients supplied by Kuroshio subsurface water is much larger than that from Changjiang Diluted Water, regardless of the nutrient element ratio and spatial variation (Chen & Wang, 1999;Zhang et al., 2007). Taking into account the strong stratification that suppresses vertical diffusion in summer, Hung and Gong (2011) proposed that typhoon passage can bring nutrient-rich water up, thus potentially acting as a significant nutrient supply pathway, particularly in summer. However, the role of submesoscale processes in the high productivity in summer is not well understood.
We hypothesize that SVA constitutes a significant nutrient supply pathway for the euphotic layer on the broad continental shelf, particularly in a stratified summer environment. Xuan et al. (2019) have shown that the enhanced vertical velocity induced by submesoscale processes, which can penetrate the upper mixed layer, is up to an order of magnitude larger than the mean vertical velocity in the southern ECS. This suggests that SVA might be more efficient in terms of vertical nutrient supply to nutrient-depleted surface water than mixing and mean vertical flow, particularly in summer when the mixed layer depth is shallower than the euphotic layer depth (ELD) (Guo et al., 2019;Lévy et al., , 2018Mahadevan, 2016;Zhang, Qiu, et al., 2019). This study investigates the role SVA plays in the nutrient and phytoplankton dynamics in the southern ECS with a focus on summer.

Summer 2013 Cruise
To understand the large-scale hydrographic conditions in the southern ECS, a multidisciplinary survey was conducted. The in situ data and water samples were collected by a CTD-rosette assembly with Niskin bottles (Instrument: Sea-Bird 911plus); temperature, salinity, photosynthetically active radiation (PAR), and fluorescence data were collected (Sea-Bird Electronics). Seawater samples (0.2 to 1 L) for Chl-a analysis were filtered through GF/F filters. The filters were extracted with 90% acetone at −20°C in the dark for 24 h. Chl-a concentrations were determined by using a Turner Design (Model II) fluorometer that was calibrated with pure Chl-a (Sigma). The Chl-a values converted from the CTD fluorescence measurements during this cruise were in good agreement with the Chl-a analyses of water samples. Thirty-five stations were implemented (see Figure 1b) from 16 to 27 August 2013.
The analysis region in this study for the southern ECS is confined within the latitudinal range of 25-30°N, the longitudinal range of 120-128°E, and the depth range of 0-200 m. Spatially averaged analyses (e.g. time series) are conducted for that region.

Three-Dimensional Physical-Biogeochemical Model 2.2.1. Physical Component
The physical component of the coupled three-dimensional (3D), physical-biogeochemical model used in this study is the Finite Volume Coastal Ocean Model (Chen et al., 2003). The model domain covers the Bohai Sea, Yellow Sea, ECS, and part of the Pacific Ocean, the longitudinal and latitudinal ranges of which are approximately 118°E to 136°E and 21°N to 41°N, respectively (Xuan et al., 2019). We established a large area model to have a lateral boundary that is far from the study area to reduce the error from open boundaries. The general bathymetric chart of the oceans (Smith & Sandwell, 1997) provided the high-resolution (0.5 km) bathymetric data. A sigma-stretched coordinate system was used to specify 40 vertical layers with a vertical resolution of approximately 1-5 m in the mixed layer on the shelf. The vertical diffusion coefficient was calculated according to the K-profile parameterization scheme proposed by Large et al. (1994). The model is driven by the daily mean heat fluxes from objectively analyzed air-sea fluxes (Yu & Weller, 2007) along with 3-hourly wind stress data from the ERA-Interim reanalysis data set (Dee et al., 2011). The driving forces also include tidal forces from TPXO.7.0 (Egbert et al., 1994;Egbert & Erofeeva, 2002) and Changjiang River discharge from the Datong Hydrological Station (http://yu-zhu.vicp.net). The temperature, salinity, and currents at open boundaries were derived from the Hybrid Coordinate Ocean Model (Bleck, 2002;Chassignet et al., 2007).
The model was set up with two horizontal resolutions, 0.5 and 10 km, which are referred to as the high-resolution model (HRM) and low-resolution model (LRM), respectively. The gridded depth data were interpolated from the bathymetric chart following the volume conservation principle, thus enabling the coincidence of grid cell volumes from both the HRM and the LRM. The resolution-dependent horizontal diffusion coefficient was calculated according to the canonical parameterization scheme proposed by Smagorinsky (1963). The submesoscale activities are resolved in the HRM but not in the LRM. The relative influence of submesoscale processes on the hydrographic environment preserved within the HRM can be distinguished by comparing the HRM and LRM output.
Hourly outputs of sea surface height, temperature, salinity, and velocity for 2013 and 2014 were produced following 3 years of spin-up from initial conditions. Since our analysis focuses on daily data, a tidal filter (Godin, 1972) was used to remove tidal signals. Simulation output of both the HRM and the LRM was used to force the 3D biogeochemical model.

Biogeochemical Component
We used a modified version of the model of Hu et al. (2004), which is a nutrient-phytoplankton-zooplankton-detritus-type model (Riley, 1949). We removed the silicon cycle in their biogeochemical model, because the main limiting nutrients are nitrogen and phosphorus in the study area (Chen et al., 2001;Hung et al., 2000;Wang et al., 2015). The module for PAR as well as some parameter settings were adapted for the ECS (see Appendix A). The ecological state variables are the phytoplankton concentration (P), zooplankton concentration (Z), organic detritus concentration (D), dissolved nitrate concentration (N n ), and dissolved phosphate concentration (N p ). Chl-a is the most common light-absorption pigment within phytoplankton, and it is known to be a good proxy for phytoplankton biomass (Boyer et al., 2009;Huot et al., 2007;Van De Poll et al., 2013). In this study, the phytoplankton concentration is converted to Chl-a concentration via a constant N-to-Chl-a ratio (c n ; units of mmol N/(mg Chl-a)) of the phytoplankton nitrogen content to Chl-a. Gowen et al. (1992) used a large data set to determine the range of c n between 0.23 and 4 mmol N/(mg Chl-a); Liu et al. (2010) applied the range of 0.2857 to 1 mmol/mg in a biogeochemical model for the ECS. Based on survey data for the ECS during different seasons, Liu et al. (2019) found the C-to-Chl-a ratios (weight : weight) within the plume and shelf waters to be relatively stable in both summer and winter, resulting in c n ≈ 0.377 mmol N/(mg Chl-a) via the Redfield ratio (Redfield, 1934). In this study, c n is set to 0.3 mmol N/(mg Chl-a), which is consistent with that used by Wan et al. (2000) and Hu et al. (2004). Processes formulated in the model include the following: growth, respiration, grazing, mortality of phytoplankton, excretion, mortality of zooplankton, sinking, and remineralization of organic detritus. The phytoplankton growth is limited by either the availability of nutrients or the availability of PAR by introducing a limiting coefficient.
The governing equations, parameterization schemes, and parameters are elaborated in Appendix A. The initial conditions and open boundary conditions for N n and N p were extracted from the WOA2013v2 data set (Garcia et al., 2014). Zero-gradient boundary conditions were applied for other variables. The monthly data of nitrate and phosphate concentrations at Changjiang river mouth were observed by the Second Institute of Oceanography, Ministry of Natural Resources. The simulation time step of the biogeochemical model is 90 s, and daily mean outputs of biogeochemical state variables and fluxes are available.

Validation
The circulation model has been validated for the ECS in terms of the temperature-salinity field, tides, and currents (Xuan et al., 2016(Xuan et al., , 2017(Xuan et al., , 2019. Based on observations in 2013, the mean temperature deviation in the southern ECS was <1°C, and the standard deviation of the observational data ranged from 0.  Figures 2a and 2b, respectively. The density in summer monotonously increases with depth, indicating stable stratification ( Figure 2a). There are two dense cores, one nearshore and one offshore. The pycnocline deepens with depth between 0-and 150-km offshore and shoals from 200-km offshore. Both the HRM and the LRM reproduce the stratified hydrographic environment. Looking into the simulated Chl-a at sections A1-A7 of the cruise (Figure 2b), neither the HRM nor the LRMs captures the high Chl-a at the surface on the inner shelf. Both (HRM and LRM) reproduce the subsurface chlorophyll maximum (SCM; Cullen, 2015) characteristic on the middle shelf. The HRM results agree better with cruise data than do the LRM results in terms of the higher vertical extent of the SCM layer and its magnitude. The survey stations are too sparse to reveal the submesoscale spatiotemporal variation. Figure 3 shows the scattered profiles of Chl-a and PAR as measured by the summer cruise on the inner, middle, and outer shelves. The SCM can be found but is not a common characteristic on the inner shelf ( Figure 3a). The coastal region (inner shelf) is the most (directly) anthropogenically affected. The interplay among nutrient-rich riverine water, high turbidity, and strong tidal mixing results in the high spatiotemporal variation of phytoplankton (Hao et al., 2019). Our model does not fully capture the high variability on the inner shelf. There is an SCM at almost all stations on the middle shelf and outer shelf (Figures 3b  and 3c). The simulated results fall in the variation range of the Chl-a profiles observed on the middle shelf ( Figure 3b). Figure 3c shows that the HRM and LRM reproduce the SCM characteristic on the outer shelf in terms of the magnitude of SCM. Although the HRM performs better than the LRM, the depth of the SCM is not well simulated by either. Mismatch on the outer shelf probably originates from the constant ratio of Chl-a to nitrogen used in this model (Steele, 1962(Steele, , 1964. Figures 3d-3f indicate that the modified PAR model results in higher PAR attenuation rate in coastal water and slower attenuation rate offshore, which agrees with the field measurements. The cruise data also illustrate that the actual PAR attenuation rate on the inner shelf is higher than the simulated rate unlike the good agreement on the middle and outer shelves. The performance of the HRM and LRM is further assessed in terms of the Chl-a horizontal distribution. The vertically averaged Chl-a within the euphotic layer (ELD mean Chl-a) based on cruise data is shown in Figure 4a, and the 12-day composite GOCI image for Chl-a during the cruise is shown in Figure 4b. The ELD is defined as the depth at which PAR reaches a value equal to 1% of its surface value (Morel, 1988). The high horizontal variation of Chl-a on the inner shelf found by the survey (Figure 4a) cannot be seen in the remote sensing data (Figure 4b) or the model outputs (Figures 4c and 4d). The remote sensing data overestimate Chl-a on the inner shelf, which is in part caused by interference with the suspended  Figure 4b are compared with outputs from the HRM and LRM point-to-point. Since Chl-a on the inner shelf is overestimated by the remote sensing product (Figures 4a and 4b), samples are selected from the middle and outer shelves. As shown in Figure 4e, the coefficient of determination R 2 for the HRM data (95% confidence interval is 0.58-0.61) is much higher than that for the LRM data (95% confidence interval is 0.36-0.40), suggesting that the patchy Chl-a distribution is better captured by the HRM than by the LRM. Note that the LRM underestimates the Chl-a at many positions, while the HRM overestimates the Chl-a at many positions ( Figure 4e). Liu et al. (2019) found that the variation range of Chl-a in the southern ECS in summer is up to 100-fold. The model relatively well reproduces the horizontal distribution of Chl-a over most of the southern ECS in summer despite not fully capturing the highly variable Chl-a on the inner shelf. By integrating, the HRM estimates that there is 7.70 × 10 12 mg Chl-a in the analysis region, which is approximately 40% higher than that determined by using the LRM (5.56 × 10 12 mg). The Chl-a on the middle shelf accounts for 72% of all Chl-a in the analysis region.
The observed growth and grazing rates in the ECS are further compared with those in the model. Zheng et al. (2015) investigated the specific growth and specific grazing rates in the ECS via experiments during cruises. In summer, the mean specific growth rate and specific grazing rates at the surface layer are 0.77 ± 0.53 and 0.69 ± 0.42 day −1 , respectively, based on data from 26 stations; in winter, the mean

10.1029/2019JG005540
Journal of Geophysical Research: Biogeosciences specific growth rate and specific grazing rates at the surface layer are 0.39 ± 0.18 and 0.21 ± 0.08 day −1 , respectively, based on data at 24 stations. The growth and grazing rates in the model are not constant but vary over space and time, and the cruise data set is not large enough for comprehensive comparison. We average the values from surface to 20-m depth during summer and winter, respectively, to ascertain whether they are on the same order of magnitude as the cruise data. The mean growth and grazing rates in the HRM for summer are 0.79 ± 0.46 and 0.42 ± 0.37 day −1 , respectively; the mean growth and grazing rates in the HRM for winter are 0.59 ± 0.24 and 0.29 ± 0.23 day −1 , respectively. Given the large standard deviation of both observed and simulated rates, the growth and grazing rates determined by the model can be assumed consistent with the cruise data. As compared with the standard product, the data set calibrated by the ship-measured Chl-a reduced the bias from 35% to 8% and the root-mean-square log error from 46% to 32%. The daily product from 1998 to 2013 was averaged to obtain the climatological monthly mean Chl-a areal distribution. The climatological remote sensing data reveal the general seasonal variation pattern of surface Chl-a. The simulated ELD mean Chl-a concentrations averaged over the analysis region over four seasons are compared with those from the remote sensing data in Figure 5. In spring, summer, and autumn, both the HRM and LRM averages fall within one standard deviation of the remote sensing data, and the HRM mean is closer to the remote sensing mean. In contrast, ELD mean Chl-a in winter simulated by either the HRM or the LRM is lower than the remote sensing result by more than one standard deviation. In winter, the Chl-a is more concentrated near the surface due to low PAR, and the monthly averaged Chl-a at the surface in the analysis region can reach as high as 0.4 mg/m 3 , which is very close to the value extracted from the remote sensing data. Nevertheless, caution is necessary when interpreting the simulated winter results on account of this mismatch.

Weak SVA Activity in Summer
Seasonal variation in the SVA intensity in terms of magnitude and occurrence frequency in the southern ECS has been comprehensively studied by Xuan et al. (2019). They also diagnostically evaluated the contributions to SVA by examining four mechanisms: frontogenesis, mixed layer instabilities, turbulent thermal wind, and topographic wakes. Although these mechanisms are complex, it was found that a deeper mixed layer depth generally correlated with more active SVA. Four HRM snapshots of the vertical velocity at 20-m depth, each corresponding to one of the four seasons, are chosen to illustrate the SVA activity in the upper layer. To distinguish SVA from the mean flow field, a low-pass space filter on the vertical velocity field is applied with a cut-off scale of 50 km (Sasaki et al., 2014), and the short-wavelength part is taken to represent SVA, as shown in Figure 6. The long-wavelength part, designated as mean vertical flow, is shown in Figure 7. Overall, the maximum magnitude of SVA is 1 order of magnitude larger than the mean vertical flow in the southern ECS. Notably, SVA occurs least frequently in summer (Figure 6c).   PP in this paper is defined as the carbon produced by phytoplankton minus its respiration per day. It is conventionally referred to as net PP, which indicates the rate of available energy fixed by marine autotrophs for higher trophic levels. The simulated nitrogen within phytoplankton is converted to carbon biomass for the calculation of PP via the Redfield ratio (Redfield, 1934). Monthly and spatially averaged PP in the analysis region, simulated by the HRM and LRM, is plotted in Figure 8. To determine whether the phenomenon is unique in 2013, we also run the models for 2014. The seasonal variation of PP simulated by the HRM is consistent with previous studies, exhibiting a high value in summer. By subtracting the PP simulated by the LRM from the HRM result, it is clear that the most significant difference occurs in summer. The annual mean PP simulated by the HRM is higher than that simulated by the LRM by approximately 17%. In particular, PP in summer simulated by the HRM is approximately 40% higher than that simulated by the LRM. Since the HRM is better at resolving submesoscale processes than is the LRM, the impact of SVA on PP can be studied by comparing the HRM and LRM results. It is of interest to understand the potential relationship between weak SVA activity and high PP in summer.

SVA as an Important Vertical Nutrient Supply Pathway in Summer
The HRM yields a quite different vertical flow field by resolving SVA when compared with that of the LRM result (Figures 9a and 9b). SVA establishes many upwelling and downwelling vertical water exchange pathways with the horizontal extent of O(10) km at section A on a summer day. The pathways can reach both the surface and bottom mixed layers. The maximum value of the total vertical advection in the HRM is twice that in the LRM and an order of magnitude larger than the mean vertical flow. Figures 9c and 9d show the vertical nitrate transport on the same day simulated by the HRM and LRM, respectively, at section A by multiplying the vertical velocity and nitrate concentration. The influence of SVA on nutrient dynamics depends on the vertical gradient of the nutrient profile. The HRM can resolve many submesoscale vertical nitrate transport pathways that penetrate from the surface mixed layer to the bottom mixed layer. In contrast, the vertical advection of nitrate in the LRM is very weak. The maximum daily-averaged vertical nitrate transport in the HRM is approximately 8 times that in the LRM. Figures 9e and 9f show PP at section A on that day in the HRM and LRM, respectively. We shall look into potentially different response of PP to SVA from  (Figures 9a and 9b); the nutrient transport is also little influenced by SVA (Figures 9c and 9d); therefore, the PP is almost unchanged in the presence of SVA (Figures 9e and 9f). From 50 to 200 km, SVA gradually becomes more and more active (Figures 9a and 9b), which is also reflected in the fluctuation of the isopycnals (Figure 9a); the vertical nutrient transport becomes noticeable in the HRM (Figures 9c and 9d); however, the PP simulated by the HRM is not quite different from that in the LRM (Figures 9e and 9f); that is probably due to the high turbidity in the coastal water. From 200 to 400 km, various submesoscale vertical circulations result in either small fluctuations or sharp oscillations of the isopycnals (Figures 9a and 9b); nutrient transport pathways are established on the basis of both the SVA and the vertical nutrient gradient, and many of them can reach as deep as the bottom boundary layer (Figures 9c and 9d); big difference in the PP simulated by the HRM and LRM (Figures 9e and 9f) appears to be associated with SVA reflected in the oscillatory pattern in Figure 9e; the subsurface high PP band widens downward in the HRM as compared with that in the LRM (Figures 9e and 9f), and it is probably stretched by SVA (Figures 9a and 9b). The overall PP in the HRM is higher than that in the LRM, suggesting a positive contribution by SVA. The impact of SVA on nutrient dynamics relating to PP in summer is further analyzed quantitatively in the following.
Nutrient dynamics are governed by not only biogeochemical processes but also multiscale physical processes including mean circulation, mesoscale eddies, submesoscale processes, and microscale mixing. The evolution of nutrient content governed by those processes can be represented by the Reynolds-averaged equation by decomposing a state variable into a mean part that is averaged over a spatiotemporal scale that is much larger than the submesoscale field and an eddy fluctuation part.  omitted negligible terms including the large-scale upwelling/downwelling term, mesoscale vertical advective term, submesoscale

10.1029/2019JG005540
Journal of Geophysical Research: Biogeosciences horizontal advective term, and horizontal mixing term to quantify the relative contributions of multiscale physical processes to the nutrient dynamics. Equation 1 includes the large-scale upwelling/downwelling term, comparing to that put forward by : (1) where t is time, z is vertical coordinate, u is the horizontal velocity containing the zonal and meridional components, w is the vertical velocity, k z is the vertical diffusion coefficient, ∇ H denotes the operation to obtain the horizontal gradient, the overbar denotes the spatiotemporal averaging operation, and the prime denotes the eddy fluctuation. The long-term change of local nitrate concentration (Local term) is approximately balanced by transport by the large-scale circulation system (Mean term), mesoscale horizontal advections (Meso term), SVAs (SVA term), vertical diffusion due to microscale turbulent mixing (Diffusion term), and biogeochemical processes (Bio term).
To quantify the relative contributions of several physical nutrient pathways to the surface water, the Reynolds-averaged equation for nitrate is integrated over the depth range of 0-20 m from 1 to 31 August 2013, as in Equation 2: (2) The integrated terms based on the HRM and LRM results are plotted in Figures 10a-10l. Figures 10a and  10b suggest that the surface water of the analysis region slightly accumulates nitrate in August. The HRM and LRM exhibit similar patterns for transport of nitrate due to the large-scale circulation system (Figures 10c and 10d). Figures 10c and 10d suggest that the eutrophic Zhe-Min Coastal Current injects nitrate into the inner shelf, whereas the oligotrophic Taiwan Warm Current takes nitrate away from the middle and outer shelves. Figures 10e and 10f suggest that the mesoscale horizontal advections make minor contributions for nitrate transport in both the HRM and the LRM. Figure 10g suggests that the SVA contribution to the nitrate supply is predominant in the HRM particularly on the middle and outer shelves. On the contrary, Figure 10h shows the SVA contribution to the nitrate budget is minor in the LRM. The SVA contribution in the HRM is an order of magnitude larger than that in the LRM (compare Figures 10g and 10h). As another vertical nutrient supply pathway, the diffusion in the HRM (Figure 10i) provides slightly less nitrate to the surface water than does that in the LRM (Figure 10j). Therefore, in the HRM, SVA is the major vertical nitrate supply pathway, which contrasts the LRM that determines the major contribution to be diffusion. The amount of nitrate removed from the surface water by the biogeochemical processes in the HRM (Figure 10k) is much greater than that in the LRM (Figure 10l). This suggests that SVA is a significant factor for the high PP in the analysis region in summer as it relieves nutrient limitation; this mechanism is predominant on the middle and outer shelves.

Implication of Model Results and Limitations
The main purpose of our study is to underline the importance of SVA for ecosystem dynamics. Although both the HRM and the LRM underestimate the stratification (Figure 2a), the HRM performs much better than the LRM for simulating Chl-a (Figures 2b and 3). One major cause for the differences in simulating biogeochemical processes between the LRM and the HRM is the difference in nutrient supply to the nutrient-depleted surface (in "upwelling" SVA) and the downward transport of biomass (in "downwelling" SVA). Different availability of nutrients results in up to 40% difference in PP at a regional scale. Therefore, our numerical experiments demonstrate the necessity to include the impact of SVA in modeling ecosystem dynamics of coastal shelf seas.

Journal of Geophysical Research: Biogeosciences
A way to improve the accuracy of the LRM is to optimize subgrid-scale process parameterizations. For example, both Schrum (1997) and Xuan et al. (2019) incorporated submesoscale effects in coarse-resolution modeling to improve the simulated thermohaline structure in shelf seas. This paper reveals further challenges on the parameterization work related to ecosystem dynamics. Xuan et al. (2019) found a larger impact of SVA on the interior thermal structure of the southern ECS in autumn than in other seasons, whereas this study found that even weak SVA in summer can exert great impact on the PP of the southern ECS, which means the impact on each state variable should be carefully examined. The parameterizations of biogeochemical processes during the submesoscale circulation should also be taken into account when incorporating the SVA term into other terms, for example, the Diffusion term.
The comparisons for Chl-a, growth rate, and grazing rate in section 3.1 suggest that the results of our model can be regarded as a first-order estimate of the ecosystem dynamics. Several potential model limitations of the biogeochemical component due to simplifications are worth noting. First, the phytoplankton adaption itself cannot be depicted by the constant nitrogen-to-Chl-a ratio. This probably is the reason for the shallower depth of the SCM in the simulations than the observation (Figure 3c). Second, the model is incapable of simulating complicated and intertwined community structure change from summer to winter Zheng et al., 2015) with only one phytoplankton species and one zooplankton species. This probably results in the underestimation of Chl-a by simulations in winter ( Figure 5). Third, the phytoplankton self-shading effect was not taken into account. From Figure 3d, it can be seen that the PAR is overestimated on the inner shelf, which is likely partially owed to this simplification. Light limitation would be stronger if self-shading was included in the HRM, which would likely offset part of the positive impact of SVA on PP. More work, for example, introducing self-shading, more nutrients, and more species, needs to be done to reveal the SVA impact on the more complex ecosystem dynamics.

SVA Impact on the Phytoplankton Dynamics in Stratified Shelf Seas
As introduced, the main purpose of this study is to investigate whether SVA can relieve the nutrient limitation for PP in the southern ECS in summer. The simulated results support that hypothesis. It is worthy to further discuss how SVA achieves the overall positive impact on PP via coupled physical-biogeochemical processes and the distinct features compared to those in an open ocean context. SVA first influences the hydrographic environment of the southern ECS in summer. In Figure 11a, the density profiles averaged on the middle shelf simulated by the HRM and LRM deviate more in the interior and the bottom boundary layer than in the upper layer. The lower vertical density gradient in the HRM suggests that SVA increases the mixing near the bottom layer. The change in hydrographic environment (density) suggests the further influence on the nutrient transport. The vertical nutrient transport capacity by SVA is also controlled by the background distribution of the nutrient. Figure 11b shows that the biggest deviation, in terms of nitrate or phosphate profiles simulated by the HRM and LRM, occurs in the interior rather than the surface or bottom boundary layers.

Journal of Geophysical Research: Biogeosciences
Whether the SVA exerts a positive or negative net effect on the PP may be site and time specific. SVA is often induced by the submesoscale front via frontogenesis (McWilliams, 2016). The submesoscale fronts due to inhomogeneous mixing in the upper ocean are generally ephemeral and shallow (Fox-Kemper et al., 2011). The SVA structure is often featured by an upwelling branch and a downwelling branch within a cross-frontal ageostrophic secondary vertical circulation (Lévy et al., 2018). Looking into one submesoscale ageostrophic vertical circulation, the upwelling branch pumps nutrient up to the euphotic layer to facilitate phytoplankton growth, but the downwelling branch transports phytoplankton downward and reduces effective light availability of phytoplankton. Therefore, the upwelling branch imposes a positive effect on phytoplankton growth, while the downwelling branch imposes a negative effect. The background vertical mean flow is then important for the sign of the relative impact of SVA. In regions characterized by strong large-scale upwelling, nutrient supply by the upwelling branch of SVA is minor compared to that of the mean flow; subduction of phytoplankton to deeper water with less PAR by the downwelling branch of SVA becomes prominent in reducing the PP (Lathuilière et al., 2010). In contrast, SVA would be a very important vertical nutrient supply pathway to the euphotic layer in regions characterized by large-scale downwelling or weak upwelling (Lathuilière et al., 2011). Figures 7c, 9a, and 9b suggest that the mean vertical flow in the southern ECS is weak in summer. Combining Figures 8, 10g, and 10h, it is concluded that the vertical nutrient supply by SVA results in an overall positive impact on PP. From Figures 9a, 9b, 9e, and 9f, the downward widening

10.1029/2019JG005540
Journal of Geophysical Research: Biogeosciences of the subsurface high PP band should be due to the downward SVA motion. It is also reflected in the Chl-a profiles on the middle shelf as shown in Figure 11a. The higher PP resulting from SVA also supports more zooplankton and produces more detritus than circumstances without SVA (Figure 11c). Taking the circumstance on the middle shelf as an example, we can, to some extent, partially unravel why weak SVA in summer exerts significant impact on the PP of the southern ECS in summer: There exists large vertical nutrient gradient there, and SVA can reach the interior and even the bottom boundary layer to supply nutrient-rich deep water to the nutrient-depleted surface layer; in addition, PAR is sufficient. Therefore, in a strongly stratified shallow shelf sea, even weak SVA could exert tremendous positive impact on the PP. The impact of SVA on the middle shelf of the southern ECS in summer is conceptually sketched in Figure 12.
This study and other recent studies (Dauhajre et al., 2017(Dauhajre et al., , 2019Xuan et al., 2019) suggest that SVA may change the interior and even to the bottom boundary layer structure on the continental shelf. Shelf seas with active submesoscale processes and high PP can be also found elsewhere, for example, in the North Sea (Sherman, 2019;Zhao, 2019;Zhao et al., 2019), South China Sea , and Iberian shelf (Zhang et al., 2016), and reflected in the heterogeneity of organic matter deposition on surface sediments (Zhang & Wirtz, 2017;Zhang, Wirtz, et al., 2019). The significance of SVA for the phytoplankton dynamics is therefore believed to be typical in stratified shelf sea environments.
Many previous studies on the contribution of SVA to PP as well as relevant biogeochemical fluxes were conducted in an open ocean context. They mainly focuses on the upper-ocean dynamics when involving SVA (e.g., Guo et al., 2019;Lévy et al., 2001;Mahadevan & Archer, 2000;Spall & Richards, 2000;Uchida et al., 2020). For example, SVA plays a role in regulating the upper-ocean mixed layer depth. SVA always tends to flatten the density isosurfaces in the front and subsequently reduce the mixed layer depth via an inverse cascade of turbulent energy (Thomas & Ferrari, 2008). When the mixed layer depth is deeper than the ELD, the reduced mixed layer depth increases the light exposure time of phytoplankton and promote early local phytoplankton blooms (Lévy et al., 2005;Mahadevan et al., 2012). In another scenario where the mixed layer depth is shallower than the ELD, the SVA that penetrates the mixed layer becomes important for the vertical nutrient supply other than mixing. However, whether the SVA can effectively transport the nutrient depends on the relative depth of the lower extent of the ageostrophic secondary vertical circulation and the nutricline (Lévy et al., 2018). Although the submesoscale front can reach as deep as 100 m in the open ocean, it does not always extend to the nutricline to establish an effective nutrient supply pathway . Liu and Levine (2016) questioned the importance of SVA on the phytoplankton dynamics in summer. They found negligible impact of SVA on Chl-a concentration within the North Pacific subtropical gyre in summer. Lévy et al. (2018) argued that "submesoscale dynamics might not be as important as previously thought for marine productivity." Alternatively, persistent fronts associated with intensified large-scale currents like KC can reach the deep nutrient pool and establish a persistent nutrient supply pathway . With increasing computing capacity, the resolution of global ocean numerical models certainly will increase and both submesoscale fronts and persistent fronts will be better resolved. Lévy et al. (2018) reminds that the submesoscale fronts might not be an important vertical nutrient transport pathway to drive changes in global marine productivity in higher-resolution models, but the strengthening of persistent fronts should be also paid attention to. The present study suggests that SVA impact on PP is significant in shelf seas when stratification is strong, because SVA can regulate the nutrient distribution within the whole water column due to the shallow depth of shelf seas.
The discussion in this study is mostly confined in the vertical transport since we do find that SVA is a predominant vertical nutrient supply pathway within the analysis region in summer as illustrated in Figures 10g and 10i. However, it only reveals the tip of the iceberg for comprehensive understanding of the "role of submesoscale currents in structuring marine ecosystems" (Lévy et al., 2018). For example, a study has suggested that the mean circulation system could be modified by resolving the submesoscale currents . A clue for modification of the circulation system can be also found in Figures 10c and 10d by decomposing the nitrate transport. The seasonal variation revealed by Figures 5-8 is also of interest. Therefore, much more work should be done to complete the entire picture.

Conclusions
To understand the mechanism for high PP in the southern ECS in summer, the impact of SVA was studied based on a 3D physical-biogeochemical model. Two runs were carried out with different resolutions; the HRM was able to resolve SVA, whereas the LRM could not. Comparison between the simulation results and the summer 2013 cruise and remote sensing data show that the HRM reproduces the Chl-a distribution better than does the LRM, particularly in summer. The SVA activity is relatively weak in summer in terms of magnitude and occurrence frequency, as compared to that in other seasons. However, the HRM yielded an approximately 40% higher PP than that of the LRM in summer. By analyzing the nutrient dynamics in summer based on the Reynolds-averaged equation, we found that SVA, rather than diffusion, acts as the most important vertical nutrient supply pathway for the nutrient-depleted surface water in summer on the middle and outer shelves. The impact of SVA on the shelf is unique compared to the open ocean in that it efficiently enhances vertical supply of nutrient-rich subsurface waters to the nutrient-depleted surface layer. The mechanism for regulating phytoplankton dynamics by SVA in a stratified shelf sea is as follows: (1) the upwelling branch induced by SVA transports nutrient-rich bottom water to the euphotic zone to promote PP and (2) downward motion associated with SVA results in a reduced entrained phytoplankton growth rate due to light limitation. However, the effect of (2) is minor compared to that of (1) in the southern ECS in summer when not taking into phytoplankton self-shading. Therefore, SVA is critical in supporting the high PP of the southern ECS, especially in summer. Our finding suggests an added value by resolving SVA in numerical modeling for a better assessment of PP in coastal shelf seas.

Appendix A: Governing Equations of the 3D Biogeochemical Model
Bottom-up control for phytoplankton growth is expressed by the Michaelis-Menten equation (Dugdale & Goering, 1967). However, the phytoplankton growth rate may be limited due to the absorption inefficiency of any given macronutrient. The nutrient-limiting coefficient is obtained by taking the minimum of the limiting factors in terms of N n and N p , i.e., where p 1 and p 2 are the half-saturation coefficients for nitrate uptake and phosphate uptake, respectively.
The PAR model was customized. Rather than model the complex dynamics of suspended particulate matter that influences turbidity, we used remote sensing retrieval data to estimate surface turbidity directly. The diffuse attenuation coefficient for downwelling irradiance at a wavelength of 490 nm, i.e., Kd490, can be readily derived from the ocean color satellite and has been adopted as a turbidity index for coastal waters (Shi & Wang, 2010). In this paper, the climatological monthly mean Kd490 retrieved from the MODIS Aqua mission is interpolated to the daily surface turbidity (turb) (Morel et al., 2007); turb is then used to obtain the spatially varied transparency (k) by k ¼ 4turb 0.5 /3. PAR sufficiency in the water column at depth z is expressed as I ¼ I 0 e −kz , where I 0 is the solar irradiance at the sea surface provided by the ERA-Interim reanalysis (Dee et al., 2011). The light limitation factor of phytoplankton growth can be written as (Platt et al., 1980) f L ¼ e −p 3 I 1 − e −p 4 I À Á where p 3 and p 4 are constant coefficients.
According to Liebig's law of the minimum (Von Liebig, 1840), growth is dictated by the scarcest resource, and the limiting factor for phytoplankton is given as Taking the influence of temperature (T) into account, the growth rate of phytoplankton is defined as μ ¼ p 5 e p 6 T f ; where p 5 is the maximum growth rate at T ¼ 0 and p 6 is the coefficient of the temperature exponent.

Journal of Geophysical Research: Biogeosciences
Other physiological processes in the model are expressed by simple linear or exponential relations. The maximum respiratory rate of the phytoplankton at T ¼ 0 is p 7 . The natural mortality rate is p 8 . Phytoplankton grazing by zooplankton is modeled by the function z 2 e z3T 1 − e −z4P ð Þ ZP, where z 2 is the maximum growth rate of zooplankton at T ¼ 0, z 3 is the coefficient of the temperature exponent, and z 4 is the grazing rate on the phytoplankton. Considering the aforementioned factors and processes, the coupled physical-biogeochemical governing equation of phytoplankton can be written as follows: ∂ t P þ u∇ H P þ w∂ z P − ∂ H k H ∂ H P − ∂ z k z ∂ z P ¼ μP − p 7 e p 6 T P − p 8 P − z 2 e z3T 1 − e −z4P À Á ZP; where t is time, z is vertical coordinate, u is the horizontal velocity containing the zonal and meridional components, w is the vertical velocity, k H is the horizontal diffusion coefficient, k z is the vertical diffusion coefficient, and ∇ H denotes the operation to obtain the horizontal gradient. The ingestion, excretion, and mortality processes of zooplankton were also taken into account. The zooplankton food assimilation rate is denoted by z 1 . The maximum zooplankton excretion rate is z 5 , and the natural mortality rate is z 6 . Then, the governing equation of zooplankton biomass is given by The source and sink terms for nutrient dynamics are introduced as follows. Taking the Redfield ratio (Redfield, 1934) and the different recycling speeds into account, the nutrient source/sink terms for phosphorus and nitrogen have a ratio of p 2 /p 1 . The physical transport and diffusion coefficients for the two nutrients are the same as that for phytoplankton in Equation A5 and zooplankton in Equation A6. The major loss of dissolved nutrients is caused by uptake by phytoplankton. Part of the mass loss of the phytoplankton, zooplankton, and detritus is assumed to return to the dissolved nutrient pool. The detritus remineralization rate is denoted by d 1 . The governing equations for N n and N p are written as ∂ t N n þ u∇ H N n þ w∂ z N n − ∂ H k H ∂ H N n − ∂ z k z ∂ z N n ¼ −μc n P þ p 7 c n e p 6 T P þ z 5 Z þ d 1 D and ∂ t N p þ u∇ H N p þ w∂ z N p − ∂ H k H ∂ H N p − ∂ z k z ∂ z N p ¼ −μc n P þ p 7 c n e p 6 T P þ z 5 Z þ d 1 D À Á p 2 =p 1 ; respectively.

Table A1
The Parameters in the Biogeochemical Model Detritus, as modeled, is a simplified concept of a relatively massive organic aggregation that includes the excreta and corpses of phytoplankton and zooplankton. It sinks faster than P and Z, which is reflected in the additional detritus sinking rate w d . Only some of the detritus returns to the matter cycle in the euphotic layer via remineralization. Therefore, the control equation for the detritus can be written as The parameters are listed in Table A1.

Data Availability Statement
Observational data and key model data can be accessed via the Mendeley Data repository with the https:// doi.org/10.17632/5ynjsyw5jt.2.