Our study area is located in the NSCS and covers the Peral River plume ( Figure 1 ). As the third-largest river in China, the Pearl River discharges into the NSCS and delivers approximately 3.3 × 10¹¹ m³ of freshwater per year (Li et al., 2017). The NSCS is under the influence of the East Asia Monsoon, which makes our study area experience frequent cyclonic and anti-cyclonic circulation, resulting in diluted Pearl River water reaching our study area through the continental shelf (Morimoto et al., 2000; Lin et al., 2010). Thus, the area-of-interest is jointly affected by the water masses of SCSWC, SCSKB, and a combination of the diluted Pearl River water and the GDCC. Under these circumstances, the nutrient distribution in our study area can be influenced by shelf water in some specific periods.

Two cruises were conducted in the NSCS during June 2015 and March 2017 while onboard the R/V Nanfeng. Surface (~1 m) Ra samples of approximately 200 L were collected using a submersible pump at 21 stations and 16 stations in June 2015 and March 2017, respectively ( Figures 2B, C ), while subsurface samples (~100 L) were taken directly from an onboard Conductivity–Temperature–Depth (CTD) rosette. After collection, water samples were immediately passed through a column that was filled with approximately 20 g of MnO₂-impregnated acrylic fiber at a flow rate of 0.5 L min⁻¹ to enrich Ra isotopes (Moore & Reid, 1973). The seawater temperature and salinity were measured in situ using a Sea-Bird CTD (SBE 911plus, Sea-Bird Electronics, Inc., USA). Notably, the Ra samples in the Pearl River plume (K1–K5) were collected while the vessel was in transit, so the seawater temperature and salinity were measured by a portable salinometer with multiple parameters (Germany, multi350i). For each Ra sample, approximately 60 ml of samples were collected for the dissolved nutrients after being filtered through a 0.4 µm pore-size polycarbonate filter (Whatman, USA), and then stored frozen for laboratory analysis.

Upon returning to the laboratory, the Mn fibers were ashed at 800°C for 8 h, homogenized, and loaded into a plastic vial sealed with an epoxy sealant for measurement. The details are shown in Liu et al. (2021a). Dissolved nutrient concentrations were determined using a QuAAtro Continuous–Flow Automatic Analyzer (SEAL Analytical GmbH, Norderstedt, Germany), and the analytical precision of NO 2 − , NO 3 − , NH 4 + , DIP, and DSi were all better than 5%, and the detection limits were 0.01, 0.01, 0.02, 0.01, and 0.04 μmol/L respectively (Wu N. et al., 2021). The concentration of dissolved inorganic nitrogen (DIN) is the sum of NO 3 − , NO 2 − , and NH 4 + .

The distributions of surface salinity and temperature in both two seasons are shown in Figure 2 . During the observation period of June 2015, surface salinity ranged from 31.42 to 34.10 with the lowest values occurring in a transect of stations J9–J12, and the highest salinity was observed in the western and eastern parts of our study area ( Figure 2A ). Temperature ranged from 29.36 to 31.08°C with no clear distribution characteristic ( Figure 2B ). While in March 2017, surface salinity ranged between 33.33 and 34.70 and showed a significant differentiation trend, in which high salinity stations were all located on the eastern side and lower salinities were observed on the western side ( Figure 2C ). Temperature ranged from 24.38 to 26.38°C and had the opposite distribution trend as salinity, with the exception of station M1 ( Figure 2D ).

The distributions of surface ²²⁶Ra, ²²⁸Ra activities, and ²²⁸Ra/²²⁶Ra ratios in our study areas are shown in Figure 3 , with obvious spatial variations. Actually, surface ²²⁶Ra, ²²⁸Ra activities, and ²²⁸Ra/²²⁶Ra ratios show similar patterns in both seasons. Specifically, in June 2015, our observed activities of ²²⁶Ra and ²²⁸Ra ranged from 77 to 120 dpm/m³ and from 111 to 305 dpm/m³, respectively, while the ²²⁸Ra/²²⁶Ra ratios were between 1.4 and 2.6. The lowest values of ²²⁶Ra, ²²⁸Ra, and the ratio occurred in the eastern region, which corresponded to the highest salinity. Similar to salinity and temperature, obvious differences in ²²⁶Ra, ²²⁸Ra and the ratios distributions were also observed between the eastern and western regions, with low values occurring on the eastern side and high values occurring on the western side.

