A nitrate budget of the Bohai Sea based on an isotope mass balance model

The Bohai Sea (BHS) is a semi-closed marginal sea impacted by one of the most populated areas of China. The supply of nutrients, markedly that of reactive nitrogen, via fluvial and atmospheric transport has strongly increased in parallel with the growing population. Therefore, it is crucial to quantify the reactive nitrogen input to the BHS and understand the processes and determine the quantities of nitrogen eliminated in and exported from the BHS. The nitrogen budget and in particular the internal sources and sinks of nitrate were constrained by using a mass-based and dual stable isotope approach based on δ¹⁵N and δ¹⁸O of nitrate (NO${}_{\mathrm{3}}^{-}$). Samples of water, suspended matter, and sediments were taken in the BHS in spring (March and April) and summer (July and August) 2018. The Yellow River (YR) was sampled in May and July to November, and Daliao River, Hai River, Luan River, and Xiaoqing River were sampled in November of 2018. In addition to nutrient, particulate organic carbon, and nitrogen concentrations, the dual isotopes of nitrate (δ¹⁵N and δ¹⁸O), δ¹⁵N of suspended matters, and sediments were determined. Based on the available mass fluxes and isotope data an updated nitrogen budget is proposed. Compared to previous estimates, it is more complete and includes the impact of interior cycling (nitrification) on the nitrate pool. The main nitrate sources are rivers contributing 19.2 %–25.6 % and the combined terrestrial runoff (including submarine fresh groundwater discharge of nitrate) accounting for 27.8 %–37.1 % of the nitrate input to the BHS, while atmospheric input contributes 6.9 %–22.2 % to total nitrate. An unusually active interior nitrogen cycling contributes 40.7 %–65.3 % to total nitrate via nitrification. Nitrogen is mainly trapped in the BHS and mainly removed by sedimentation (70.4 %–77.8 %), and only very little is exported to the Yellow Sea (YS) (only 1.8 %–2.4 %). At present denitrification is active in the sediments and removes 20.4 %–27.2 % of nitrate from the pool. However, a further eutrophication of the BHS could induce water column hypoxia and denitrification, as is increasingly observed in other marginal seas and seasonally off river mouths.

1 Introduction

Reactive nitrogen (N_r) is an essential nutrient of life on earth, but only a few organisms can fix nitrogen directly from the large atmospheric dinitrogen pool. Since the invention of the Haber–Bosch process the amount of fixed nitrogen (120 Tg N yr⁻¹) has constantly grown, and since 2010 it has exceeded the natural terrestrial sources of reactive N of 63 Tg N yr⁻¹ (Fowler et al., 2013). Hotspots of agricultural N fertilizer application shifted from the US and western Europe in the 1960s to eastern Asia in the early 21st century (Lu and Tian, 2017). In China, the Haber–Bosch process produces 37.1 Tg N yr⁻¹, which is almost 3 times that of the biological N fixation of 12.0 Tg N yr⁻¹. An estimated 32.0 Tg N yr⁻¹ is produced for fertilizer (Gu et al., 2015), and China accounted for 29 % of the global ammonium production in 2018 (IFA, 2019). The leakage and volatilization of this man-made reactive nitrogen have strongly impacted limnic and marine ecosystems in China. Riverine reactive nitrogen discharged to the ocean from China was estimated at 5.4 Tg N yr⁻¹ (Gu et al., 2015). The total load of Chinese major estuaries to coastal seas was about 9 % of the global river load for dissolved inorganic nitrogen (DIN) and 1.5 % of the global phosphate load (Smith et al., 2003; Liu et al., 2009).

The Bohai Sea (BHS) is a semi-enclosed basin with a surface area of 77 × 10³ km² and an average depth of 18 m (Chen, 2009; Su, 2001) that is heavily impacted by human activities in one of the most densely populated terrestrial catchments of the world. It exchanges salt water with the Yellow Sea (YS) through Bohai Strait, and the Yellow River (YR) is a major source of freshwater to BHS (Chen, 2009). During last 50 years, rising anthropogenic activity in the catchment induced severe environmental changes in the BHS, including increasing salinity, temperature, concentrations of DIN, and changes in stoichiometric nutrient ratios (Zhao et al., 2002; Zhang et al., 2004; Wang et al., 2019; Ning et al., 2010; Xin et al., 2019). DIN concentrations increased from 0.30 to 3.55 µmol L⁻¹ in the time from 1982 to 2009, while phosphate (from 0.76 to 0.31 µmol L⁻¹) and silicate (26.6 to 6.60 µmol L⁻¹) concentrations significantly decreased so that $\mathrm{N}/\mathrm{P}$ increased about 30 times (Zhang et al., 2004; Liu et al., 2011). Phytoplankton nutrient limitation in the BHS switched from nitrogen to phosphorus in the period of the 1980s to the 1990s, and this limitation pattern persists until the present day (Xu et al., 2010; Liu et al., 2009; Wang et al., 2019).

The total annual water discharge of rivers into BHS is about 68.5 × 10⁹ m³ yr⁻¹, of which the YR accounts for more than 75 % (Liu et al., 2011). The water exchange time of the YR estuary is only 0.1–0.2 d (Liu et al., 2009), which implies a fast transfer of nutrients into the open BHS, and many of these are trapped in Laizhou Bay (Fig. 1) (Zhang et al., 2004). The atmospheric deposition of nitrate (3.42 × 10⁹ mol yr⁻¹) in BHS was modeled to be less than riverine nitrate (7.25 × 10⁹ mol yr⁻¹), while more ammonium was supplied from atmospheric deposition (6.15 × 10⁹ mol yr⁻¹) than from riverine input (0.93 × 10⁹ mol yr⁻¹) in the 1990s (Zhang et al., 2004). BHS nitrate budgets reported during the last two decades were not completely constrained because crucial data, such as groundwater discharge or nitrification, were not available (Zhang et al., 2004; Liu et al., 2003, 2009, 2011). There are few published nutrient data from the BHS over the last decade, and the terms in the N_r budget of BHS concerning the quantities of N_r generated or eliminated by biogeochemical cycling within the basin have not been addressed.

