The model domain covers the area of the Bohai Sea, the Yellow Sea, and the ECS, ranging from 24 to 41° N and 117.5 to 131.5° E with a horizontal resolution of 1/18° ( Figure 1 ). The model used here is the same as the one used by Zhang et al. (2019) and Zhang et al. (2021). It consists of two parts. The first part is a physical-biological coupled module. The physical module is based on the Princeton Ocean Model (Guo et al., 2003; Guo et al., 2006). The biological module is based on NORWECOM (Skogen et al., 1995; Skogen and Søiland, 1998). The model is driven by monthly climatological forcing involving wind, freshwater flux, heat flux, nutrient deposition at the sea surface, ocean currents, temperature, salinity, and nutrient concentrations along lateral open boundaries. The main state variables in the module include DIN, dissolved inorganic phosphorus, two kinds of phytoplankton (diatoms, DIA and flagellates, FLA) that are collectively referred to as chlorophyll (CHL), and detritus (DET). The zooplankton is not included in the model because of the bottom-up control in the low-trophic ecosystem of the ECS (Noman et al., 2019). In comparison to zooplankton fecal pellets and other particles, phytoplankton contributes substantially more to the particulate organic carbon (POC) pool in the ECS (Qiu et al., 2018). Biogeochemical processes in the water column involve photosynthesis, respiration, mortality of phytoplankton, and remineralization of detritus, while those in the benthic layer include remineralization and denitrification.

The second part is a tracking module following the methods of Ménesguen et al. (2006). Biological state variables from different external sources are separately calculated in this module rather than as a whole. Only DIN-related variables are tracked here. Because, in contrast to nearshore ECS, offshore ECS is the key region for export production and DIN is a limiting nutrient there. In addition, the input fluxes of the DIN from various sources in the ECS are comparable. Since the Kuroshio contributes a significantly greater amount to DIP than rivers and the atmosphere, DIP-related tracking may result in an expected Kuroshio-dominated result. The DIN from any external sources has the same physical and biogeochemical processes as the total DIN. The governing equations in the tracking module are the same as those in the first module. A complete subset of equations (Equations 1-4) is added to address the DIN-related tracking variables DIN, DIA, FLA, and DET for each source:

where the subscript i represents the DIN from the ith source. The adv and diff represent the physical terms advection and diffusion. The resp, pp, remi, and mort indicate the biological terms respiration, primary production, remineralization, and mortality. The terms for biological processes are separated following the ratio of each source of DIN concentration to the total DIN concentration.

Four external DIN sources in the ECS are considered here, namely, the Kuroshio (K), the Taiwan Strait (T), atmospheric deposition (A), and rivers (R). The nitrogen cycle of each source is processed independently by independent equations. For example, the DIN_K denotes DIN from the Kuroshio. The phytoplankton assimilates DIN_K and produces CHL_K. Then CHL_K becomes DET_K after mortality. DET_K is further decomposed and regenerates DIN_K.

The DIN concentration, phytoplankton, and detritus of the four external sources are set to zero at the beginning of the model run. The input fluxes of DIN from ten main rivers were from published data (Zhang, 1996; Liu et al., 2009). The total riverine DIN flux into the ECS is 1.24 kmol s^-1. The DIN flux from the Changjiang River into the ECS is 0.89 kmol s^-1, which accounts for approximately 80% of the total and plays a dominant role among all the rivers. The climatology data of atmospheric wet and dry depositions was from observations (Zhang et al., 2011). The nutrient concentrations needed for the lateral boundary conditions are climatological mean of observation data for the Taiwan Strait (Huang et al., 2019) and from the Japan Meteorological Agency for the Kuroshio east of the Taiwan Island.

The tracking module is run for over 5 years until it reaches equilibrium, and the outcomes in the final year are analyzed here. The physical-biological coupled model has been fully verified that it can reproduce the climatology features of physical and biological variables in the ECS (Zhao and Guo, 2011; Wang et al., 2019; Zhang et al., 2019; Shen et al., 2021).

We compared the modelled surface climatology temperature and salinity with Chen (2009). The model results are similar with those in Chen (2009). The circulation of the ECS from the model results can demonstrate the most well-known features in current field (Zhang et al., 2019). We also made comparison between the modeled and observed DIN and CHL concentration distributions at PN line and northern ECS (Wang et al., 2019; Zhang et al., 2019). Similar features were identified for DIN and CHL in both our model results and observations. The tracking module has been applied to evaluate the roles of multisource nutrients in primary production (Zhang et al., 2019) and the riverine nitrogen budgets in the ECS (Zhang et al., 2021). These factors are critical for the reliability of the tracking module.

