The GBA, with a total land area of ∼56,000 km² (including nine cities in Guangdong province and two Special Administrative regions: HongKong and Macao) and a coastal length of ∼1,480 km, is located in the northern South China Sea (Figure 1). There are one important estuary in southern China (Pearl River Estuary) and many important bays including Daya Bay, Dapeng Bay as well as Tolo Harbor in HongKong. The Pearl River, with an annual runoff of 3.26 × 10¹¹ m³/yr and suspended particles matter (SPM) load of 5.0 × 10¹⁰ g/d, is the largest river into the South China Sea through eight outlets (Wang et al., 2021). During the dry winter season, the coastal current flows from northeast to southwest, and the Pearl River flow turns to the west. Rainfall in this region is abundant and occurs mostly (80%) during the wet summer season (from April to September).

Fieldwork was carried out during the dry winter season (January 6–13, 2020, Supplementary Table 1). The sampling sites within 30 seawater from three transects (A1, A2, and A3), 24 coastal groundwater, and 7 rivers are shown in Figure 1. Water samples for radium extraction were collected in a large volume (60 L for seawater, 30 L for river water, ∼5 L for coastal groundwater) and drained slowly (<1 L/min) through manganese-coated acrylic fiber (Mn-fiber) that absorb dissolved radium from water (Moore, 1976). Nitrate samples were filtered through 0.45 μm filters and collected with 45 ml Nalgene sampling bottles. DIC and DOC samples were filtered through 0.22 μm filter membranes, poisoned with HgCl and H₃PO₄ solution, respectively, and stored in 45 mL sampling vials. These samples were refrigerated at 4°C before analysis. The values of pH, oxidation-reduction potential (ORP), and salinity in water were immediately measured in situ using a multi-parameter water quality analyzer (HI9829T, HANA).

²²⁴Ra, ²²³Ra, and ²²⁸Ra were detected by radium delayed coincidence counting system (RaDeCC; Moore, 2008). The uncertainties for radium measurements were 7% for ²²⁴Ra and ²²⁸Ra, and 12% for ²²³Ra. Analyses for nitrate were carried out colorimetrically using a flow injection analyzer (FIA) with a detection limit of 5 μg/L (HACH QC8500, American). DOC and DIC were determined by a Total Organic Carbon/Nitrogen analyzer (Multi N/C 3100 Analyzer, Germany). The uncertainties of DOC and DIC concentration were better than 3%.

Radium isotopes (²²⁴Ra, ²²³Ra, ²²⁸Ra, and ²²⁶Ra), with the half-life from 3.6 days to 1,600 years, are widely used to trace SGD fluxes at different scales because of their enrichment in SGD (Rodellas et al., 2015; Wang et al., 2015; Zhang et al., 2016). Garcia-Orellana et al. (2021) summarized the application of radium isotopes tracing SGD and associated solute fluxes. Three models, radium mass balance model, radium endmember mixing model, and eddy diffusive mixing model, are usually applied to estimate SGD fluxes. The basic strategies of these models to estimate SGD is to determine first the radium flux supplied by SGD, and subsequently convert it into SGD flux by characterizing the groundwater radium endmember. Radium mass balance model which considers all the potential radium sources and sinks is the most widely used (Moore et al., 2008; Wang et al., 2018). The eddy diffusive mixing model used in this study is suitable for open coastal systems.

The long-lived ²²⁸Ra was used to assess SGD flux into the GBA. The total input of ²²⁸Ra in the system is mainly controlled by SGD and river discharges; the input from sediments is negligible due to the long half-life of ²²⁸Ra. Thus SGD-derived ²²⁸Ra (FSGDRa-228) can be calculated by subtracting riverine ²²⁸Ra (FRRa-228) from the total input (FTotalRa-228) (Eq. 1a). The ²²⁸Ra flux attributed to river was calculated as the product of the activity of ²²⁸Ra (CR-iRa-228) in the river and the corresponding river discharge (Q_R–i) (Eq. 1b). The total ²²⁸Ra input to the system was balanced by the offshore export of ²²⁸Ra, which was calculated as the product of the eddy diffusion coefficient (K_h), the offshore gradient of ²²⁸Ra (gSeaRa-228), and the cross-sectional area (A) of the system (Eq. 1c).

where K_h can be determined by an advection-diffusion model of short-lived ²²⁴Ra and ²²³Ra; n is the number of rivers in the region. The ²²⁴Ra and ²²³Ra distribution can be described as follows if ignoring the advection (Hancock et al., 2006; Wang et al., 2021):