Figure 1 Sampling sites in Bohai Sea and rivers (Xiaoqing River not shown). Open red circles stand for sampling sites in spring, and open blue triangles stand for sampling sites in summer; names of the sites are marked nearby. Black arrows stand for the most significant current flows in and out of the BHS. Dashed blue lines stand for two main sections, and the dashed black line stands for the boundary of our study area. (Base map sourced by Ocean Data View, https://odv.awi.de, Schlitzer, 2021)

The nitrogen budget and in particular the internal sources and sinks of nitrate can be constrained with a mass-based and dual stable isotope approach based on δ¹⁵N and δ¹⁸O of nitrate. The combination of data permits tracking of nitrate and quantification of internal cycling of inorganic nitrogen (Sigman et al., 2005; Wankel et al., 2006; Sugimoto et al., 2009; DiFiore et al., 2006; Montoya et al., 2002; Emeis et al., 2010). Stable isotopes of reactive nitrogen have been used to explore nitrogen sources in the eastern Chinese seas (Umezawa et al., 2014; Wang et al., 2016; S. M. Liu et al., 2017; Wu et al., 2019; Liu et al., 2020; Yu et al., 2021; Y. Li et al., 2019) and the Changjiang Estuary (Yu et al., 2015; Wang et al., 2017; Yang et al., 2018; Chen et al., 2013).

For this study, we analyzed water, suspended matter, and sediments in the Bohai Sea sampled during the spring and summer seasons for nutrient concentrations, carbon and nitrogen contents, dual isotopes (δ¹⁵N and δ¹⁸O) of nitrate, and δ¹⁵N of particulate nitrogen. The aim of the study is to characterize and quantify N_r sources and sinks, in particular those from internal cycling processes that have not been included in previous budgets (Zhang et al., 2004; Liu et al., 2003, 2009, 2011), to track the fate of N_r in the present Bohai Sea. The observation data presented here are the basis for a combined mass and isotope balance model, results of which will be a basis for future studies on the rising impact of fast-growing mega cities in the BHS catchment and their possible impact on the adjacent YS.

2 Materials and methods

2.1 Sample collection

Research cruises were carried out on board R/V Dongfanghong 2 in spring and summer 2018 with 24 sampling sites in April and 25 sites in August, respectively (Fig. 1). Water samples were taken from several depths by 12 L Niskin bottles attached to a CTD rosette (911plus, Seabird, USA). The water samples for nutrients and isotope analysis were filtered using nucleopore polycarbonate filters (0.4 µm) with plastic Nalgene filtration units. The filtered water was collected in Falcon PE tubes (45 mL) as duplicates for measurements of nutrients and nitrate isotope analysis, frozen immediately (−20^∘C), and kept frozen until analyses in the home laboratory in Germany. Between 1 to 8 L of water was filtered through pre-weighted GF/F filters (0.7 µm, Φ=47 mm, Sigma Aldrich) which had been pre-combusted at 450^∘C for 4 h. The filters were subsequently dried on board under 45^∘C for 24 h. Surface sediments were taken with a box corer, and surface samples were transferred into plastic bags with a metal spoon, frozen at −20^∘C, and kept frozen until later analysis in the home lab.

YR water samples were taken from the Kaiyuan floating bridge in Lijin, located 44 km upstream of the river mouth. Samples were taken in the middle of the river course with a plastic reversing water sampler at 1 m under the surface. The water samples were filtered immediately for nutrient analysis and collection of suspended particles, and subsequently they were stored frozen until delivered to the home laboratory. Samples were taken monthly in May and July to November of 2018 from the YR and in November of the same year from Daliao River, Hai River, Luan River, and Xiaoqing River (Fig. 1).

2.2 Measurements of nutrients and nitrate isotopes

Nutrient concentrations were measured with an AutoAnalyzer 3 system (Seal Analytics) using standard colorimetric methods (Grasshoff et al., 2009). The relative error of duplicate sample measurements was below 1.5 % for NO_x and phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$) concentrations, as well as below 0.3 % for ammonium. The detection limit was <0.05 µmol L⁻¹ for NO_x, >0.01 µmol L⁻¹ for PO${}_{\mathrm{4}}^{\mathrm{3}-}$, and >0.013 µmol L⁻¹ for ammonium.

Both δ¹⁵N and δ¹⁸O of nitrate (δ¹⁵N = [(¹⁵N $/$ ¹⁴N)_sample $/$ (¹⁵N $/$ ¹⁴N)_standard−1] × 1000 ‰, δ¹⁸O = [(¹⁸O $/$ ¹⁶O)_sample $/$ (¹⁸O $/$ ¹⁶O)_standard−1] × 1000 ‰) were determined following the denitrifier method (Sigman et al., 2001; Casciotti et al., 2002). Only the samples with nitrate concentrations >1.7 µmol L⁻¹ were analyzed, and δ¹⁵N and δ¹⁸O were analyzed in one sample run. Water samples were injected into a suspension of the denitrifier Pseudomonas aureofaciens with injection volumes adjusted to yield 10 nmol N₂O. The N₂O gas was purged by helium into a GasBench II (Thermo Finnigan) for purification. Afterwards the N₂O gas was analyzed by a Delta V Advantage and a Delta V Plus mass spectrometer. Samples were measured in duplicate, and the two international standards IAEA-NO${}_{\mathrm{3}}^{-}$ (δ¹⁵N–NO${}_{\mathrm{3}}^{-}=+\mathrm{4.7}$ ‰, δ¹⁸O–NO${}_{\mathrm{3}}^{-}=+\mathrm{25.6}$ ‰) and USGS-34 (δ¹⁵N–NO${}_{\mathrm{3}}^{-}=-\mathrm{1.8}$ ‰, δ¹⁸O–NO${}_{\mathrm{3}}^{-}=-\mathrm{27.9}$ ‰) and an internal potassium nitrate standard were measured in each batch. The data were corrected by applying a bracketing correction (Sigman et al., 2009), and the standard deviations of the international and in-house standards were found to be ≤ 0.2 ‰ for δ¹⁵N and ≤ 0.5 ‰ for δ¹⁸O. The standard deviations of duplicate samples were in the same range. Nitrite affects the results and was removed following the protocol of Granger and Sigman (2009) whenever [NO${}_{\mathrm{2}}^{-}$] exceeded 5 % of the NO_x pool. In all other cases, we report combined (nitrite + nitrate) values.

2.3 Measurements of suspended matters and sediments

The tared GF/F filters were weighed to calculate the amount of suspended particulate matter (SPM) per liter of water. Total carbon and nitrogen concentrations in SPM and sediment samples were measured by a Euro EA 3000 (Euro Vector SPA) elemental analyzer, and SPM samples with high carbon and nitrogen contents and sediments were acidified to measure organic carbon content. The precision of total and organic carbon determination is 0.05 %, that of nitrogen is 0.005 %, and the standard deviations are less than 0.08 for total and organic carbon and 0.02 for nitrogen. Nitrogen isotope ratios were determined with a FlashEA 1112 coupled to a MAT 252 (Thermo Fisher Scientific) isotope ratio mass spectrometer. The precision of nitrogen isotope analyses is better than 0.2 ‰, and the standard deviation is less than 0.03 ‰.

2.4 Measurements of dissolved oxygen

The dissolved oxygen (DO) samples were collected, fixed, and titrated on board following the Winkler procedure at an uncertainty level of <0.5 %. A small quantity of NaN₃ was added during subsample fixation to remove possible interferences from nitrite (Wong, 2012). The DO saturation (DO %) was calculated from field-measured DO concentration divided by the DO concentration at equilibrium with the atmosphere which was calculated from temperature, salinity, and local air pressure, as per the Benson and Krause (1984) equation.

2.5 Hydrodynamic model of nutrient export from BHS to the Yellow Sea