From the Pearl River plume to the slope of the NSCS, values of surface salinity, Ra, and nutrients showed remarkable variations ( Figure 4 ). Salinity significantly increased from 17.40 to 34.1, and ²²⁶Ra and ²²⁸Ra decreased from 388 to 96 dpm/m³ and from 610 to 210 dpm/m³, respectively. Before entering our square sampling area, good relationships with the distances from the coasts supported the distributions of ²²⁶Ra and ²²⁸Ra. The concentrations of DIN ranged from 0.14 to 99 μmol/L, PO 4 3 − concentrations ranged from 0.03 to 0.45 μmol/L and DSi concentrations ranged from 1.8 to 81 μmol/L. Nutrient concentrations both showed clear decreasing trends with the gradients (Δnutrient/Δx) of −0.56 ± 0.30, −0.0022 ± 0.0013, and −0.41 ± 0.27 μmol/L/km for DIN, PO 4 3 − and DSi, respectively. Along this transect, values of surface Ra and nutrients changed little over our investigating area (stations J8–J5), with averages of 93 ± 8.6 dpm/m³, 190 ± 36 dpm/m³, 0.10 ± 0.015 μmol/L, 0.013 ± 0.0047 μmol/L, and 1.7 ± 0.62 μmol/L for ²²⁶Ra, ²²⁸Ra, DIN, PO 4 3 − , and DSi, respectively, which were much lower than those on the transit pathway.

The vertical profiles of hydrological parameters, ²²⁶Ra, ²²⁸Ra activities, and ²²⁶Ra/²²⁸Ra ratios are shown in Figure 5 . In the upper 400 m of station J6, temperature, density, and salinity largely varied. The activity of ²²⁶Ra decreased from 98 to 85 dpm/m³, the activity of ²²⁸Ra decreased from 177 to 43 dpm/m³, and the ²²⁶Ra/²²⁸Ra ratio decreased from 1.8 to 0.47. However, from 1,000 m to near the bottom, density and salinity stayed nearly constant, while ²²⁶Ra and ²²⁸Ra activities both showed increasing trends and maximum values occurred near the bottom. Especially for ²²⁶Ra, the highest activity was observed near the bottom. The ²²⁶Ra/²²⁸Ra ratios below 1,000 m were almost constant at 0.019 but were noticeably lower than those observed in the upper 400 m. Similar patterns were also distributed at station J8, which suggested that the temperature, density, and salinity changed very little below 70 m, and ²²⁶Ra, ²²⁸Ra activities, and ²²⁶Ra/²²⁸Ra ratios showed downward trends from the bottom up.

When plotting ²²⁸Ra/²²⁶Ra ratios on the potential temperature–salinity diagram, it showed that the mixing of the water masses in the NSCS distributed their unique signals of ²²⁸Ra/²²⁶Ra ratios and salinity and has obvious seasonal variations ( Figure 6 ). Specifically, relatively low salinity occurred in June 2015, indicating that the Shelf Water (SHW) had a significant impact on the region of interest, which may be driven by the joint contribution of Pearl River freshwater and mesoscale eddies (He et al., 2016; Zhang et al., 2019a). Besides, our study areas were affected by the SCSWC and SCSBK in both seasons, as illustrated in Figure 6 .

Due to fact that the effect of biogenic particles on Ra activity can be neglected, surface ²²⁶Ra and ²²⁸Ra activities were controlled only by water mixing and radioactive decay (G. Wang et al., 2021b). As a result, a mixing model of multi-end members based on Ra isotopes could be set up to quantify the proportion of water masses in the NSCS. The contributions of SHW, SCSWC, and SCSBK to the surface water in the area of our interest in June 2015 can be estimated by solving the following equation:

where f refers to the fraction of each water mass, S is salinity, ²²⁶Ra is the activity of the indicated Ra isotopes, and the subscripts SHW, SCSWC, SCSBK, and obs represent the SHW end-member, the SCSWC end-member, the SCSBK end-member, and the observed values of an individual sample, respectively. Here, we took the station K3 as the SHW end-member, in which the salinity was 30.30 and the ²²⁶Ra activity was 134 ± 13 dpm/m³, so ²²⁶Ra could be totally desorbed at such high salinity. Besides, the salinity and ²²⁶Ra activity in the SCSWC end-member were 33.81 and 117 ± 6.9 dpm/m³, respectively, and in the SCSBK end-member were 34.69 and 42 ± 2.7 dpm/m³ (Nozaki et al., 1989), respectively, and they were both in the respective transit pathways of the two water masses. Thus, the fraction of each water mass in June 2015 can be obtained using Eq. (1).

While in March 2017, signals of two water masses were observed, and we applied a two-end-member mixing model developed by Moore et al. (1986) to access the fractions, which can be written as follows:

where AR_obs denotes the ²²⁸Ra/²²⁶Ra ratio in our observed samples. Note that even if there were only two water masses, salinity was still used to correct for evaporation and precipitation. Meanwhile, we chose the ²²⁸Ra/²²⁶Ra ratio to build the mixing model rather than the ²²⁶Ra or the ²²⁸Ra alone, because it was expected to reduce the effects of biological usage and physical interaction of a single isotope on the calculation (Kawakami & Kusakabe, 2008; Chen et al., 2010). For the SCSWC end-member during the period, we applied the values of 33.70 for salinity, 76 ± 3.7 dpm/m³ for ²²⁶Ra and 197 ± 11 dpm/m³ for ²²⁸Ra as collected in similar seasons (Tan et al., 2018). Because the ²²⁸Ra/²²⁶Ra activity ratio in the Kuroshio water stayed relatively stable over time (Wang et al., 2021b), the Ra activities in the SCSBK end-member were used as previously described. Combining the salinity and the ²²⁸Ra/²²⁶Ra ratios measured during our observation, the individual fraction of each water mass can be estimated from Eq. (2).

The fraction of the SHW in the study area of interest ranged from −0.01 to 0.73, with an average of 0.23 ± 0.26. The pattern of the SHW fraction was similar to salinity, with a high fraction corresponding to low salinity and vice versa ( Figure 7 ). The maximum fraction of the SCSBK in June 2015 was observed in the southeast of the study area and a considerable fraction also appeared on the western side, indicating that the Kuroshio water could reach as far west as 115°E in the surface water of the NSCS under the possible influence of cyclonic circulation (Gan et al., 2016). The averaged SCSBK fraction was estimated to be 0.25 ± 0.16, which was highly in agreement with the result of 0.23 ± 0.11 derived by Wang et al. (2021b), which was both obtained in the summer and based on Ra isotopes. The SCSBK fraction in March 2017 ranged from 0.01 and 0.91, with an average of 0.57 ± 0.32 and was significantly higher than that in summer. The comparison followed a typical pattern because the strongest intrusion was usually shown to occur in the wintertime and weakened toward summer (Hsin et al., 2012; Nan et al., 2015). The western parts of the study area were heavily affected by the SCSWC in both two seasons, and the mean fractions were 0.52 ± 0.27 and 0.43 ± 0.32 for June 2015 and March 2017, respectively. Generally, the formation and expansion of the SCSWC are dynamically related to the SCSBK (e.g., Xue et al., 2004; Yu et al., 2021), so the intrusion of the SCSBK into the NSCS can influence the SCSWC proportion, which resulted in that higher fraction of the SCSWC occurring in June 2015. Overall, the fractions of the multi-water masses distributed significant seasonal variations.