Here, we compare the annual-mean density, velocity normal to the section, DIN and DIP concentrations, and DIN and DIP fluxes at the PN section in Figure 2 with the measurements described in (Guo et al., 2012, Figure 3 ) to validate the model’s capability. The location of the PN section is shown in Figure 1 . The potential density increases from around 23 σ θ at the surface layer to over 27 σ θ at the bottom ( Figure 2A ). Both the range and isopycnal patterns match the measurement. The maximum velocity over 0.8 m s^-1 is located at the surface layer of the shelf break ( Figure 2B ) which is the same as the observation. The velocity core demonstrates the main axis of the Kuroshio, and the velocity decreases from the core to its flanks and depth. The DIN and DIP concentrations generally rise with depth ( Figures 2C, D ). While the concentrations in the upper 200 m are not as low as those from the measurements, the ranges of DIN (5~35 mmol m^-3) and DIP (0.5~2.5 mmol m^-3) concentrations are comparable to the observations. The product of the DIN or DIP concentration and the velocity normal to the section is used to determine the DIN and DIP fluxes ( Figures 2E, F ). In both model and observation results, the largest DIN and DIP fluxes are seen at a depth of around 400 m. The maximum DIN and DIP fluxes (>8.0 and >0.6 mmol m^-2s^-1) have magnitudes that are in agreement with the measurements ( Figure 3 in Guo et al., 2012).

The inventory is defined as the volume-integrated values for a variable concentration (Equation 5). Here we define the area shown in Figure 3 as our target region and calculate the monthly inventories of DIN, CHL, and DET, whose temporal variations are affected by biological and physical processes (Equation 6). Biological processes (denoted as Bio) are responsible for the transformation of DIN, CHL, and DET. The Bio term of DIN includes the remineralization of detritus, phytoplankton photosynthesis, and respiration (Equation 7). The Bio term for CHL involves respiration, mortality, and photosynthesis (Equation 8). The Bio term for DET includes its remineralization and phytoplankton mortality (Equation 9). The physical fluxes (denoted as Phy) include river input or atmospheric deposition (only for DIN), lateral transport across boundaries, and water–sediment flux at the sea bottom. A positive (negative) value of flux represents an increase (reduction) in the variable inventory. The time variation term (denoted by Tendency) refers to the sum of the physical and biological fluxes, which corresponds to the temporal change in inventory. The export production is based on the output flux of PON, which is referred to here as the sum of CHL and DET.

The detailed distributions of the depth-integrated PON concentrations from different sources are shown in Figure 3 . The horizontal distributions of the depth-integrated DIN, CHL, and DET are shown in Supplementary Figures 1–3 .

The PON_R is concentrated along the inner shelf (0~50 m) and plays a dominant role. The higher concentration of PON_R appears in the coastal areas north of the Changjiang Estuary in summer and autumn as the Changjiang Diluted Water ( Figures 3C, D ). The seasonal patterns of PON_R mainly follow those of DET_R ( Supplementary Figure 3 ). The concentration of PON_A is higher on the northern middle shelf (50~100 m) and then decreases southward. It expands further in spring and summer ( Figures 3F, G ). Both PON_A and PON_R show the lowest concentrations in winter ( Figures 3A, E ).

The highest concentration of PON_K appears on the outer shelf (100~200 m) in winter ( Figure 3I ). In spring, the high-concentration area of PON_K narrows down to the southwest of Kyushu ( Figure 3J ). Unlike in winter and spring, the high-concentration area in summer appears northeast of Taiwan Island where the upwelling DIN_K can be converted to CHL_K and the DET_K effectively ( Figure 3K and Supplementary Figures 2, 3 ). The obvious difference between the distributions of PON_K in winter and summer largely corresponds to seasonal variations in the Kuroshio intrusion, which flows onshore from the surface in winter and from the bottom in summer (Yang et al., 2012; Zhang et al., 2017). The pattern of PON_K in autumn is similar to that in summer but with a lower concentration ( Figure 3L ).

The concentration of PON_T is lowest among the four kinds of PON except in summer when the CHL_T and DET_T are largely produced ( Supplementary Figures 2, 3 ) and has a maximum value on the southern middle shelf ( Figure 3O ). The high-concentration areas of PON_T mainly follow the Taiwan Warm Current and its branches (Liu et al., 2021).