The regional three-dimensional hydrodynamic Hamburg Shelf Ocean Model, HAMSOM (Backhaus, 1985), was applied in the East China Sea (23–45^∘ N, 117–131^∘ E) to calculate the water and nutrient transport through the Bohai Strait for the year 2018. The HAMSOM model has been applied to investigate the Bohai Sea physical circulation for several decades now and has been extensively validated in the Bohai Sea (Jia and Chen, 2021; Hainbucher et al., 2004; Huang et al., 1999). The spatial resolution of the model is 2^′ (approx. 3.7 km) with 20 layers in the vertical direction, while the calculation time step is 3 min. The upper 50 m of the HAMSOM model are resolved by layers of 5 m thickness. The topography data (resolution of 2^′) were obtained from marine navigation charts. The meteorological forcing was derived from an hourly ERA5 dataset with a spatial resolution of 0.25^∘ (Copernicus Climate Change Service, 2017). The open-boundary SSH (sea surface height) and the boundary T and S data and for the initial T and S fields were extracted from the daily Mercator Ocean dataset ($\mathrm{1}/\mathrm{12}$^∘) (Lellouche et al., 2019). A total of 13 partial tides derived from the TPXO8-atlas v1 were superimposed onto the SSH along the open boundary (Egbert and Erofeeva, 2002). The observed monthly river discharges were available for the two largest rivers, i.e., Changjiang and YR (MWR, China, 2015–2018), while the inputs for the remaining rivers were derived from the WaterGAP dataset (0.5^∘, monthly climatology) (Müller Schmied et al., 2014). The spin-up period of this model is 1 year.

Four sites on a north–south section through the Bohai Strait, i.e., B33, B34, B35, and B36, have been selected to represent the open boundary of the BHS (Fig. 1). The simulated SSH and current velocities (west–east component) were extracted along this section. In addition, nutrient concentrations were interpolated from the observed data at the four sites to the grid of the hydrodynamic model along the Bohai Strait section. Since the observational data just include spring and summer values, the mean value of nitrate in spring and summer had to be extrapolated to an entire year.

3 Results

3.1 Hydrological properties and nutrients

3.1.1 Hydrological properties

Averages of salinity and temperature in spring were 32.3 ± 0.5 (n=72) and 4.7 ± 0.8 ^∘C (n=72), respectively, and the water column was vertically mixed (Figs. 2 and 3). The YR discharged relatively warm water, and the lowest salinity was observed in the southeast of the YR estuary (site B68, T>6^∘C, S<31). Thus, the Yellow River diluted water (YRDW) is here defined as the water off the YR estuary with salinities lower than 31. In summer, averages of salinity and temperature were 31.6 ± 0.8 (n=88) and 22.4 ± 4.2 ^∘C (n=88), respectively, and the surface layer was stratified. The YRDW extended to an even larger area than in spring caused by high river discharge. The YRDW turned northeast towards Liaodong Bay into the central BHS in the surface layer (T>27^∘C, S<31).

Figure 2 Temperature (^∘) and salinity of section 1 for spring (a, c) and summer (b, d).

Figure 3 Temperature (^∘) and salinity of section 2 for spring (a, c) and summer (b, d).

The water column oxygen concentrations (see Supplement S1) in the study area in spring and summer were 10.27–11.47 and 3.84–8.86 mg L⁻¹, respectively, and thus much higher than the threshold for water column denitrification (0.15 mg L⁻¹) (Altabet, 2006). The detailed results of DO and other parameters are shown in Supplement S2.

3.1.2 Nutrient concentrations and distributions

Nutrient concentrations in spring were almost vertically uniform, consistent with temperatures and salinities, and no distinct nutricline was observed (Figs. 4 and 5). The prominent plume of low NO${}_{\mathrm{3}}^{-}$ of the middle layer along stations B45 and B63 may mark the influence of water advected from more offshore sites with low NO${}_{\mathrm{3}}^{-}$ concentration. The YR discharge is indicated by a nitrate maximum at station B68 (31.0 µmol L⁻¹) and the lowest observed salinity (Fig. 5). Concentrations of NH${}_{\mathrm{4}}^{+}$ were higher than 1 µmol L⁻¹, PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations were lower than 0.4 µmol L⁻¹ adjacent to the YR estuary, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ increased towards the central BHS (along section 1) and north of Shandong Peninsula (along section 2). The average concentrations of NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ of all spring samples were 6.5 ± 5.8, 0.2 ± 0.2, 0.8 ± 0.5, and 0.4 ± 0.2 µmol L⁻¹ (n=72), respectively.

Figure 4 Nutrients of section 1 (µmol L⁻¹) for spring (a, c, e) and summer (b, d, f).

Figure 5 Nutrients of section 2 (µmol L⁻¹) for spring (a, c, e) and summer (b, d, f).