The fractions of water masses obtained by mixing models (Eqs. 2, 3) are generally sensitive to end-member variations. In this study, we also conducted an uncertainty analysis to evaluate our end-member choices. In June 2015, the water mass fractions were the most sensitive to ²²⁶Ra activity variations in the SCSWC end-member. Specifically, a 10% increase in the activity of ²²⁶Ra would increase in the fraction by 7.9% (SHW) to 28% (SCSBK), while 10% increases in the SHW and SCSBK end-members would only cause fraction variations by 5.0% (SHW) to 18% (SCSBK) and 1.7% (SHW) to 6.3% (SCSBK), respectively. In March 2017, similarly, a 10% increase in the activities of ²²⁶Ra or ²²⁸Ra in the SCSWC end-member would cause a fractional increase of 11% (SCSBK) to 17% (SCSWC), but only 1.1% (SCSBK) to 4.2% (SCSWC) caused by the activity variations of ²²⁶Ra or ²²⁸Ra in the SCSBK end-member. The quantitative uncertainty analysis indicated that our considered end-members were suitable to build the mixing models. Additionally, ²²⁶Ra and ²²⁸Ra measurement errors can also be involved in our mixing models and lead to variations in water mass fractions. Here, the measurement errors of ²²⁶Ra and ²²⁸Ra were both 3.6–12%, and we used a maximum measurement error of 12% to conduct the uncertainty analysis. The results showed that a 12% increase in observed ²²⁶Ra activities would cause a change in the water mass fraction by 19% (SHW) to 68% (SCSBK) in June 2015 and 14% (SCSBK) to 18% (SCSWC) in March 2017, suggesting a significant influence of Ra measurement errors on the uncertainty of water mass fraction, and the greater measurement error would usually contribute to greater uncertainty (Wang et al., 2021b).

In surface water, horizontal mixing is typically several orders of magnitude higher than vertical mixing, indicating that the term of downward mixing for ²²⁸Ra could be neglected (Hsieh et al., 2021). Therefore, except for water mixing and radioactive decay, there was no additional input or removal of ²²⁸Ra in the transit pathway and the study area of our interest. Assuming a steady state, a common one-dimensional ²²⁸Ra advection–diffusion model was set up to measure diffusion coefficients and advection rates, and the formula was expressed as follows (Moore, 2015):

where K_x denotes the horizontal eddy diffusion coefficient, ω is the advection velocity, x is the distance from the coast, A_ex the excess ²²⁸Ra activity and λ is the ²²⁸Ra decay constant. Actually, a boundary condition was applied in Eq. (3), namely, the ²²⁸Ra activity was zero at an infinite distance (x→∞). However, this condition is not usually valid within the relatively small offshore distance scale (<50 km) (Moore, 2000), just like in our study region, not to mention that our measured surface ²²⁸Ra activity is far away from zero. As a consequence, we also followed the suggestion of Moore (2015) and Hsieh et al. (2021) to define the excess ²²⁸Ra by subtracting the background value in the center of the NSCS and then applying it to build the advection–diffusion model.

In this study, due to the significant excess ²²⁸Ra distributions before and after the shelf break ( Figure 8 ), two scenarios were considered in our horizontal ²²⁸Ra estimations, one was mixing only (ω = 0) and the other was advection only (K_x = 0). Therefore, with the boundary conditions of A _{ex_0} = A ₀ − A _bg at x = 0 and A _ex = 0 at x = ∞, Eq. (3) can be solved for diffusive mixing only as:

And for advection only as:

where A _{ex_0} is the ²²⁸Ra_ex activity at x = 0, X _1/2 is the distance at which A _ex = 0.5A _{ex_0} and T _1/2 is the half-life of ²²⁸Ra.