The seasonal variations of inventories of DIN, CHL, and DET are quite different between different sources ( Figure 4 ). They are affected by biological and physical processes (bars in Figure 4 , Equation 5). The inventory of DIN_R shows a peak value in March and a valley value in July ( Figure 4A ). The Bio term controls the variation of DIN_R and makes the inventory decrease in spring and summer. The enhancement of DIN_R from July to the next March is attributed to the reduced consumption by biological processes and gain from the river input in summer. The inventory of DIN_A shares similar seasonal patterns with DIN_R ( Figure 4D ). The correlation coefficient between DIN_A and DIN_R inventories is 0.97, which indicates a significant relationship. There was no significant correlation between the other DIN inventories.

The inventory of DIN_K plays a leading role in the quantity and seasonal amplitude ( Figure 4B ). It is highest in winter and lowest in April-May. The Bio term of DIN_K shows a negative value all year round, suggesting that the total biological processes spend the DIN_K. The seasonal variations of DIN_K follow its Bio term. The DIN_K inventory decreases from January to April due to the Bio term, then it increases sharply after May resulting from the physical flux.

The inventory of DIN_T is relatively low and changes weakly ( Figure 4C ). The DIN_T inventory reaches its maximum in October and its minimum in May-June. Unlike the Bio terms of the other sources with only one peak in spring, the Bio term of DIN_T has two peaks in April and July while the latter one is larger. The primary production supported by DIN_T mainly happens in the middle shelf resulting from the large volume transport and corresponding DIN_T input from the Taiwan Strait in summer (Zhang et al., 2019). The large DIN_T input in July makes its tendency term positive and increases the DIN_T inventory, even though the DIN_T consumption by biological processes also peaks.

The inventories of CHL_R and CHL_A change in phase with respect to DIN_A and DIN_R. Both of them peak in May and dissipate away in winter with a correlation coefficient of 0.95 ( Figures 4E, H ). The strong Bio terms of CHL_R and CHL_A dominate the time variations of their inventories. The CHL_K is still the largest among all the CHL at any time with a maximum in March and a minimum in autumn ( Figure 4F ). However, unlike CHL_R and CHL_A, the seasonal variations of the CHL_K inventory do not fully follow the Bio term, which is positive the whole year, even in winter. The strong Bio term of CHL_K gives a positive tendency term from January to March. But after April, the combined effects of enhanced Phy term and weakened Bio term make the tendency term of CHL_K below zero. The CHL_T inventory also shows low value in winter, and peaks in August ( Figure 4G ). It is significantly correlated to CHL_R and CHL_A with correlation coefficients of 0.91 and 0.77. The Bio term of CHL_T dominates the seasonal variations of tendency term, unlike the DIN_T situation.

The Bio term for DET comprises the mortality of phytoplankton and remineralization to DIN (Equation 9). The DET Bio term is positive most of the time, its peak time agrees with that of mortality of CHL (figure is not shown). The DET inventory from a specific source is markedly larger than the CHL inventory from the same source. The amplitudes of the four kinds of DET inventories are not much different. The superiority for Kuroshio source in DIN_K and CHL_K no longer exists in DET_K ( Figure 4J ). In March, the inventory of DET_K reaches its maximum, while the other three inventories of DET are the lowest. In the following June, July, and August, the maximum of DET_A, DET_R, and DET_T inventories appear one after another. The DET_R inventory is significantly large in summer and autumn due to the Bio term and then drops down with gradually increasing Phy term ( Figure 4I ). The Bio term of DET_R has the largest value in May, which is one month after the peak of CHL_R. In wintertime, the Bio term of DET_R is negative, indicating a stronger remineralization process than mortality. The Bio term of DET_T is largest in July when the CHL_T are produced in large numbers and a high temperature causes more mortality of phytoplankton ( Figure 4K ). The tendency term is therefore significantly above zero in spring and summer, which causes the DET_T inventory to increase and then drops to below zero in autumn. Comparing of these four kinds of DET, it is found that DET_R, DET_A and DET_T are correlated two by two, while DET_K is not correlated with any of the remaining three DET.

The annual means of each flux and inventory are shown in Supplementary Table 1 . In the climatological model result, the inventories of the variables change little and the physical and biological fluxes are equivalent on a yearly basis. Over such a time scale, physical fluxes raise DIN inventories while biological activities reduce them. The CHL and DET, on the other hand, are produced by biological activities and exported by physical fluxes. Returning to seasonal variations, the peak times of Phy terms of CHL and DET lag those of Bio terms by one or a few months. That is because the outputs of the CHL and DET occur after a certain amount of accumulation of their inventory.