In summer, in the central BHS (section 1), the upper 10 m water layer was nutrient depleted, while concentrations were high in the lower layer (Fig. 4). NO${}_{\mathrm{3}}^{-}$was almost depleted in the BHS except for stations B66 and B67 (included in section 2) located in the tongue of YRDW. NO${}_{\mathrm{2}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ concentrations were higher than 2 and 5.5 µmol L⁻¹, respectively, in the same area (Fig. 5). The average concentrations of NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ of 1.9 ± 2.7, 0.8 ± 1.1, 1.6 ± 1.9, and 0.1 ± 0.1 µmol L⁻¹ (n=85), respectively.

3.2 Dual isotopes of nitrate

In spring, the average values of δ¹⁵N and δ¹⁸O were 7.8 ± 1.4 ‰ and 12.0 ± 3.3 ‰ (n=52), and δ¹⁵N and δ¹⁸O ranged between 5.6 ‰ and 10.5 ‰ and 5.0 ‰ and 20.0 ‰, respectively. δ¹⁵N and δ¹⁸O were both vertically homogenous, corresponding to the vertical profile of salinity, temperature, and nitrate. Relatively high δ¹⁵N values were observed in the southwestern part of BHS including the areas adjacent to Bohai Bay and Laizhou Bay, which is consistent with high nitrate concentrations. The δ¹⁸O was inversely related to δ¹⁵N with low values near the YR estuary and high values in the northeast of the central BHS and north of Shandong Peninsula (Figs. 6 and 7).

Figure 6 Nitrate isotopes (‰) of section 1 for spring (a, c) and summer (b, d).

Figure 7 Nitrate isotopes (‰) of section 2 for spring (a, c) and summer (b, d).

Due to the low nitrate concentration in summer, only a subset of samples could be analyzed, and most of these are from the YRDW that had [NO${}_{\mathrm{3}}^{-}$] > 1.7 µmol L⁻¹. The average values of δ¹⁵N and δ¹⁸O were 9.9 ± 3.5 ‰ (n = 23) and 8.7 ± 3.3 ‰ (n=23), and they ranged between 3.5 ‰ and 23.9 ‰ and 3.1 ‰ and 18.4 ‰, respectively. The mean value of δ¹⁵N was higher than that of spring samples, whereas the δ¹⁸O value was lower. Relatively high δ¹⁵N and δ¹⁸O values of 23.9 ‰ and 18.4 ‰ were registered in the surface water of the YRDW (site B62) and decreased with water depths at increasing nitrate concentrations.

3.3 Suspended particulate matter

In spring suspended particulate matter (SPM) concentrations were mostly vertically homogenous along both transects with high values (>15 mg L⁻¹) close to the YR mouth (see Supplement S3). C_org % and N % are inversely correlated with SPM concentrations, and high values occurred in the central BHS and north of Shandong Peninsula. In summer SPM concentrations were significantly higher than in spring, and maxima occurred in deep water off the YR (>30 mg L⁻¹) and in the west part of the BHS. C_org % and N % maxima occur in surface waters in the eastern Laizhou Bay and the central BHS.

The average δ¹⁵N of SPM in spring was 4.8 ± 0.9 ‰ (n=14), with a maximum (>6.10 ‰, n=3) south of Luan River. The lowest values were observed in the southern Bohai Strait and northeast of the YR estuary. The other samples varied in a narrow range of 3.9 ‰–4.7 ‰ (n=11). In summer, the average of δ¹⁵N was 5.7 ± 0.8 ‰ (n=34) and ranged from 3.9 ‰ to 7.2 ‰. Systematic variation in δ¹⁵N of SPM was barely discernable and only exhibited a weak decline from the YR mouth into the northeastern BHS (section 1) and into the Bohai Strait (section 2) (see Supplement S3), thus tracking the salinity dilution gradient in the surface layer.

3.4 The discharge of the Yellow River

The water discharge of the year 2018 determined at the Lijin hydrography station was 333.8 × 10⁹ m³, which was by 14 % higher than the multi-year average of 292.8 × 10⁹ m³ (1952–2015) (MWR, China, 2019). The monthly mean discharge was 27.80 ± 20.21 × 10⁹ m³ per month (n=12), which was higher than the multi-year average value by 14 %–51 %, indicating that in the YR basin 2018 was a flood year. The water discharge maximum was from July to October (Fig. 8). During May to November in 2018 (no value for June), the mass fluxes of nitrate and phosphate increased with the water discharge of the YR, while the monthly fluxes of nitrite and ammonium exhibited an opposite trend. δ¹⁵N and δ¹⁸O of YR nitrate ranged from 9.1 ‰ to 10.9 ‰ and 1.1 ‰ to 3.0 ‰, respectively. The monthly mass-weighted average values of δ¹⁵N and δ¹⁸O were 9.9 ‰ and 1.9 ‰, respectively. δ¹⁵N and δ¹⁸O were positively correlated (δ¹⁸O = 1.04 × δ¹⁵N − 8.43, R²=0.85), as were δ¹⁵N and the monthly mass flux of nitrate (R²=0.63).

Figure 8 Monthly variation in water flux (shadowed bar), nitrate concentration (dashed line), mass flux of nitrate (solid line), and dual nitrate isotopes (δ¹⁵N in upper solid red line and δ¹⁸O in lower solid blue line) of the YR in 2018. Water samples were collected monthly in May and July to November in Lijin, nitrate concentrations (µmol L⁻¹) and dual isotopes (δ¹⁵N and δ¹⁸O, in per mil) were measured in home laboratory, and monthly water flux data (10⁸ m³) are according to the YR Sediment Bulletin 2018 (http://www.yrcc.gov.cn/nishagonggao/2018/index.html#p=20, last access: 28 April 2022).

3.5 Nitrate exchange with the Yellow Sea

Based on current velocities and nutrient concentrations along the section crossing the Bohai Strait, the annual water and nitrate export from BHS to the YS in the year 2018 was calculated to 1.26 × 10⁻³ Sv (1 Sv = 10⁶ m³ s⁻¹) and 0.9 × 10⁹ mol yr⁻¹, respectively. In this study, the exported nitrate is assumed to have the average isotope values of the BHS (δ¹⁵N = 8.9 ‰ and δ¹⁸O = 10.4 ‰).

Making use of the three-dimensional model (HAMSOM) results, it is also possible to determine a spatial distribution of the annual nutrient flux through the Bohai Strait section (Fig. 9). Positive values represent a nutrient flux out of the Bohai Sea, while negative ones indicate a flux into the Bohai Sea. The strongest nutrient export occurs at the southern part of the Bohai Strait, while the major import takes place in the upper 15 m in the northern Bohai Strait (Fig. 9).

Figure 9 Simulated net fluxes of nitrate (mol yr⁻¹) in the Bohai Strait for the year 2018. White dots are the sampling sites, the color bar is the flux of nitrate (10⁷ mol yr⁻¹), and positive values stand for the export of nitrate from BHS to YS.

4 Discussion

4.1 The hydrographic and nutrient characteristics in spring and summer

The sampling in early spring occurred during a season of low biological activity with the result that nutrients behaved almost conservatively (Wei et al., 2003; Tan and Shi, 2006). The YR discharged 333.8 × 10⁹ m³ water and 8.0 × 10⁹ mol nitrate to the Bohai Sea in 2018, accounting for 85 % and 84 % of water and nitrate discharge of all large rivers in the Bohai Sea, respectively. As is indicated by the negative correlations of NO${}_{\mathrm{3}}^{-}$ ($r=-\mathrm{0.78}$, p<0.01) and NH${}_{\mathrm{4}}^{+}$ ($r=-\mathrm{0.79}$, p<0.01) with salinity, the YR is one of the major sources of these nutrients in the BHS, whereas PO${}_{\mathrm{4}}^{\mathrm{3}-}$ is contributed by the inflow of saline waters from the YS as indicated by the positive correlation with salinity (r=0.43, p<0.01). Concentrations of nitrate were relatively high in the southern Bohai Strait but low in its northern part, suggesting that in spring nitrate-rich water flows out of Bohai Strait along the northern shore of Shandong Peninsula in the LCC (Lubei Coastal Current), while nitrate-depleted water flows in from the northern YS via the northern strait. In summer, the water is stratified with the thermocline at about 8 m water depth and coinciding with halo- and nutriclines. Nutrients are depleted to trace amounts above the thermocline. In contrast to the other nutrients, phosphate concentrations did not increase with depth in the southwestern part of BHS (i.e. Bohai Bay and Laizhou Bay). Similar to the spring situation, salinity was weakly positively correlated with PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (r=0.29, p<0.05) and NO${}_{\mathrm{2}}^{-}$ (r=0.32, p<0.05), whereas it was negatively correlated with NO${}_{\mathrm{3}}^{-}$ ($r=-\mathrm{0.69}$, p<0.01) and NH${}_{\mathrm{4}}^{+}$ ($r=-\mathrm{0.37}$, p<0.01), respectively. The average $\mathrm{N}/\mathrm{P}$ ratio in BHS in spring and summer was 28.2 ± 38.2 and 86.9 ± 126.3, respectively, implying that productivity in BHS was phosphorus limited. Thus, diazotrophic N₂ fixation is excluded as a significant input of N_r because of high $\mathrm{N}/\mathrm{P}$ ratios.

4.2 The main sources and sinks of nitrate in the BHS

Most of the external and internal sources of nitrate to the BHS are characterized by distinct dual isotope values. These fingerprints combined with mass flux estimates are in the following used to constrain the mass and isotope budget of nitrate in the BHS. Specifically, the role of internal cycling processes can thus be quantified, which was lacking in previous budgets. In the following, each of the sources and sinks is described, along with isotope composition or isotope fractionation associated with cycling processes.

4.2.1 Riverine inputs

The main input from this source of BHS N_r is from the YR, and we calculated a nitrate input of the YR of 7.95 × 10⁹ mol yr⁻¹, based on the annual average discharge (Yellow River Conservancy Commission of MWR, China, 2019) and the load-weighted nitrate concentration during our sampling period during the year 2018. The water discharge of the other seven important rivers (Hai River, Liao River, Daliao River, Luan River, Xiaoqing River, Daling River, and Xiaoling River) sums up to 59.86 × 10⁸ m³ yr⁻¹ (MWR, China, 2019; Ma et al., 2004; Yu et al., 2018; Zhang et al., 2004). Nitrate fluxes of these rivers were calculated based on water discharge and our nitrate measurements in the YR, Hai River, Daliao River, Luan River, and Xiaoqing River. Owing to a lack of data on nitrate concentrations of Liao River, Daling River, and Xiaoling River and considering that these rivers drain basins adjacent to Daliao River, we assume that nitrate concentrations of these rivers were the same as those of the Daliao River. The total riverine input of nitrate summed up to 9.49 × 10⁹ mol.

The mass-weighted average annual values for δ¹⁵N and δ¹⁸O of nitrate in these rivers were 10.0 ‰ and 1.3 ‰, respectively, taken here to represent the river nitrate isotopic composition discharged into BHS.

4.2.2 Submarine groundwater input

The DIN supplied to BHS by submarine groundwater discharge (SGD) flux has been estimated to be 2–10 times the YR discharge (Luo and Jiao, 2016; Peterson et al., 2008; Wang et al., 2015). These fluxes of SGD are a mixture of submarine fresh groundwater discharge (SFGD) and recirculated saline groundwater discharge (RSGD) (Liu et al., 2011; Peterson et al., 2008; J. Liu et al., 2017), but only the freshwater component is relevant as a source for the budget.

The latest estimates of SGD for the BHS are 10.0 × 10⁹ m³ yr⁻¹ (J. Liu et al., 2017) and 19.1 × 10⁹ m³ yr⁻¹ (Wang et al., 2015). The nitrate concentration of SFGD is not documented, but the nitrate concentration of groundwater in the YR delta was 304.2 ± 254.2 µmol L⁻¹ (Liu et al., 2011). For reducing the error of these indirectly measured data, we decided to use the value of 4.25 × 10⁹ mol yr⁻¹ for the nitrate flux of SFGD, which is the averaged products of SFGD water fluxes and nitrate concentrations shown above. This value is only approximately 10 % of previous estimates of the input nitrate for SGD into BHS (Liu et al., 2011) due to the exclusion of RSGD.

Because of pollution and denitrification processes in soils, aquifers, and groundwater (Zhang et al., 2013; Chen et al., 2007; Soares, 2000), the values of δ¹⁵N and δ¹⁸O of nitrate in SFGD are more enriched than those of river runoff, and this is illustrated by the observed δ¹⁵N value of 20.2 ± 9.0 ‰ (n=19) of on-land groundwater near the YR delta (Chen et al., 2007). As there are as yet no reported δ¹⁵N values of SFGD and RSGD inputs, we decided to take this value as the signal of nitrate δ¹⁵N imported by SFGD into the BHS. There are no data available for δ¹⁸O of nitrate of SFGD, and we will discuss possible constraints in the box model discussion (Sect. 4.3.2).

4.2.3 Atmospheric deposition

Combined atmospheric input by wet and dry deposition ranged from 3.14 × 10⁹ to 3.42 × 10⁹ mol yr⁻¹ (Zhang et al., 2004; Liu et al., 2003). We adopted the annual mass of NO_x deposition for China of 6.2 Tg yr⁻¹ (Zhao et al., 2017) and related this value to the area of the BHS, which results in an annual deposition of 3.6 × 10⁹ mol yr⁻¹. Owing to a lack of directly measured data for atmospheric NO_x in the BHS, we adopt 3.42 ± 2.29 × 10⁹ mol yr⁻¹ (Zhang et al., 2004) as the atmospheric nitrate flux.

The nitrate δ¹⁵N values of PM_2.5 (fine particulate matter suspended in the air) ranged from 3.5 ‰ to 17.8 ‰ in northern China (Fan et al., 2019; Zong et al., 2017; Song et al., 2020; Zhang et al., 2019), whereas nitrate δ¹⁵N of precipitation ranged from −2.5 ‰ to +0.9 ‰ (Zhang et al., 2008; Chang et al., 2019; Z. Li et al., 2019; Kim et al., 2019). Assuming that PM_2.5 is the main component of dry deposition we use a δ¹⁵N value of 8.2 ‰ reported from Beihuangcheng island in Bohai Strait (Zong et al., 2017), while the wet deposition has a δ¹⁵N value of −2.35 ‰ (Chang et al., 2019). Wet deposition in BHS was estimated as 54 %–68 % of total deposition (Liu et al., 2003; Zhang et al., 2008; Zhao et al., 2017), resulting in a mass-weighted average of BHS atmospheric deposition ranging from 1.03 ‰ to 2.56 ‰, of which the arithmetic mean value of 1.80 ‰ is adopted here.

The nitrate δ¹⁸O of PM_2.5 in BHS ranged from 65.0 to 88.1 ‰ (seasonally) (Zong et al., 2017), the value in Beijing is 88.3 ± 6.9 ‰ (Song et al., 2020), and it is 57.80 ± 4.23 ‰ in cloud samples of Shandong (Chang et al., 2019). For nitrate δ¹⁸O of dry deposition, we thus assume a mean value of 80.5 ‰ and that wet deposition has a δ¹⁸O of 57.8 ‰. These two estimates combined give a ratio of dry and wet deposition in the range of 65.8 ‰–68.0 ‰. The arithmetic mean value is 66.9 ‰, which we take to represent the δ¹⁸O of nitrate deposited from the atmosphere to the BHS.

The ammonium deposited from the atmosphere is assimilated by phytoplankton and is subsequently entrained into the N cycle via remineralization and nitrification or is nitrified directly in the water. Thus, the nitrified atmospheric ammonium is included here as a source bearing on δ¹⁵N of nitrate in the seawater. The ammonium deposition in BHS was 6.15 × 10⁹ mol yr⁻¹, which is more than the nitrate deposition of 3.42 × 10⁹ mol yr⁻¹ (Zhang et al., 2004).

The atmospheric ammonium has low δ¹⁵N values of −6.53 ‰ to −1.2 ‰ (Zhang et al., 2008; Chang et al., 2019). Given that the δ¹⁵N value of ammonium of the North China Plain is −1.2 ‰ ± 4.5 ‰ (Zhang et al., 2007) and that there is no obvious accumulation of ammonium in the surface layer in the observations, we assume that this isotope value is identical to the δ¹⁵N value of nitrified atmospheric ammonium.

The δ¹⁸O of nitrate from nitrification is roughly 1 ‰ higher than that of ambient H₂O (Sigman et al., 2009; DiFiore et al., 2009; Casciotti et al., 2007; K. L. Casciotti et al., 2008). The δ¹⁸O of H₂O in the BHS was reported as −0.67 ‰ ± 0.25 ‰ (n=10) (Kang et al., 1994; Wu, 1991), and thus the δ¹⁸O of the newly nitrified nitrate should be approximately 0.3 ‰.

4.2.4 Benthic fluxes

Nitrate diffusing from water to sediment

A latest number of benthic reactive nitrogen loss including denitrification and anammox for the BHS and northern YS is 3.5 × 10⁹ t N yr⁻¹ (Zhang et al., 2018). Combining the area of the BHS and assuming that 82 % of benthic nitrogen loss was by denitrification, the denitrification flux calculates to 10.1 × 10⁹ mol yr⁻¹. Globally, the sediment denitrification rate varies in the range of approximately 0.5 to 2 mmol m⁻² d⁻¹ (Devol, 2015), which is equivalent to 14.1 × 10⁹ to 28.2 × 10⁹ mol yr⁻¹ in the BHS. We assume that diffusion is not accompanied by isotope fractionation (Devol, 2015) with the result that δ¹⁵N and δ¹⁸O of nitrate diffusing into the sediment are the same as the nitrate pool in BHS (8.9 ‰ and 10.4 ‰, respectively).

Ammonium diffusing from sediment to water

The processes of nitrogen cycling in sediments are complex and variable (Lehmann et al., 2004). The degradation of organic matter, nitrification, and assimilation are acting under aerobic conditions, whereas denitrification, anammox, and dissimilatory nitrate reduction to ammonium (DNRA) are observed under anaerobic conditions. When organic matter is degraded in the surface sediments, part of the produced ammonium diffuses into the overlying bottom water and subsequently is nitrified to nitrite and nitrate under aerobic condition. For our purpose only the nitrified ammonium bears on the seawater nitrate pool. The mean δ¹⁵N value of sediment in BHS was 5.4 ‰ (n=20), and according to the fractionation factor during organic matter remineralization of 2 ‰ (Möbius, 2013) and subsequent nitrification (see above), the δ¹⁵N and δ¹⁸O of nitrate efflux from the sediment are assumed to be 3.4 ‰ and 0.3 ‰, respectively.

4.2.5 Sedimentation

The mass flux of N_r sedimentation is unknown. In terms of the effects of N_r sedimentation on nitrate dual isotopes, phytoplankton organisms that assimilate nitrate from the dissolved phase are the main source of sinking particles with the result that the N and O will be removed from the nitrate pool following the assimilation fractionation factor. Sinking particles in the BHS have a δ¹⁵N of 5.2 ‰ (δ¹⁵N_sink) corresponding to the average values of spring and summer, which integrates multiple processes such as photosynthesis of phytoplankton, heterotrophic synthesis of bacteria, and heterotrophic degradation (remineralization).

There are no observed data of δ¹⁸O of nitrate removed from the pool during assimilation (δ¹⁸O_sink), but this value can be estimated by the assimilation fractionation factor (¹⁸ε). The per mil fractionation factors ε of N (¹⁵ε) and O (¹⁸ε) in nitrate during assimilation are generally assumed to be around 5 ‰ so that ¹⁵ε : ¹⁸ε=1 : 1. Here we adopt the average of ¹⁵ε and ¹⁸ε as 5 ‰ (Granger et al., 2010; DiFiore et al., 2009; S. M. Liu et al., 2017; Wu et al., 2019; Umezawa et al., 2014; Wang et al., 2016) so that the δ¹⁸O of nitrate removed from the pool during assimilation (δ¹⁸O_sink) should be 5.0 ‰ according to the δ¹⁸O (10.0 ‰) of the dissolved nitrate pool.

4.3 The nitrate budget in the BHS

A box model of the nitrate budgets for the Bohai Sea following the LOICZ approach (Zhang et al., 2004) balanced sources and sinks of nitrate in BHS and was updated by several other nitrate budgets for the BHS during the last two decades (Zhang et al., 2004; Liu et al., 2003, 2009, 2011). All were, in general, not completely constrained because of a lack of data on some important source or loss terms. We here associate the nitrate isotope compositions of pools, sources, and sinks of nitrogen with a box model of the BHS nitrate in order to improve the understanding of nitrate cycling in the BHS. Finally, based on the combined mass and isotope box model informed by new data on the isotopic composition of nitrate, surface sediment, and suspended particulate nitrogen in the water column discussed above, we propose an updated N budget that is internally consistent.

4.3.1 The nitrate budget based on mass fluxes and corresponding δ¹⁵N values

The sources of nitrate for BHS are river inputs, submarine fresh groundwater input, atmospheric deposition, and remineralization. The most important sinks are net export to the YS, sediment denitrification, and particulate matter sedimentation. Assuming the mass and N isotope of nitrate in the BHS are in steady state, the sources and sinks of nitrate follow Eqs. (1) and (2):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& \left(m_{\mathrm{atm}}+m_{\mathrm{r}}+m_{\mathrm{N}}+m_{\mathrm{SFGD}}+m_{\mathrm{ntr}}\right)\\ & -\left(m_{\mathrm{net}}+m_{\mathrm{sink}}+m_{\mathrm{denitr}}\right)=\mathrm{0},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& ({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{atm}}{m}_{\mathrm{atm}}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{r}}{m}_r+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{N}}{m}_{\mathrm{N}}\\ & +{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{SFGD}}{m}_{\mathrm{SFGD}}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{ntr}}{m}_{\mathrm{ntr}})\\ & -({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{net}}{m}_{\mathrm{net}}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{sink}}{m}_{\mathrm{sink}}\\ & +{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{denitr}}{m}_{\mathrm{denitr}})=\mathrm{0},\end{array}\end{array}$$