where CSeaRa is ²²⁴Ra or ²²³Ra activity in seawater (dpm/100 L), x and D are offshore distance and seawater depth (m), respectively; l is the decay constant of ²²⁴Ra or ²²³Ra (d^–1), and B is the sedimentary radium flux [dpm/(m²d)]. The boundary conditions, the radium fluxes across the coastline (x = 0) to be F₀ and at the shelf edge (x = x_e) to be limx→xe∂⁡CSeaRa∂⁡x=0, were given to solve this equation. Integration of Eq. 2a with respect to x yields:

The values of K_h can be determined when we obtained optimal agreement between measured and modeled radium activities (²²⁴Ra and ²²³Ra) for each transect. SGD-derived nitrate or carbon flux (F_N) is commonly calculated by multiplying the SGD flux (Q_SGD) by the concentration of nitrate or carbon in coastal groundwater (Wang et al., 2020). It can be represented by the following equation:

where C_N and CGWRa-228 are solute (nitrate, DOC, and DIC) concentration and ²²⁸Ra activity in coastal groundwater, respectively.

The common method to estimate SGD-derived solute (nitrate and carbon) fluxes is to calculate first the SGD flux by dividing SGD-associated radium flux by radium activity in groundwater (i.e., radium endmember), and then multiply the SGD flux by solute concentration in groundwater (i.e., solute endmember). During this process, the determination of the two endmembers is usually major source of uncertainty for the final SGD-derived solute fluxes because of the spatial and temporal variability of radium and solutes in groundwater (Michael et al., 2011; Waska et al., 2019; Wang et al., 2020). To reduce the uncertainty of SGD-derived nitrate/carbon flux estimations, in this study, the term CN/CGWRa-228 in Eq. 3 was regarded as a unified parameter to discuss later.

The spatial distributions of radium (²²⁴Ra, ²²³Ra, and ²²⁸Ra) and salinity in seawater and coastal groundwater are shown in Figure 2. In general, the high radium activity and low salinity occurred in the nearshore area, whereas the low radium activity and high salinity occurred in the offshore area. Seawater salinity ranged from 14.49 to 34.99 PSU (Practical Salinity Unit) and it was very low in the western nearshore area resulting from the Pearl River freshwater dilution (Figure 2D). The radium activity of seawater ranged from 11.72 to 113.72 dpm/100 L for ²²⁴Ra, from 0.15 to 4.76 dpm/100 L for ²²³Ra, and from 8.04 to 80.41 dpm/100 L for ²²⁸Ra (Supplementary Table 1). Coastal groundwater, however, had very high radium activity within ²²⁴Ra: 0.81–119.2 dpm/100 L, ²²³Ra: 0.05–4.10 dpm/100 L, and ²²⁸Ra: 0.17–34.46 dpm/100 L (Supplementary Table 2). As shown in Figure 3, ²²⁴Ra had a strong linear relationship with ²²³Ra and ²²⁸Ra in both coastal groundwater and seawater. The coastal groundwater with high Ra activities and ratios (²²⁴Ra/²²³Ra and ²²⁴Ra/²²⁸Ra) was the main source of radium in the ocean.

Solute concentrations of coastal groundwater ranged from 0.005 to 3.17 mg/L with a median of 0.19 mg/L for NO₃^–, from 1.00 to 7.83 mg/L with a median of 1.34 mg/L for DOC, and from 4.58 to 31.49 mg/L with a median of 13.56 mg/L for DIC. For the river water samples, solute concentrations ranged from 0.69 to 4.36 mg/L with a median of 1.22 mg/L for NO₃^–, from 1.08 to 1.82 mg/L with a median of 1.23 mg/L for DOC, and 9.50 mg/L for DIC (Figure 4 and Supplementary Table 3). The highest concentration of NO₃^– (4.36 mg/L) and DOC (1.82 mg/L) occurred in the Dan’ao River (R7) which suffered serious environmental pollution from high-intensity human activities (Wang et al., 2021). Solute concentrations in both coastal groundwater and river water are ranked as DIC> DOC> NO₃^–. Nitrate concentrations in coastal groundwater are less than that in river water, while carbon contents in coastal groundwater are slightly greater than that in river water (Figure 4).