The seasonal variations of each physical flux are shown in Figure 5 . The physical fluxes include the DIN from rivers (Riv) and atmospheric deposition (Atm), the benthic fluxes (Bot), and the lateral transports across the Taiwan Strait (TAS), the 200-m isobath (200m), the Tsushima Strait (TUS) and the 34.7° N section (34N). The return of DIN from the sediment to the water column denotes the Bot flux of DIN, while that for CHL and DET is the net flux of resuspension and sinking.

The main physical flux of DIN_R is river input ( Figure 5A ). It is highest in summer like the river discharge. While benefiting from the large input from the river, the total inventory of DIN_R increases slowly after spring ( Figure 4A ). Meanwhile the high export of DIN_R to Tsushima Strait begins. The output flux of DIN_R through the Tsushima Strait is larger than other fluxes in the summer half year. The patterns of DIN_A fluxes are similar to those of DIN_R ( Figure 5D ) but it has a larger flux across 200-m isobath. The variation ranges of DIN_K fluxes are the largest among all the DIN sources ( Figure 5B ). The fluxes through the 200-m isobath and Tsushima Strait represent the main pathway of DIN_K, which intrudes across the 200-m isobath and exits via the Tsushima Strait. Therefore, their seasonal variations are almost in phase. Both reach a maximum in autumn when the Kuroshio onshore intrusion is highest (Guo et al., 2006; Zhao and Guo, 2011). The most prominent fluxes of DIN_T are input from the Taiwan Strait ( Figure 5C ).

From June to December, the CHL_R flux through the Tsushima Strait is distinct with a maximum in autumn ( Figure 5E ). The Bot flux of CHL_R is largest in summer with weak seasonal variation. The fluxes of CHL_K across the 200-m isobaths and the Tsushima Strait reach a maximum in April after the peak of the CHL_K inventory in March ( Figures 4F , 5F ). The flux through the Tsushima Strait is a little larger than that across 200-m isobaths throughout the year. The input of CHL_T from the Taiwan Strait is the dominant flux like the DIN_T situation ( Figures 5C, G ). It has two peaks in May and August. The maximum physical flux of CHL_A in summer results from the export through Tsushima Strait and 200-m isobath ( Figure 5H ).

The Bot terms show their important roles in all DET fluxes. The striking flux of DET_R is in the Bot terms ( Figure 5I ). The Bot term of the DET_R increases from April when the DET_R greatly increases. The largest two fluxes of DET_A are the Bot term and export to the Tsushima Strait ( Figure 5L ). Both the Bot term and export to the Tsushima Strait reach the maximum in June. The peak months of Bot fluxes of DET_R and DET_A are the same as those of their inventories. The Bot term of DET_T also plays a dominant role and it reaches the maximum one month after DET_T inventory peaks ( Figure 5K ). The seasonal variations in the DET_K fluxes are similar to those in the CHL_K fluxes except for a larger Bot flux ( Figure 5J ). The export fluxes of DET_K and CHL_K through three different exits (Bot, TUS, and 200m) show a comparable peak time in March or April.

In terms of long-term export production, the biogenic particles fluxes that are removed through the sediment and lateral boundaries to the open ocean are characterized as export production on the continental shelf. The main pathways of the export production over the ECS shelf are via the sediment and the Tsushima Strait to the Japan Sea and across the 200-m isobath to the northwest Pacific Ocean.