where the terms m with different subscripts refer to the corresponding nitrogen mass fluxes, m_atm refers to atmospherically deposited nitrate, m_r refers to river nitrate, m_N refers to nitrified ammonium deposited from the atmosphere, m_SFGD refers to nitrate in submarine fresh groundwater discharge, and m_ntr refers to nitrification in the water column. In terms of sinks, m_net refers to the mass fluxes associated with net export of nitrate from BHS to the YS, m_sink refers to nitrate sedimenting from seawater as particulate N, and m_denitr refers to denitrification in the sediment. The unit of the mass fluxes is 10⁹ mol. The “δ¹⁵N” refers to the δ¹⁵N value of the N mass flux with the same subscripts. As mentioned previously, the mass fluxes for m_N, m_ntr, m_sink, and δ¹⁵N_ntr are unknown and need to be constrained.

The range of δ¹⁵N_ntr can be constrained by a simplified interior nitrate cycling model. Ammonium links particles and nitrate in this interior cycling, and there are two different sources of remineralized ammonium in seawater. One is ammonium diffusing from the sediment, and the other is the ammonification of PN in the water column. The ammonium from remineralization through both processes is then nitrified in the water column and is a source of nitrate. The average δ¹⁵N values of PN and sediments in the BHS in our study are 5.2 ‰ and 5.4 ‰, respectively. The fractionation factor of ammonification of PN and sediment as the first step of generating recycled nitrate is estimated at 3 ‰ (Sigman and Fripiat, 2019) and 2 ‰ (Möbius, 2013), respectively. The remineralized ammonium from PN and sediments thus should have δ¹⁵N values between 2.2 ‰ and 3.4 ‰.