Therefore, when plotting the surface excess ²²⁸Ra activities versus the distance from the coast, the best exponential fit was observed before the shelf break, with a gradient of 0.018 ± 0.0063 ( Figure 8 ). So based on Eq. (4), the horizontal diffusion coefficient K_x was estimated to be (1.2 ± 0.79) × 10⁵ cm²/s. The result was significantly lower than that obtained in the Pearl River estuary region (<80 km from the coast) of 4.7 × 10⁶ cm²/s (Wang et al., 2021c), but was a bit higher than that mentioned in Liu et al. (2020) of 5 × 10⁴ cm²/s obtained by numerical modeling. In this region, we only use the mixing model because of the Pearl River intrusion that results in a strong mixing signal. For the advection after the shelf break, according to Eq. (5) and the distinct linear fit between excess ²²⁸Ra activity and distance from the coast ( Figure 8 ), the X _1/2 was 445 ± 288 km (from station J8), thereby accessing advection velocity ω of 0.25 ± 0.16 cm/s. Our ²²⁸Ra-derived advection velocity was lower than the result from the track of the drifter during the same cruise and distributed a consistent direction (Chen et al., 2016), but was comparable to other surface regions of the NSCS (e.g., Chou et al., 2006; Liu et al., 2020). However, the surface water of the NSCS is occasionally controlled by eddies, which may result in considerable spatio-temporal variations of the advection velocity values.

In the vertical profiles, ²²⁸Ra activities were distributed, with ²²⁸Ra concentrations in surface and near-bottom waters considerably higher than those in intermediate waters, indicating that the surface-mixed layer could supply the ²²⁸Ra down into the subsurface and that ²²⁸Ra in near-bottom water could be transported upward under the influence of sediment diffusion (Cai et al., 2002; Van Beek et al., 2007). Considering that ²²⁸Ra for vertical mixing calculations is mixed horizontally far away from the coast, the one-dimensional model can also be applied to obtain the vertical diffusion coefficients (K _z) downward and upward and was expressed as (Hsieh et al., 2021):

where A _ez is the activity of excess ²²⁸Ra in the depth profile and z is the water depth. Like the horizontal mixing, we used the similar boundary conditions of A _{ez_0} = A ₀ − A _bg at x = 0 and A _ez = 0 at x = ∞, then solved Eq. (6) as follows:

where A _{ez_0} is the ²²⁸Ra_ex activity at z = 0 (the mixed layer or near-bottom water). For a conservative calculation, the lowest ²²⁸Ra activity in the vertical profile of the station J8 was the background value, which was 22 ± 2.6 dpm/m³. So we plotted the ²²⁸Ra_ex activity versus water depth, as shown in Figure 9 , and both showed significant exponential fitting coefficients in the upper and lower layers. Note that the average activity and depth of the upper two samples of Figure 5 were used in the mixed layer.

In the upper 1,000 m, the best-fit exponential curve gradient (a) was 0.0094 ± 0.0036, leading to a vertical diffusion coefficient (K _z) of 0.43 ± 0.33 cm²/s down to the subsurface. Our result was highly consistent with the value of 0.23 cm²/s using ²²⁸Ra, estimated by Cai et al. (2002) in the upper 300 m, and within the range of 0.037 to ~10 cm²/s obtained from the results reported in the upper layer of the South China Sea (Li et al., 2016; Shang et al., 2017; Shih et al., 2020). As mentioned above, Ra diffusing across the sediment–water interface usually dominates Ra activity in the lower layer of the water column, accompanied by good exponential regression for ²²⁸Ra_ex and simulated a fitting curve of 0.0014 ± 0.00038 ( Figure 9 ). So the upward vertical diffusion coefficient was revealed as 18 ± 9.9 cm²/s according to Eq. (7). To our knowledge, our estimate was the first exploration of the Ra-derived upward vertical diffusion for the South China Sea but was comparable with the results in other Chinese seas (Liu et al., 2010; Su et al., 2013b). As a whole, in the upper layer, the vertical diffusion downward was 4–5 orders of magnitude lower than the surface horizontal mixing, confirming our above assumption when calculating horizontal mixing. While in the water column, the upward vertical diffusion from the bottom sediment was significantly higher than the vertical diffusion down to the subsurface.