The nitrate and carbon fluxes via rivers into the GBA were calculated as the product of the concentrations of nitrate and carbon in the river and the corresponding river discharge. In the GBA, the Pearl River is the dominant river with a flux of 8.92 × 10⁸ m³/d and a mean concentration of 1.01 mg/L for NO₃^–, 1.18 mg/L for DOC, and 9.50 mg/L for DIC. Thus, the riverine NO₃^–, DOC, and DIC inputs were estimated to be 8.82 × 10⁸ g/d, 1.05 × 10⁹ g/d, and 8.48 × 10⁹ g/d, respectively. Numerous studies have well-investigated cycles of nitrogen and carbon and their seasonality within the Pearl River and adjacent waters (Lu et al., 2009; He et al., 2010; Liu et al., 2020; Xuan et al., 2020; Ye et al., 2021). Liu et al. (2020) investigated the spatial distribution of organic carbon and nitrogen in the Pearl River. They found that both DOC and total organic nitrogen (TON) were increased from upstream to downstream, with an increase of five times for TON and two times for DOC in the dry season. The nitrogen and carbon fluxes from the current study are consistent with previous findings.

Wang et al. (2021) quantified nitrate fluxes from SGD and rivers into the GBA in the wet summer season with values of (5.28 ± 0.73) × 10⁸ g/d and 1.32 × 10⁹ g/d, respectively. It can be found that nitrate fluxes from SGD and river to the GBA were greater in the wet summer season than in the dry winter season. This is the less discharge of river and groundwater during the dry season, which can be supported by the seawater salinity and pH distributions in the seawater. Figure 7 shows seawater salinity and pH in the dry and wet seasons. Both salinity and pH were less in the wet summer season than in the dry winter season, which indicated directly the decrease of terrestrial freshwater including river water and fresh groundwater.

The eddy diffusion coefficients for the GBA in the summer season were estimated by Wang et al. (2021) using the same method. The summer eddy coefficient was 21.04 km²/d for A1 transect, 88.4 km²/d for A2 transect, and 11.76 km²/d for A3 transect (Wang et al., 2021). The winter eddy coefficient determined in this study is about 4–10 times greater than the summer eddy coefficient. During the winter, the strong current will enhance horizontal mixing and result in relatively higher values for eddy mixing. Among three transects (A1–A3), A2 has the highest values of eddy diffusion coefficient during both summer and winter seasons. Transect A2 is the longest one and is more influenced by terrestrial freshwater (Pearl River). The force of Pearl River discharge could enhance horizontal mixing and thus create a relatively higher value of eddy diffusion coefficient. In addition, differences in scale can produce some of the variances in the eddy diffusion coefficient. The eddy diffusion term usually becomes larger at greater length scales (Moore, 2000).

Regarding carbon flux, it can be found that SGD-derived DOC and DIC fluxes are ∼2 times as great as riverine input, but SGD-derived NO₃^– flux is one-fourth of the riverine input. Anthropogenic input of nitrogen is the primary determinant of NO₃^– accumulation in the Pearl River (Xuan et al., 2020). DIC is the dominant form of carbon in SGD, and DOC concentration is below one-tenth of DIC in most groundwater samples (except for two sites, Figure 8A). The relationship between DOC and the residual NO₃^– is shown that DOC concentrations decreased while dissolved NO₃^– concentrations increased (Figure 8B). This observation revealed that DOC could be consumed by biological production when they transport from land to sea. Waska et al. (2021) also found that the subterranean estuary is a net sink of terrestrial organic carbon. SGD-derived DIC and DOC have been reported in some coastal systems worldwide such as Okatee Estuary, South Carolina, United States (Moore et al., 2006), Southwest Florida Shelf, United States (Liu et al., 2014), and Northern South China Sea (Tan et al., 2018). These findings suggested that SGD is a dominant source of DIC and DOC to coastal waters and may contribute significantly to coastal carbon budgets.

Submarine groundwater discharge-derived and riverine NO₃^– into the GBA can stimulate new primary production and enhance biological pump efficiency. Using the Redfield ratio of C/N = 106:16 to derive primary production in terms of carbon uptake, estimated potential primary production by SGD and river were 1.17 × 10⁹ g C/d and 5.01 × 10⁹ g C/d, which were 7.1 and 59.1% of SGD-derived and riverine DIC, respectively. The new primary production rates supplied by SGD showed a similar range to the study in the Daya Bay [i.e., 54–73 mg C/(m²d)], an important bay in the GBA (Wang et al., 2018).