The ammonium concentrations in the BHS are low in the water column in spring and in the surface layer in summer, indicating that the ammonium from PN mineralization is most likely completely converted to nitrate so that there is no fractionation effect for this step. Thus, the δ¹⁵N value of newly nitrified nitrate from complete nitrification of ammonium generated by PN mineralization is 2.2 ‰. In the case of incomplete nitrification, especially under the thermocline in summer, the newly nitrified nitrate has a δ¹⁵N of 0.2 ‰, given a net fractionation factor of nitrification (ammonium to nitrate) of 2 ‰ (Sigman and Fripiat, 2019). In our model below, the assimilation of ammonium originating from SPM remineralization was not included as its proportion is unknown in the BHS. The simplified model may thus underestimate the input of ¹⁵N-depleted nitrogen into the nitrate pool.

Ammonium diffusing out of the sediment will either be mixed into the euphotic layer and subsequently assimilated by the phytoplankton or nitrified in the water column. The accumulation of ammonium beneath the thermocline was a significant process in summer, as shown by the high nitrite and ammonium concentrations beneath the thermocline. The ratio of these two branching processes is not known in the BHS. If all of the ammonium from sediment is nitrified, the produced nitrate will have a δ¹⁵N of 3.4 ‰ (see above). If the ammonium is only partially nitrified (especially in summer beneath the thermocline), the produced nitrate will have a δ¹⁵N of 1.4 ‰ at a fractionation factor of nitrification of 2 ‰ (Sigman and Fripiat, 2019). Thus, the δ¹⁵N value of the nitrate produced by nitrification (δ¹⁵N_ntr) of ammonium from sediment is in the range of 1.4 ‰–3.4 ‰.