Because of water mixing, the surface waters on the slope of the NSCS receive nutrient supplies from multiple sources. For example, shelf water surprisingly intruded into the study area of interest in June 2015, and it was bound to carry nutrients from the coastal zones. Combined with the estimated horizontal diffusion and gradients of nutrient concentrations ( Figure 4 ), the horizontal nutrient fluxes were assessed to be (5.6 ± 4.9) × 10² mmol/m²/d for DIN, 2.2 ± 2.0 mmol/m²/d for DIP, and (4.1 ± 3.9) × 10² mmol/m²/d for DSi. Previously, Tan et al. (2018) reported that SGD delivered nutrient fluxes of 7.9–19 mmol DIN/m²/d, 0.037–0.079 mmol DIP/m²/d, and 12–28 mmol DSi/m²/d into the NSCS, which were much lower than the ²²⁸Ra-derived horizontal diffusive nutrient fluxes. Meanwhile, our estimated horizontal DIN flux was much higher than that obtained by Li et al. (2018) of 0.2–3.6 mmol/m²/d (only for NO 3 − ). In fact, our estimates usually provide integrated fluxes of nutrients, not only considering all possible inputs, such as the Pearl River, SGD, and sediments from the coastal zones to the slope of the NSCS, but also the eddy-entrained Pearl River plume into the NSCS, which was also observed during our investigation (He et al., 2016), suggesting that the horizontal diffusion mixing could carry a large amount of nutrients to our study area of interest. While in the slope region after the shelf break, the mixing process also occurred for nutrients due to advection. The term was obtained by multiplying the advection velocity with the concentrations of nutrients in the initial advective waters (Hsieh et al., 2021). Here, we used nutrient concentrations at station J8 as the initial advective water, which were 0.12, 0.02, and 0.69 μmol/L for DIN, DIP, and DSi, respectively. Thus, the estimated ²²⁸Ra-derived advective nutrient fluxes in the slope region were 25 ± 16 mmol/m²/d for DIN, 4.2 ± 2.7 mmol/m²/d for DIP, and 146 ± 195 mmol/m²/d for DSi. For the vertical water column, the nutrient concentrations showed decreasing trends from near-bottom to the surface waters, and the upward vertical diffusion plays an important role among it. Considering the decreasing trends with the gradients of 0.017 ± 0.0047 μmol/L/m for DIN, 0.0011 ± 0.00032 μmol/L/m for DIP, and 0.095 ± 0.014 μmol/L/m for DSi, giving the upward vertical nutrient fluxes of 2.7 ± 1.6, 0.18 ± 0.11, and 15 ± 8.4 mmol/m²/d for DIN, DIP, and DSi, respectively, in which the vertical upward DIN flux was much lower than the result of 45 mmol/m²/d during the impacts of mesoscale eddies in the South China Sea (Guo et al., 2015), but was consistent with the previous estimate of 0.66 mmol/m⁻²/d (only for NO 3 − ) in Cai et al. (2002) and the ranges of 0.11 to 1.54 reported in Shih et al. (2020) and the references therein. The mixing associated DIP and DSi fluxes were scarcely reported in previous studies of the NSCS, however, our comparable DIN fluxes with others also gave us confidence in the DIP and DSi fluxes.

To better understand the nutrient dynamics in the NSCS, a very simple diagram of nutrient sources was built in the region of interest, as shown in Figure 10 . It is obvious that the dominant nutrient source in the study area was identified as horizontal diffusion, namely, the input from the coastal zone by shelf water, without considering internal regeneration. Besides, vertical upward mixing appeared to be a more important source of supplying nutrients to the upper layer compared to atmospheric deposition (Gao et al., 2020b; Wu et al., 2018). As mentioned above, because of the multiple sources and eddy-entrained Pearl River plume, shelf water unusually transported materials to the slope of the NSCS, not only the nutrients raised in this study, but also the carbon (Zhang et al., 2020) and trace metals (Wang et al., 2019; Zhang et al., 2019b), suggesting that the intrusion of shelf water played a vital role in the slope of the NSCS during our cruise. In our nutrient diagram, even though advection generally does not provide nutrients in the study area, it transports nutrients and mixes them with other waters. Note that we did not consider the biological effects when conducting our estimates of nutrient fluxes, and further studies are needed to understand the combined effects of water mixing and biological uptake on nutrient cycles in the ocean.