In the ECS, the particulate organic carbon to PON ratio approaches the Redfield ratio of 6.63, with a mean values ranging from 4 to 8 (Hung et al., 2000; Zhu et al., 2006). The PON fluxes may be utilized to estimate the organic carbon flux in this situation. The overall export production over the ECS is 6.83 kmol N s^-1 (=17.16 Tg C yr^-1) using a C/N ratio of 6.63, and ranging from 10.35-20.71 Tg C yr^-1 when considering a C/N ratio of 4~8. The export production is composed of net-Bot, TUS, and 200 m fluxes ( Table 1 ). This value is slightly lower than the atmospheric CO₂ absorption of 23.3 ± 13.50 Tg C yr^-1 (Guo et al., 2015), suggesting that export to the open ocean and seafloor can reflect efficient carbon sequestration. The benthic flux of 1.40 kmol N s^-1 (44.15×10⁹ mol yr^-1) is less than the burial flux of organic nitrogen at 75×10⁹ mol N yr^-1 (Chen and Wang, 1999) because the benthic flux here considers only ocean-generated nitrogen and excludes the input of terrestrial organic nitrogen, because they only comprise small portions of the total Changjiang River nitrogen input flux (Zhang et al., 2003; Kwon et al., 2018). The 200 m flux of 1.95 kmol N s^-1 (=4.89 Tg C yr^-1) is between the reported POC output flux of 0.25~1.7 Tg C yr^-1 (Deng et al., 2006; Jiao et al., 2018; Yuan et al., 2018) and 12.61 Tg C yr^-1 (Chen and Wang, 1999). The net-Bot flux and the 200 m flux explain 20% and 29% of the total export production, respectively, while the TUS flux accounts for 51%. Consequently, approximately 80% of the total export production leaves by lateral transport. This value is consistent with the percentages of off-shelf transport and sediment burial of 80% and 20%, respectively, in eastern North America (Najjar et al., 2018), 70% and 30%, respectively, on the whole North American shelf (Fennel et al., 2019), and 60~100% and 40~0%, respectively, in northwest European continental shelf seas (Legge et al., 2020). However, this ratio varies among the PONs supported by four sources of nutrients, which is covered in the following subsection.

The e-ratio is defined as the ratio of export production to primary production, which indicates the export efficiency of organic matter output from the generated carbon. The total primary production and export production are 39.9 and 6.83 kmol N s^-1, respectively, forming an e-ratio of 0.17. Our value is comparable to the result from the carbon budget based on a box model used by Chen and Wang (1999). They predicted an f-ratio of 0.15, which was supposed to balance the e-ratio in a steady state. Zuo et al. (2016) found that approximately 12% of the generated particulate nitrogen was exported between May and October in the ECS, which is less than our annual-averaged e-ratio due to the relatively high temperature and remineralization during their observation period. Based on corrected trap-collected POC fluxes, Hung et al. (2016) estimated an e-ratio of 0.59 nearshore and 0.16 offshore. The high e-ratio in the nearshore area demonstrated a significant vertical flux across the euphotic layer depth rather than efficient carbon export.

The export production values of PON from various sources have a wide range of features. The export production of PON_K is 4.40 kmol N s^-1, accounting for almost 64% of the total export production. The export production of PON_A is 1.10 kmol N s^-1, yet it accounts for 16%. The export production of PON_T accounts for 12%, while that of PON_R makes up only 8%. The PON_K net-Bot flux (0.51 kmol N s^-1) is the largest, accounting for nearly 36% of the overall net-Bot flux. The net-Bot fluxes of PON_R and PON_A are equivalent and account for approximately 23% of the total flux. The burial flux is larger in the coastal and shelf-break areas (Deng et al., 2006). Thus, the PON_T, which is mainly concentrated on the middle shelf, has a slightly lower net-Bot flux. The TUS flux of the PON_K has a fraction of 64%, whereas that of the PON_R has a fraction of 5%. The PON_K flux across the 200-m isobath accounts for 85% of the total, and the PON_R flux accounts for only 0.1%.

The primary export pathway for the PON_R is the sediment. The net-Bot flux of PON_R accounts for over 60% of the total export production of rivers. This ratio is disproportionate larger than the 20% for the total PON. The PON_R flux across the 200-m isobath is negligible. On a yearly basis, the impact of DIN_R and PON_R on the Kuroshio region is insignificant. The excess DIN_R and a small amount of PON_R may be exported through the Tsushima Strait (Isobe and Matsuno, 2008; Zhang et al., 2021). The PON_K flux is the largest through the Tsushima Strait, which benefits from the Kuroshio Branch Current west of Kyushu (Hsueh et al., 1996; Lie et al., 1998; Isobe, 2000) and the outflow through the Tsushima Strait to the Japan Sea (Isobe et al., 2002; Teague et al., 2003). The total lateral transport of PON_K (TUS and 200m) accounts for nearly 90% of the export production. The Tsushima Strait is also the principal export pathway for PON_T and PON_A, and the sediment is the second important export pathway. About 30% of the overall export production of PON_T and PON_A are made up of their net Bot fluxes, respectively. PON_T and PON_A are mostly found on the middle shelf ( Figure 3 ), demonstrating a distinction between coastal PON_R and oceanic PON_K.