Overall, combining the δ¹⁵N ranges of PN- and sediment-originated nitrate, the range of δ¹⁵N for newly nitrified nitrate is 0.2 ‰–3.4 ‰.

4.3.2 The coupled N and O budget box model of nitrate

Because the nitrate mass fluxes m_N, m_ntr, and m_sink cannot be segregated only based on N mass budgets and δ¹⁵N values, we turn to the δ¹⁸O values of the sources and sinks of nitrate for further constraints. Equation (3) applies if we assume that the oxygen isotope composition of nitrate reflects the steady-state mass fluxes, as does δ¹⁵N of the nitrate pool:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& ({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{atm}}{m}_{\mathrm{atm}}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{r}}{m}_{\mathrm{r}}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{ntr}}{m}_{\mathrm{N}}\\ & +{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{SFGD}}{m}_{\mathrm{SFGD}}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{ntr}}{m}_{\mathrm{ntr}})\\ & -({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{net}}{m}_{\mathrm{net}}+{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sink}}{m}_{\mathrm{sink}}\\ & +{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{denitr}}{m}_{\mathrm{denitr}})=\mathrm{0},\end{array}\end{array}$$

where the δ¹⁸O subscripts refer to the nitrate mass flux with the same subscripts. The δ¹⁸O of different sources and sinks are either fixed values or ranges of values in our own data, or those taken from the literature (Table 1).

According to Eqs. (1), (2), and (3), the unknown mass fluxes m_N, m_ntr, and m_sink can be solved by a set of ternary linear equations including the three unknown terms, when appropriate boundary values of δ¹⁵N_ntr and δ¹⁸O_SFGD are chosen. δ¹⁵N_ntr ranged from 0.2 ‰ to 3.4 ‰, as discussed in Sect. 4.3.1, whereas δ¹⁸O_SFGD is a crucial term without any data or literature constraint. As the only constraint, δ¹⁸O_SFGD is expected to be higher than the value of nitrate imported from the rivers (δ¹⁸O_r=1.3 ‰) due to the fractionation associated with denitrification in the anaerobic aquifers (see Sect. 4.2.2). The results can be summed up in three different cases.

1. When setting the value of δ¹⁵N_ntr to 3.4 ‰, we obtain estimates for m_N that range from 0.00 to 4.83 × 10⁹ mol yr⁻¹, for m_ntr in the range of 15.08 to 32.27 × 10⁹ mol yr⁻¹, and m_sink in the range of 26.08 to 38.45 × 10⁹ mol yr⁻¹. The corresponding values of δ¹⁸O_SFGD range from 1.3 ‰ to 15.9 ‰, and the upper range of δ¹⁸O_SFGD yields m_N=0.0 due to the assumption that any mass flux must be equal to or greater than 0.

2. When we choose a δ¹⁵N_ntr value of 3.2 ‰–2.4 ‰ to explore effects of the methodological error of δ¹⁵N for our isotope method (0.2 ‰; see Sect. 2.2), again under the premise that the mass fluxes are positive numbers, the result is m_N, m_ntr, and m_sink estimates in a narrower range than when δ¹⁵N_ntr is 3.4 ‰; these results are not shown.

3. The δ¹⁵N_ntr values lower than 2.4 ‰ result in δ¹⁸O_SFGD lower than 1.3 ‰, which is highly unlikely.

Thus, reasonable solutions only are reached when δ¹⁵N_ntr is between 3.4 ‰ and 2.4 ‰. The results of the budget are shown in Fig. 10 and Table 1, respectively.

Figure 10 The budgets and the corresponding dual isotope values of nitrate in the BHS in 2018. The terms included are discussed in the text in Sect. 4.3. The unit of mass flux is 10⁹ mol yr⁻¹.

Table 1 Sources and sinks and the corresponding δ¹⁵N and δ¹⁸O values of nitrate in the BHS.

The mass flux of nitrate (m_N) originating from the nitrification of atmospheric ammonium ranges from 0.00 to 4.83 × 10⁹ mol yr⁻¹ and accounts for up to 79 % of the total ammonium deposition (6.15 × 10⁹ mol yr⁻¹; Zhang et al., 2004). This in turn implies that most of the atmospherically deposited ammonium is directly assimilated rather than nitrified to nitrate. This agrees with phytoplankton preference to assimilate ammonium rather than nitrate (Glibert et al., 2016). The bulk of internal sources of nitrate originates from nitrification in the water column (from water column ammonification and ammonium diffusing from sediment). This single source (m_ntr) accounts for 40.7 %–65.3 % of the total sources of nitrate in BHS and appears to be much more important than in other coastal environments. For example, between 15 % and 27 % of productivity was supported by nitrified ammonium in the seawater in Monterey Bay (Wankel et al., 2007). Likewise, nitrification supplied 34 % of the surface nitrate in eastern Hainan Island, which like Monterey Bay is also an upwelling area (Chen et al., 2020). This indicates that nitrate regeneration by nitrification may play a more important role in shallow and land-input-dominated marginal seas than in upwelling-dominated marine settings.

Table 2 Results of sensitivity tests on mass fluxes and isotopes.

4.3.3 Assessment of model uncertainties

