Zhanjiang Bay is located in the northwest of the South China Sea and is a typical bay with a small mouth and large belly ( Figure 1 ). The connection between Zhanjiang Bay and the South China Sea is mainly through a narrow channel with a width of approximately 2 km (Zhang et al., 2020b). Zhanjiang Bay is located in the subtropical monsoon climate zone, with a rainy season from April to September, and less rainfall from November to February of the following year (Chen et al., 2019). There are many industrial zones, agricultural zones, aquaculture zones, ports, and densely populated areas along the coast of Zhanjiang Bay. A large amount of industrial and agricultural wastewater and domestic sewage are discharged into the bay, bringing a large amount of nutrients and organic matter to Zhanjiang Bay, which has a certain impact on the ecological environment of Zhanjiang Bay (Zhang et al., 2023). Previous studies have shown that the degree of eutrophication in the water of Zhanjiang Bay is gradually becoming severe (Zhang et al., 2020b).

The sample collection was conducted in Zhanjiang Bay in May and September of 2016, respectively, with 23 surface water samples collected in May and 29 surface water samples collected in September. The sampling stations are shown in Figure 1 . The water samples were placed in polyethylene bottles (each bottle was acid-cleaned and rinsed with ultrapure water) and refrigerated at 4°C in a refrigerator. The water samples were taken back to the laboratory for further analysis on the same day. Furthermore, a spectroradiometer (USB2000+, Ocean Optics, Inc., USA) was used to measure the remote sensing reflectance spectra above the water’s surface between 200 and 1100 nm (1 nm interval) in accordance with the protocols proposed by Mobley (Mobley, 1999). Remote sensing reflectance was determined from an above-water method with an azimuth angle of 135° from the sun and 45° from the nadir (Mobley, 1999). At every water sampling site, the radiances from the sky, the water, and the reference panel were measured. Remote sensing reflectance (R_rs(λ)) was calculated using the following Equation 1:

where λ is the wavelength, L_u(λ) is the upwelling spectral radiance, L_s(λ) is the incident spectral sky radiance and r is the proportionality coefficient, with a value of 0.025 (Yu et al., 2023). L_p(λ) is the radiance from gray reference panel, ρ_p(λ) is the known reflectance of the gray panel (Yu et al., 2023).

Glass fiber filter membranes (pre-combustion at 450 °C for 4 h, GF/F, Whatman) with a 47-mm diameter were used to filter the TSM, Chl a, and PN samples. The weight method was used to calculate TSM concentrations (Zhou et al., 2021). Chl a in the GF/F filter was extracted using 90% acetone and analyzed using the fluorometric method (Lao et al., 2021; Zhou et al., 2021). An element analyzer, coupled with a stable isotope ratio mass spectrometer (EA Isolink-253 Plus, Thermo Fisher Scientific, Inc. USA) was used to measure the concentration of PN and δ¹⁵N_PN (Chen et al., 2021). The mean standard deviation of δ¹⁵N_PN and PN concentration was ±0.3‰ and ±0.3%, respectively (Chen et al., 2021).