Submarine groundwater discharge-derived and riverine DOC, once entering the coastal sea within the environment of sufficient inorganic nitrogen, could be easily consumed by biological production and thus cannot be stored as inert organic carbon. If assuming that 80% of the flux of SGD-derived and riverine DOC will be remineralized, the remineralization process can produce CO₂ of 1.80 × 10⁹ g C/d and 0.84 × 10⁹ g C/d, respectively. The remineralization of SGD-derived DOC usually occurs at the bottom of the ocean which will potentially result in hypoxia and acidification at the bottom water body. Combined with the processes of nitrate and DOC, SGD is a potential net source of atmospheric CO₂ with a flux of 0.63 × 10⁹ g C/d, while river is a potential net sink of atmospheric CO₂ with a flux of 4.17 × 10⁹ g C/d. Previous studies in the Pearl River Estuary (Wu et al., 2017) and other river-dominated systems (Bauer et al., 2013; Guo et al., 2015) have often observed the net production of DOC in the river plume.

The input of land DIC mainly affects the balance of the carbonate system in marine water. Compared with seawater, SGD and river water have lower pH values, which suggests that much more CO₂ could be dissolved in SGD and river water. Charbonnier et al. (2022) also found that terrestrial groundwater contains a lot of CO₂ and is an important contributor to the DIC flux to the coastal ocean. When SGD and river discharge to the ocean, the balance of the carbonate system may be changed resulting in enhancement of calcium carbonate or CO₂ exchange with the atmosphere. If assuming that the free CO₂ in SGD-derived and riverine DIC can escape to the atmosphere from water body. Based on pH values (Supplementary Table 2), the free CO₂ was estimated to be 5% of DIC. Thus, SGD-derived and riverine DIC would potentially supply atmospheric CO₂ of 0.83 × 10⁹ g C/d and 0.42 × 10⁹ g C/d, respectively.

Indeed, marine carbon cycles are influenced by the combination of terrigenous NO₃^–, DOC, and DIC. Based on the above discussion, we conclude that SGD in the studied system is a potential net source of atmospheric CO₂ with a flux of 1.46 × 10⁹ g C/d and river is a potential net sink of atmospheric CO₂ with a flux of 3.75 × 10⁹ g C/d. Their coupled effect on ocean suggests that SGD and river would potentially enhance CO₂ exchange flux from atmospheric to ocean in the dry winter season and they play a significant role in the carbon cycles of GBA.

Submarine groundwater discharge and river water usually have different physical and chemical characteristics, such as the long water residence times, strong water-rock interaction, and complex biogeochemical processes in groundwater flow. Therefore, the composition and structure of chemical materials in river water and groundwater discharged into the ocean are quite different. Generally, the discharged river water has high oxygen and nutrient concentrations and forms freshwater plumes on the ocean surface. Groundwater discharged into the sea from the seabed with a low flow rate, on the other hand, has a low oxygen content and high nutrient concentrations (Slomp and Van Cappellen, 2004; Kim et al., 2017; Waska et al., 2021). Thus they present different processes in marine carbon cycles. Here two conceptual models, as shown in Figure 9, were proposed to illustrate the major potential processes of the carbon cycle induced by SGD and river discharges.

Figure 9A shows the contributions of SGD and river discharges to the carbon system on the nearshore scale with shallow oceans. In this case, the bottom water can easily exchange with upper water and air. Thus CO₂ produced from SGD processes (e.g., remineralization of SGD-derived DOC) can escape to the atmosphere. On the other hand, the consumption of CO₂ in SGD processes (e.g., new primary production by SGD-derived nitrate) can achieve from the atmosphere. Other products from SGD can participate in other processes of the carbon cycle in shallow seawater. Overall, a significant amount of nitrate, DOC, and DIC stimulate phytoplankton blooms in the upper water and consume oxygen in the bottom water, resulting in coastal hypoxia (Su et al., 2017).

Figure 9B shows the contributions of SGD and river discharges to the carbon system on the shelf scale with deep oceans. In this case, the bottom water is in a reducing environment and hardly exchange with upper water and air. Thus the high-carbon water masses usually occur in the bottom resulting from the input of SGD-derived nitrate and carbon. The carbon in these high-carbon water masses deposited or moved into deeper waters and finally deposited and sealed. This demonstrated that SGD can be regarded as a sink of carbon in a sense. In the upper water, only a fraction of the organic carbon produced by the river was transported to the deep ocean and stored there.