The uncertainties of the modeled nitrate mass and isotope budget lie essentially in some of the mass flux estimates and in possibly erroneous assumptions on values of δ¹⁵N and δ¹⁸O in some of the sources and sinks. In the sensitivity test, the three critical mass fluxes m_N, m_ntr, and m_sink and the upper limit value of δ¹⁸O_SFGD change correspondingly when only one single term or one group of terms change (i.e., mass flux or isotope value). Note that the lower limits of m_N and δ¹⁸O_SFGD are fixed as 0 ‰ and 1.3 ‰, respectively, so that only the relative deviations of their upper limits are discussed below. The results are given in relative deviation to our best estimates in the budget. In addition, the widest ranges of mass fluxes and δ¹⁸O_SFGD are covered when δ¹⁵N_ntr=3.4 ‰ in the standard budget; thus we only compare the deviations when δ¹⁵N_ntr=3.4 ‰.

River input of nitrate is 9.49 × 10⁹ mol in the year 2018 as suggested above. Considering that the river fluxes are variable annually, we also adopt the multi-year average value (1952–2015) water discharge of the YR of 292.8 km³ instead of 333.8 km³ in 2018 (MWR, China, 2019). Hence the annual nitrate discharge is 8.52 × 10⁹ mol and is 10 % less than our preferred estimate; the resulting relative deviations are shown in Table 2. With respect to the isotope composition of the river nitrate, our observations suggest that the nitrate δ¹⁵N and δ¹⁸O of the YR changed little during the sampling period. On the other hand, the average values of nitrate δ¹⁵N and δ¹⁸O for rivers are mass weighted instead of arithmetic mean values so that a change in nitrate δ¹⁵N and δ¹⁸O for rivers with low nitrate discharges would induce little change in the average values. Considering the wide variations in the riverine dual nitrate isotopes around the BHS (Yu et al., 2021; Yue et al., 2013), we tested situations in which the dual nitrate isotope values of the smaller rivers are 0.5 and 2 times our observed values (Table 2: 0.5 δ¹⁵N_r and δ¹⁸O_r, 2 δ¹⁵N_r and δ¹⁸O_r). The results suggest that the budget is particularly sensitive to higher δ¹⁵N_r and δ¹⁸O_r than assumed in our preferred model.

We consider the mass flux of atmospheric deposition of nitrate to be reliable because differences in estimates from previous studies are quite small. According to these previous estimates (Zhao et al., 2017; Liu et al., 2003), we adopt the minimum and maximum in mass flux of 3.14 × 10⁹ mol to 3.65 × 10⁹ mol for uncertainty estimates. For δ¹⁵N_atm, we test the values 1.0 ‰ and 2.6 ‰ according to the variable proportion of wet/dry deposition with different δ¹⁵N endmembers. Likewise, δ¹⁸O_atm was varied in the range of 60.3 ‰–71.6 ‰.

The flux of SFGD into the BHS varies in the range of 2.68 × 10⁹ to 5.82 × 10⁹ mol, which is mainly due to the different SFGD water fluxes assumed in our budget. Only if m_SFGD is assumed to be >2.90 × 10⁹ mol do the results match our assumption. When m_SFGD=5.82 × 10⁹ mol, the relative deviations are large (Table 2), implying that our budget is very sensitive to m_SFGD.

The flux of denitrification in the BHS has a confidence interval of 15 % (Zhang et al., 2018), resulting in the m_dentr ranging from 8.57 × 10⁹ to 11.60 × 10⁹ mol. The results of tests with these two fluxes are shown in Table 2.

Overall, these tests indicate that m_N and δ¹⁸O_SFGD vary by ± 50 % and are particularly sensitive to uncertainties of the assumed endmembers, whereas the relative deviations of m_ntr and m_sink normally vary in the range of ± 20 %. Uncertainties will be significantly reduced if any of these terms can be constrained by further empirical studies.

4.3.4 Biogeochemical implications of the box model

In other coastal eutrophic regions, such as the North Sea, a high δ¹⁵N of river nitrate is reflected in a halo of high δ¹⁵N in surface sediments in offshore areas (Pätsch et al., 2010). In the Bohai Sea, such an isotopic halo of river-borne eutrophication is not observed despite similar water exchange rates of 1–2 years (Li et al., 2015; Serna et al., 2010) and similarly isotopically enriched river inputs. We speculate that the lack of a fingerprint of river nitrate in the δ¹⁵N of sediments of BHS may be masked by active nitrification and atmospheric deposition that rapidly eradicate and homogenize spatial gradients.

Despite the uncertainties that are related to the box model approach, combining mass and isotope budgets of nitrate sources and sinks is clearly superior to solely nitrate mass balance considerations, especially when it comes to segregating the anthropogenic nitrate and the recycled nitrate inputs. It is of note that the BHS does not appear to pass on significant amounts of nitrate to the adjacent northern YS with the result that the effects of excessive loading of this shallow mixing zone between land and ocean with anthropogenic nitrogen are mitigated by internal cycling processes.

5 Conclusions

Rivers contributed 19.2 %–25.6 % of the total nitrate input to the BHS, and the combined terrestrial runoff (including submarine discharge of nitrate with fresh groundwater) accounts for 27.8 %–37.1 %. Atmospheric input contributes 6.9 %–22.2 % of nitrate to the BHS. Nitrification contributes 40.7 %–65.3 % of the total nitrate, indicating an unusually active interior nitrogen cycling in the BHS. Nitrate was mainly trapped in the BHS, and only very little was exported to the YS (only 1.8 %–2.4 %). Furthermore, nitrate was rather assimilated than exported to the YS along the main transport pathway Lubei Coastal Current, effectively retaining N_r in BHS. Sedimentation trapped 70.4 %–77.8 % of nitrate inputs, whereas denitrification was only active in the sediments that removed 20.4 %–27.2 % of nitrate from the pool. Seasonal biogeochemical variations were observed in the BHS in that dissolved inorganic nitrogen increased during summer under the thermocline, implying significant biological regeneration. If the interior cycling increases, for instance fueled by increased terrestrial and atmospheric N_r inputs, respiration coupled to organic matter and N recycling will increase, and water column hypoxia could consequently spread in the future and compromise ecosystems in the BHS. Whether this will invigorate water column denitrification to balance the additional inputs is an open question, as is the capacity of BHS as a nitrate buffer between the growing source of N_r on land and the open ocean.

Appendix A

Nitrate was reduced to nitrite with a copperized cadmium column first. The nitrite ions reacted with sulfanilamide and N-1-naphthylethylendediamine (NEDD) to form red azo dye and then was measured at 520–560 nm. Phosphate determination followed the method of Murphy and Riley (Murphy and Riley, 1962). Under acid conditions a phosphomolybdic complex was formed of orthophosphate, antimony, and molybdate ions (Wurl, 2009). Followed by the reduction of ascorbic acid, the blue color complex was measured at 880 nm. The sample with ammonium is reacted with o-phthalaldehyde (OPA) at 75 ^∘C in the presence of borate buffer and sodium sulfite to form a fluorescent species proportional to the ammonia concentration. The fluorescence is measured at 460 nm following excitation at 370 nm (Kérouel and Aminot, 1997). Silicate is reacted with ammonium molybdate to silicomolybdate and reduced in acidic solution to molybdenum blue by ascorbic acid (Grasshoff et al., 2009).