The satellite data for this study was selected from the Ocean and Land Color Instrument (OLCI) data carried by Sentinel-3. The OLCI data contains a total of 21 spectral bands, ranging from 400 to 1020 nm, including 16 water-color bands, with a spatial resolution of 300 m and a global coverage time of 1-2 days (Su et al., 2021). It can achieve global multispectral medium resolution ocean/land observation capabilities. Sentinel-3 OLCI image data can be downloaded through the European Space Agency’s official website (ESA, https://scihub.copernicus.eu/dhus/#/home). We used the C2RCC (case 2 regional coast color) algorithm integrated into Sentinel Application Platform (SNAP) software to perform atmospheric correction on Sentinel-3 OLCI data. The image data after atmospheric correction was further processed and analyzed in SNAP software.

The accuracy evaluation of this study mainly includes four indices. Specific calculation formula of Pearson correlation coefficient (r), determination coefficient (R²), mean absolute percentage error (MAPE) and root mean square error (RMSE) are following Equations 2–5:

where n is the number of samples, X i and Y i refer to the value of two variables, Z and W denote the mean value of two variables in the sample.

As is shown in Figure 3 , the TSM concentration ranged from 5.25 to 45.35 mg/L (averaged at 13.70 mg/L) in May. It should be noted that operational errors prevented us from obtaining the Chl a concentration outside the bay in May. Chl a concentration inside the bay ranged from 0.68 to 19.33 μg/L (averaged at 5.49 μg/L) in May. While in September, the concentration of TSM and Chl a ranged from 2.60 to 62.40 mg/L (with an average of 11.63 mg/L) and from 1.62 to 21.88 μg/L (with an average of 7.53 μg/L), respectively ( Figure 3 ).

During the survey period, PN concentration ranged from 0.026 to 0.135 mg/L in May, with an average value of 0.048 mg/L. In September, PN concentration ranged from 0.022 to 0.09 mg/L, with an average value of 0.043 mg/L. Overall, the PN concentration in September was slightly lower than that in May. In addition, as shown in Figure 4 , the average concentration of PN outside the bay was higher than that inside the bay in May and September. In May, the average concentration of PN outside and inside the bay was 0.062mg/L and 0.044mg/L, respectively. In September, the average concentration of PN outside and inside the bay was 0.053mg/L and 0.034mg/L, respectively. The δ¹⁵N value ranged from 5.89‰ to 10.14‰ (average 7.77‰) in May and 3.73‰ to 12.08‰ (average 7.77‰) in September, respectively. Similar to the distribution of PN concentration, the average δ¹⁵N outside the bay was higher than that inside the bay in May and September. Figure 4 presents that the average value of δ¹⁵N outside and inside the bay in May was 8.71% and 7.43%, respectively. In September, the average value of δ¹⁵N outside and inside the bay was 10% and 5.9%, respectively. The spatial distribution differences of δ¹⁵N inside and outside the bay were more pronounced in September ( Figure 4 ).

Firstly, we conducted correlation analysis using in situ R_rs(λ) and δ¹⁵N_PN, and found that single band remote sensing reflectance retrieval was not effective, with low correlations (P>0.05), which will not be presented here. In order to obtain the best band combination of retrieval bands, we designed 6 optical indicators ( Table 3 ). The method of optical indicators refers to the R_rs(λ) combination forms designed by Ling et al (Ling et al., 2020). Concretely, 240 possible combinations of R_rs(λ) with the sixteen OLCI bands (400 nm, 413 nm, 443 nm, 490 nm, 510 nm, 560 nm, 620 nm, 665 nm, 674 nm, 681 nm, 709 nm, 754 nm, 761 nm, 764 nm, 768 nm, and 779 nm) were trained by using MATLAB R2018a software to determine those optimal results for each form, X. After training, based on the correlation coefficient (r) between X and δ¹⁵N_PN, the highest value was taken, and the selection of λ₁ and λ₂ was ultimately determined ( Table 3 ). From Table 3 , it can be seen that X2 and X5 perform well, with r values of -0.78 (P<0.01) and -0.77 (P<0.01), respectively. We selected these two optical indicators as input variables for the model to construct δ¹⁵N_PN retrieval model.

We randomly divided the dataset (input and output variables) into a training set (35 samples) and a testing set (17 samples), trained the model separately, and verified its performance. As shown in Figure 5 , from the evaluation indices R², MAPE, and RMSE ( Figure 5A-F ), BPNN, RF and MLR methods perform well on the training sets (BPNN: R^2 = 0.64, MAPE = 14.93%, RMSE = 1.32‰; RF: R^2 = 0.70, MAPE = 13.84%, RMSE = 1.20‰; MLR: R^2 = 0.65, MAPE = 15.14%, RMSE = 1.29‰) and testing sets (BPNN: R^2 = 0.84, MAPE = 10.71%, RMSE = 0.99‰; RF: R^2 = 0.65, MAPE = 13.66%, RMSE = 1.22‰; MLR: R^2 = 0.84, MAPE = 11.94%, RMSE = 1.03‰), which can meet our retrieval requirements for δ¹⁵N_PN. To further validate the performance of BPNN, RF and MLR methods for δ¹⁵N_PN retrieval, we obtained quasi-synchronous Sentinel-3 data (September 20, 2016) during the sampling period, and obtained 20 satellite-ground matching points that were less affected by clouds, shadows, and solar flares. The established BPNN, RF and MLR models were applied to Sentinel-3 data. As shown in Figure 6 , the BPNN and MLR model (BPNN: R^2 = 0.63, MAPE = 20.10%, RMSE = 1.63‰; MLR: R^2 = 0.63, MAPE = 20.71%, RMSE = 1.63‰) perform slightly better than the RF model (R^2 = 0.55, MAPE = 21.99%, RMSE = 1.67‰), with points more evenly distributed on both sides of the trend line. From Figure 6 , it can be seen that the overall accuracy of the BPNN model is comparable to that of the MLR model, but the MAPE of BPNN model is slightly lower than that of the MLR model. Therefore, we used the BPNN model as the model for retrieving δ¹⁵N_PN. It is worth noting that there is a certain degree of overestimation or underestimation of satellite retrieval value of δ¹⁵N_PN ( Figure 6 ). This may be due to the fact that the sampling time and satellite transit time are not synchronized in real-time (exceeding 24 hours), and the time window is an important factor affecting the accuracy of retrieval (Fu et al., 2023). On the other hand, it may be due to errors caused by atmospheric correction (Zhao et al., 2022). Moreover, Figure 7 shows the comparison between the OLCI-derived values of optical indicators (X2 and X5) and the measured values, demonstrating acceptable performance. The band combination reduces the errors caused by atmospheric correction in the algorithm implementation process to some extent (Zhao et al., 2022).

As shown in Figure 8 , the spatial distribution map obtained by applying the BPNN model to Sentinel-3 data shows that δ¹⁵N_PN outside Zhanjiang Bay is slightly higher than inside Zhanjiang Bay. However, a few areas affected by factories, docks, and aquaculture areas (circled in Figure 8 ) (Lu et al., 2020; Zhang et al., 2020a; Zhou et al., 2022), these areas are highly susceptible to the influence of sewage or wastewater, resulting in higher dynamic changes in δ¹⁵N_PN in these areas. Therefore, the δ¹⁵N_PN retrieval results in these regions may have some differences from the measured values. But overall, the retrieval results are relatively consistent with the measured results, which also indicates that the BPNN model has a certain reliability in retrieving δ¹⁵N_PN. Simultaneously utilizing Sentinel-3 data to retrieve δ¹⁵N_PN has great potential for application.

PN concentration and δ¹⁵N_PN in the ocean are influenced by various processes such as water mass mixing, nutrient gain and loss, and phytoplankton production (Sigman and Casciotti, 2001; Dagg et al., 2004; Ye et al., 2017). For bays strongly affected by human activities, PN concentration and δ¹⁵N_PN will also be affected by terrestrial factors such as soil, land runoff, and discharge of wastewater (Cloern et al., 2002; Bristow et al., 2013; Ye et al., 2017). Under the joint action of multiple factors, PN concentration and δ¹⁵N_PN in Zhanjiang Bay exhibited different characteristics in different months and regions.

TSM is the main carrier of terrestrial particulate matter (Chester and Jickells, 2012). As shown in Figure 9 , there was a significant positive correlation between PN concentration and TSM concentration in May (r=0.746, P<0.01), but there was no significant correlation between PN concentration and TSM concentration in September. This indicated that terrestrial particulate matter in May had an important impact on PN concentration, while the terrestrial component content of PN in September was relatively low. To a certain extent, the Chl a concentration reflects the status of phytoplankton production (Luhtala et al., 2013). As illustrated in Figure 9 , there was a certain positive correlation between PN concentration and Chl a concentration in May and September (r=0.647, P<0.01 for May; r=0.476, P<0.01 for September), indicating that phytoplankton production had a certain impact on PN concentration. In general, δ¹⁵N_PN can effectively indicate the source of PN (Ye et al., 2017; Chen et al., 2021). The δ¹⁵N_PN composition of marine organic matter range from 3‰ to 12‰ (Chen et al., 2021; Huang et al., 2021), while wastewater and livestock usually have the δ¹⁵N_PN values of 10‰ to 22‰ (Huang et al., 2021). The δ¹⁵N values ranged from 3.73‰ to 12.08‰ (average 7.77‰) in this study. Therefore, the source of PN in Zhanjiang Bay may be mainly marine organic matter, but terrestrial input was also mixed in, such as wastewater. In addition, there was a significant correlation between δ¹⁵N_PN and Chl a concentration in May and September (r=0.574, P<0.05 for May; r=0.806, P<0.01 for September; Figure 9 ). Both PN concentration and δ¹⁵N_PN showed a good correlation with Chl a concentration, indicating that phytoplankton production had a significant contribution to the source of PN.

Significantly, stations with higher δ¹⁵N_PN (>10 ‰) generally had higher Chl a concentration (>10 μg/L) ( Table 4 ). The reception of hypereutrophic municipal wastewater in the bay area can easily lead to the bloom of phytoplankton, resulting in higher concentrations of Chl a (Gao et al., 2021). When the growth rate of phytoplankton accelerates, the isotopic fractionation that occurs during the rapid absorption of inorganic nitrogen by phytoplankton can lead to a heavier nitrogen isotope composition of particulate organic matter (Mariotti et al., 1984). Due to the preferential absorption of NH₄ ⁺ during the growth process of phytoplankton, the strong nitrification in coastal water and the preferential utilization of ¹⁴N in NH₄ ⁺ by phytoplankton can lead to the accumulation of residual NH₄ ⁺ in water by ¹⁵N (Cifuentes et al., 1988). When phytoplankton continue to absorb these enriched ¹⁵N in NH₄ ⁺, it will cause an increase in the δ¹⁵N_PN value of the produced particulate organic matter (Cifuentes et al., 1988; Ke et al., 2017). Moreover, due to the rapid economic development and increased human activities in Zhanjiang, the process of heterotrophic bacteria has been intensified (Li et al., 2021). The strong biodegradation process prioritizes the degradation of organic matter containing lighter isotopes, leading to the enrichment of residual organic matter with heavy nitrogen isotopes (Li et al., 2021). It is worth noting that the δ¹⁵N_PN value at station S13 and A29 was relatively high (>10 ‰), but the Chl a concentration was not high (5.34 μg/L and 6.75 μg/L, respectively), this may be due to the impact of sewage or wastewater input, as S13 and A29 are located near the factory and aquaculture industry, respectively. Previous study showed that in the region of algal uptake, sewage-derived NH₄ ⁺ and sewage-derived NO₃ ^- could raise the δ¹⁵N_PN value by 9.0-17.2‰ and 10-15‰, respectively (Estep and Vigg, 1985; Leavitt et al., 2006). In addition, Zhou et al. (2021) also indicated that during non-typhoon periods in Zhanjiang Bay, the PN heavy isotopes can be attributed to the utilization of mineralized NH₄ ⁺ from wastewater by phytoplankton. Therefore, the relatively heavy δ¹⁵N_PN component in Zhanjiang Bay and its adjacent waters can be attributed to phytoplankton production and sewage or wastewater input.

In section 4.1, we obtained that the PN source in Zhanjiang Bay is mainly phytoplankton production, and δ¹⁵N_PN had a good correlation with concentration of Chl a. Therefore, the δ¹⁵N_PN can be linked to the water color parameter Chl a, which can establish a suitable δ¹⁵N_PN remote sensing retrieval model. It is well-known that there is an absorption peak in the spectral reflectance near the wavelength of 674 nm due to the absorption of phytoplankton pigments, while there is a fluorescence peak near the wavelength of 681 nm, both of which are spectral characteristic bands specific to Chl a (Su et al., 2021). Consequently, choosing these two bands to retrieve δ¹⁵N_PN has a certain scientificity and reliability. In the RF model, due to the dominance of measured values between 5 and 11 in the training dataset, there may be biases in the decision tree constructed by decision tree learners, resulting in predicted values being concentrated within the range of 5 to 11. In addition, during prediction, each tree is given a predicted value, and the average of all predicted values is taken. This results in the predicted value of the random forest being within the range of the training sample’s predicted values, so it cannot be extrapolated. The predicted value can only be between the minimum and maximum values of the training sample’s predicted values, resulting in a set of identical estimated values between 10 and 11 in the testing set. When we need to infer independent or non independent variables that are beyond the range, the random forest does not do well (Wang et al., 2023b). The solution is to expand the scope of the dataset in the future to maintain balance. The BPNN model and MLR model performed well on both the training and testing sets, and they also performed well in the application of Sentinel-3 data. Although the multiple linear regression algorithm is simple and ease of implement, it should be noted that the regression coefficients of this model is only suitable for the Zhanjiang Bay and its adjacent sea areas. For other sea areas, parameter regionalization may be required, and the applicability of the model needs further verification. The samples in this study were collected during the rainy season. During the rainy season, increased rainfall leads to an increase in nutrients carried into the sea by land runoff, resulting in an increase in phytoplankton biomass (Baek et al., 2009). There was a good correlation between Chl a, and δ¹⁵N_PN in May and September. This provides a good foundation for the establishment of δ¹⁵N_PN remote sensing model. However, during the dry season, the decrease in rainfall leads to changes in the physical, chemical, and biological conditions of the water, as well as changes in the activity of phytoplankton. Whether Chl a and δ¹⁵N_PN still maintain a good correlation remains to be further explored. Whether the δ¹⁵N_PN retrieval model we established is still applicable depends on further sample collection and verification.

In order to better evaluate the δ¹⁵N_PN remote sensing retrieval model established in this study, we selected six widely used Chl a retrieval algorithms ( Table 5 ), including 3 empirical algorithms (Three-band algorithm: TBA; Fluorescence Line Height algorithm:FLH; Maximum Chlorophyll Index algorithm: MCI) (Gower et al., 1999; Dall'Olmo et al., 2005; Gower et al., 2005) and 3 semi-analytical algorithms (Gons; Simis; Quasi-Analytical Algorithm improved form: QAA750E) (Gons et al., 2002; Simis et al., 2005; Xue et al., 2019), to attempt to retrieve δ¹⁵N_PN. As shown in Table 5 , we substituted the in situ remote sensing reflectance into the expressions of the following six algorithms, and then perform comparison analysis between the results of the expressions with measured δ¹⁵N_PN. It is worth mentioning that the comparison analysis between the absorption coefficient of phytoplankton (a_ph(λ)) obtained by the semi-analytical algorithm and δ¹⁵N_PN was conducted. From Table 5 , it can be seen that the BPNN algorithm proposed in this study has the highest accuracy (R^2 = 0.66, MAPE=13.52%, RMSE=1.19‰), followed by the TBA algorithm (R^2 = 0.49, MAPE=14.32%, RMSE=1.46‰). The Gons and Simis semi-analytical algorithms also have a good performance (Gons: R^2 = 0.46, MAPE=16.31%, RMSE=1.49‰; Simis: R^2 = 0.42, MAPE=16.81%, RMSE=1.56‰), indicating that semi-analytical algorithms have certain application potential in retrieving δ¹⁵N_PN. In addition, MCI algorithm, FLH algorithm and QAA750E algorithm perform poorly. Therefore, the algorithms composed of more bands may lead to more indeterminacy factors introduced, which directly affects the retrieval accuracy.

Additionally, we obtained 18 measured δ¹⁵N_PN values in September 2017 and applied the BPNN model established in this study to the Sentinel-3 OLCI image data (18 September 2017). The retrieval results were shown in Figure 10 . From the perspective of retrieval accuracy ( Figure 11 ), the R², RMSE and MAPE were 0.59, 1.78‰, and 34.06%, respectively. The retrieval results can meet the requirements to a certain extent, indicating that our δ¹⁵N_PN retrieval model has certain applicability.

Isotope fractionation gives PN from different sources specific nitrogen stable isotope characteristic values, which provides the possibility of determining the source and destination of PN, making δ¹⁵N_PN a valuable tracer for tracking PN sources and understanding N cycling in water systems (Chen et al., 2021; Huang et al., 2021; Lu et al., 2021). Although traditional field surveys and chemical methods can accurately obtain δ¹⁵N_PN, they consume a lot of time, manpower, and resources. How to improve efficiency and enable us to quickly and extensively understand the dynamic changes of δ¹⁵N_PN? Satellite remote sensing has developed rapidly in recent decades, and various high-performance sensors have been developed for marine environmental monitoring, which is very conducive to our research and exploration of the ocean. Satellite remote sensing has irreplaceable advantages in large-scale spatial and long-term series monitoring. We only need to sacrifice a small amount of accuracy to obtain acceptable results, which is of great significance for the biogeochemical processes and nitrogen cycling research of PN in the ocean. According to the analysis above, water environmental pollution can also be distinguished by δ¹⁵N_PN values, which can indirectly indicate the pollution level of water bodies. This also provides a new method and strategy for traditional water quality monitoring and management.

However, the δ¹⁵N_PN remote sensing model established in this study only involves data from two cruises, and it cannot be denied that the limitations of the model exist. Meanwhile, the quality of satellite images and atmospheric correction can also bring uncertainty to δ¹⁵N_PN estimation. In the future, we will increase the sampling frequency and use more data from different seasons and regions to validate our established model, enhancing its robustness and universality.