Basin-scale control on N2O loss rate and emission in the Changjiang River network, China

Wang, Fang; Yu, Qibiao; Yan, Weijin; Tian, Siyu; Zhang, Peipei; Wang, Jianing

doi:10.3389/fmars.2022.1025912

ORIGINAL RESEARCH article

Front. Mar. Sci., 23 September 2022

Sec. Coastal Ocean Processes

Volume 9 - 2022 | https://doi.org/10.3389/fmars.2022.1025912

This article is part of the Research Topic Nitrous Oxide Production Processes and Associated Mechanisms in Estuarine and Coastal Ecosystems View all 5 articles

Basin-scale control on N₂O loss rate and emission in the Changjiang River network, China

Fang Wang^1†

Qibiao Yu^2†

Weijin Yan^1*

Siyu Tian^1,3

Peipei Zhang^1,3

Jianing Wang⁴

¹Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Science and Natural Resources Research, Chinese Academy of Sciences, Beijing, China
²College of the Environment and Ecology, Xiamen University, Xiamen, China
³College of Resource and Environment, University of Chinese Academy of Sciences, Beijing, China
⁴Chinese Academy for Environmental Planning, Beijing, China

Global riverine N₂O emissions have been made by several studies with great uncertainty. However, the regional N₂O budgets and patterns in large river networks is still unclear, due to the lacking understanding of in-river N₂O emission rate and well-classified river network water areas. Furthermore, the mass ratio of N₂O emission against nitrogen(N) load in river networks remains controversial. Here we report N₂O emissions from the largest river of China, the Changjiang River network, emphasizing the basin-scale control on riverine N₂O loss rate in response to increasing N loads and river size. We find the N₂O emission rate is negatively related to Strahler river orders, and positively related to N loading. The velocity (V_f) of N conversion into N₂O was 0.131-0.436 m yr^-1, and N₂O loss rate (ζ) was 0.27-37.64 ×10^-4 d^-1 and declined exponentially with water discharge. Both the loss rate and the mass ratio of N conversion into N₂O varied significantly at basin scale, indicating the diminishing capacity of river ecosystems to convert excess DIN into N₂O when N load increased as a direct result of human activities. Our study shows N₂O emission was 0.66 Gg N₂O-N (1Gg=10⁹g) in 1986 and increased to 10.3 Gg N₂O-N in 2014 for the whole Changjiang River network. We identified the headwater streams are hotspots of N₂O emission across the headwater stream to the estuary aquatic continuum. N₂O emission was about 0.82% - 5.31% of global riverine N₂O budget during 2010-2014. Our study suggested that an integrated approach in view of the riverine N loads and river hydrology is needed to improve estimates of riverine N₂O emissions.

Introduction

Nitrous oxide (N₂O) emitted from rivers and streams is a key component of both the global and regional N₂O budgets (Beaulieu et al., 2011; Borges et al., 2015; Turner et al., 2015; Hu et al., 2016) owing to their favorable hydrological and biogeochemical conditions in conversion of bioavailable nitrogen (bio-N, defined as dissolved inorganic nitrogen, DIN, mainly including ammonium and nitrate) into reduced N through the processes as nitrification, denitrification, or/and nitrifier-denitrification (Seitzinger et al., 2006; Mulholland et al., 2008; Kool et al., 2011), during which N₂O is produced. Some studies showed that bio-N inputs to rivers and streams significantly increased in the last few decades (Mayorga et al., 2010; Seitzinger et al., 2010; Yan et al., 2010; Nevison et al., 2016; Battye et al., 2017). For example, DIN input increased from 3.01 million t/yr in 1970 to 13.14 million t/yr in 2003 in Changjiang River watershed (Yan et al., 2010), from 1.49 million t/yr in 1986 to 3.58 million t/yr in 2015 in Zhujiang watershed (Cui et al., 2020), and from ~6000 t/yr in 1980 to ~20000 t/yr in 2010 in Jiulong River watershed (Yu et al., 2015). Thus, there is the potential to markedly stimulate N₂O production in river networks (Seitzinger, 1988; Bouwman, 1996; Yan et al., 2012). However, the response of N₂O emissions to increasing DIN load is still unclear within a large river network. Although the Intergovernmental Panel on Climate Change (IPCC) methodology has recommended an emission factor (EF) to facilitate estimation of N₂O emissions from rivers and streams (EF_5-r) (IPCC, 2006), uncertainty still remains (Beaulieu et al., 2011; Hinshaw and Dahlgren, 2013; Turner et al., 2015; Wang et al., 2015a).

In river networks, N₂O can be produced in situ through nitrification and denitrification in river sediments and water column (Wang et al., 2017; Li et al., 2019). Nitrification, the oxidation of ${NH}_{4}^{+}$ to nitrate ${NO}_{3}^{-}$ ), releases N₂O as a by-product, whereas denitrification ( ${NO}_{3}^{-}$ reduction) can either produce N₂O or transform it to N₂ in a final reaction favored under highly reduced anoxic conditions. These processes may be influenced by river hydrological characteristics, such as river depth, width, and river retention time (Seitzinger et al., 2006; Zarnetske et al., 2011). As river networks receive increasing DIN loads, the mass of DIN converted to N₂O emissions from water per unit area of riverbed per unit time, responds differently to DIN increase and hydrological conditions. Here, we define conversion velocity (V_f, m yr^-1) as the speed of DIN converted through N₂O emissions from water per unit area of riverbed per unit time. We further define the N₂O loss rate (ζ) as combination of the conversion velocity and river depth, and we hypothesis that both V_f and ζ can change in response to DIN loads and river sizes. For example, small streams generally have high ratios of riverbed area to water volume, and can efficiently convert DIN to N₂O, whereas, larger rivers have lower ratios of riverbed area to water volume making them less efficient in converting DIN to N₂O. There are limited studies of riverine N₂O emission in relation to hydrological processes in a large river network (Marzadri et al., 2014; Turner et al., 2015). Until now, researches have focused on a single river/stream or relatively small watersheds (Rosamond et al., 2012; Yan et al., 2012; Hinshaw and Dahlgren, 2013; Venkiteswaran et al., 2014) where river hydrological and watershed land use conditions are relatively homogeneous (Turner et al., 2016). However, little is known how N₂O emissions vary across a range of river sizes (Borges et al., 2015). Furthermore, quantifying N₂O emissions from large river networks is challenging due to the lack of data on three key variables; (1) in-stream N₂O production, (2) water surface area of the river network, and (3) N₂O transfer velocity.

Quantifying the amount of water surface area for a large river network is critical to evaluate riverine N₂O emission. River length and width can be obtained from the combination of publicly available geographic data and hydraulic geometry laws (Raymond et al., 2013), and river surface area can be calculated by the product of river length and width. Gas transfer velocity (k) is an important factor in N₂O emissions calculations. Models based on wind speed are widely applied, especially in open waters such as wide rivers, lakes, and estuaries (Cole and Caraco, 2001; Borges et al., 2004). Large river systems are not suited to wind-based models due to the challenge of measuring relative parameters across large scales. In contrast, measures of channel slope and flow velocity are easily obtained (Raymond et al., 2012), having been widely used in large space river gas emission (Butman and Raymond, 2011; Raymond et al., 2013; Borges et al., 2015).

Here, field sampling was taken for spatial, seasonal and diurnal study, along the main stem and two branches of the Changjiang River Network (CRN). We then estimated N₂O production through a Michaelis-Menten equation that described the nonlinear relationship between riverine N₂O concentration ([N₂O]) and DIN concentration ([DIN]) based on the field observation in the CRN. We also calculated water surface area of river network k based on the River Strahler model by combining with hydraulic geometry laws and geographic information system tools (Strahler, 1957). We estimated N₂O emissions for the temporal (period between 1970 and 2014) patterns from the entire CRN by combination of Global NEWS2-DIN (Yan et al., 2010). Finally, we used “nitrogen spiraling” theory to scale our analysis from field observations to the entire CRN (Workshop, 1990; Ensign and Doyle, 2006), and thereby obtained the temporal change of ratios of N₂O:DIN, representing dynamic response of EF at any given [DIN] in linking to DIN loads and hydrological processes at watershed scale and calculated EF values.

Materials and methods

Study area and sampling

The study was conducted in the Changjiang River basin (24°27’-35°54’N, 90°13’-122°19’E), the third largest river in the world, with an average annual discharge of 10× 10¹¹ m³. The Datong hydrological station (DHS), located at the upper limit of river mouth, is free from the influence of both tidal cycles and industrial waste from nearby cities. The stretch of the Changjiang River above DHS drains 94% of the whole basin and delivers more than 95% of the water to the estuary (Yan et al., 2003).

Field sampling was taken from the main stem and its tributaries in the Changjiang River network (CRN, Figure 1 and SI Table S1), including two sub-water systems in the lower reach of the CRN, the Poyang Lake-Ganjiang River water system and the Chaohu Lake water system. A total of 12 sampling sites were chosen to take samples in the CRN of which 4 sites along the main stem of the CRN; 4 sites located in the Poyang Lake-Ganjiang River water system; 2 sites in the Chaohu Lake water system; 1 site in the Jialing River branch, and 1 site in the Hanjiang branch (Figure 1 and SI Table S1). Spatial samples were taken from October to November 2011 at 4 hydrological stations along the main stem (Cuntan, Wanxian, Hankou and Datong hydrological stations, representing the 8^th order river) and 2 branch stations at the Beibei Hydrological station (Jialingjiang Branch River) and Xiantao hydrological station (the Hanjiang Branch River), representing the 7^th order river. The diurnal samples at the DHS were collected on 2 occasions in September 2015 and January 2016, respectively. In addition, diurnal sampling was also taken at 4 sites in the Poyang Lake-Ganjiang River water system during the summer season 2013, of which 3 sites located in the west, middle and east branch of the Ganjiang River (6th order river) and 1 site in the Xiushui River (5^th order river) (Figure 1 and SI Table S1). Monthly observations were sampled at 2 sites in the Chaohu Lake water system from June to October in 2009, 1 site located at the Pai River (2^nd order river) and the other one at the Hangbu River (3^rd order river). All the diurnal samples were collected at an interval of 6 h.

FIGURE 1

Figure 1 Location and river network of the Changjiang River basin, (A) showing the river network and hydrological stations, (B) all sampling sites, (C) sampling sites in Chaohu Lake watershed, and (D) sampling sites in Poyang Lake watershed.

At each sampling, surface water (~0.5 m) was collected in triplicate by filling in 60-mL glass serum bottles from water sampler, and water flow rate, water and air temperature, wind speed, and dissolved oxygen (DO) were measured in-situ. Samples were preserved by adding a small crystal of KOH or HgCl₂ solution (5% w/v) to each bottle immediately before sealing (Yan et al., 2010; Wang et al., 2015a). Surface water was also reserved in 200-mL glass bottle after filtrating through a 0.7-μm glass microfiber filter (GF/F Waterman, UK) for testing dissolved inorganic nitrogen (DIN, referring to nitrate and ammonium) concentrations. Dissolved N₂O concentrations were measured using headspace equilibrium method, that is, in headspace, water was firstly displaced with high pure Nitrogen (N₂) or Helium (He) gas with syringe, shaking for >1 hour, and then tested N₂O concentration of headspace by gas chromatograph (GC-2014, Shimadzu, Kyoto, Japan) or Agilent 7890A (Agilent, US) equipped with electron capture detector (ECD), with an analytical precision of 0.1 ppm (Wang et al., 2015b). Finally, dissolved N₂O concentration was calculated by using Weiss equation (Weiss and Price, 1980). Samples for nitrate and ammonium concentrations were immediately determined in lab by using Flow Injection Analyzer (FIA-3100, Beijing, China). In this study, nitrite concentration was neglected, because of this kind of DIN was relative lower comparing to nitrate and ammonium. Water temperature and other parameters were measured using a portable meter (YSI 550A, Yellow Springs, US).

Estimation of the surface area for the CRN

River surface area was computed by the product of width and length. The length of the rivers in the CRN was extracted from a 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM data set (Source: www.usgs.gov) in Arcgis 10.2 (ESRI, Inc, Redlands, USA) by Strahler river order method (Strahler, 1957). Besides, watershed area of each river was extracted. As a result, seven river orders were extracted directly, whereas first order streams were not captured due to the limit of resolution of the DEM data set. The length of smallest stream, first order river, for each subbasin was interpolated by empirical models, which were established by the regression relationship based on the extracted second to eighth orders.

River width was calculated from discharges according to hydraulic geometry scaling theory (Leopold and Maddock, 1953). Thus, velocity (V), width (W) and depth (D) are related to discharge (Q) with the following power relationships:

\begin{array}{l} V = a \times Q^{b} & (1) \end{array}

\begin{array}{l} W = c \times Q^{d} & (2) \end{array}

\begin{array}{l} D = e \times Q^{f} & (3) \end{array}

where a, c and e are hydraulic coefficients and b, d and f are hydraulic exponents. According to Leopold & Maddock (Leopold and Maddock, 1953), discharge is the product of velocity, width and depth. In this study, exponents and coefficients were fitted by hydrological station observation data of the selected 141 hydrological stations within the Changjiang River Basin (Figure 1, SI Figure S1, and SI Table S2).

The discharge of each river was computed from river watershed area and discharge yield, as:

\begin{array}{l} Q_{c a l} = Y_{Q y i e l d} \times S_{c a t c h m e n t} & (4) \end{array}

where Q_cal is the calculated discharge (m³ s^-1), Y_Qyield is the annual discharge yield (m yr^-1), and S_catchment is the river catchment area (km²).

To determine annual discharge yield of each subbasin, annual precipitation data from national meteorological stations (Source: http://data.cma.cn) and annual discharge yields from hydrological stations within the Changjiang River basin were collected. We took the mean precipitation and annual discharge yield for each subbasin and found that there exited a good regression relationship between them (R² = 0.68, n = 288, SI Figure S2a). Then we used this relationship to regress Y_Qyield from precipitation.

Determination of nitrous oxide gas transfer velocity

Nitrous oxide gas transfer velocity (k_N2O) determines the exchange volume between river and surface air. In this study, k_N2O is calculated from k₆₀₀ (coefficient of diffusion from water to the air for gas at 20°C), based on river channel hydraulic and basin characteristics according to Raymond et al., (Raymond et al., 2012), as:

\begin{array}{l} K_{600} = S_{s l o p e} \times V \times 2841 + 2.02 & (5) \end{array}

where S_slope is the slope of the river channel, and V is river flow velocity (m s^-1).

Channel slope was calculated as the difference in elevations divided by the river channel length between two hydrological stations located at the same river and river order of CRN. The channel elevation and river length data were collected from the Annual Hydrological Report (China, 1970-2012). A total of 120 pair fitted stations varying from the second- to the eighth- order were obtained (SI Figure S3a). For the smallest rivers, without available direct data, channel slope was extrapolated from the regression relationship among the second- to the eighth- order rivers (SI Figure S3a). River flow velocity was obtained from equation (1).

k_N2O between river surface water and the above air could be then determined by using:

\begin{array}{l} K_{N_{2} O} = K_{600} \times {(\frac{S c_{N_{2} O}}{600})}^{- \frac{2}{3}} & (6) \end{array}

where, k₆₀₀ is the coefficient of diffusion for N₂O at 20°C, Sc_N₂O is the Schmidt number for N₂O related to water temperature (Wanninkhof, 1992).

We first established the regression between river water temperature and air temperature observed from hydrological stations and national meteorological stations (R² = 0.92, n = 1224, SI Figure S2b). Then, annual river water temperature was interpolated by this regression relationship from air temperature.

Calculation of N₂O emissions

The air-water N₂O emission rate was calculated as

\begin{array}{l} N_{2} O_{E R} = ([N_{2} O_{w}] - [N_{2} O_{e q u}]) \times k_{N_{2} o} & (7) \end{array}

Where N₂O_ER is the N₂O emission rate (mg N m^-2 h^-1), [N₂O_w] is the dissolved N₂O concentration (μg N L^-1), [N₂O_equ] is the equilibrium N₂O concentration (μg N L^-1), ([N₂O_w] - [N₂O_equ]) can usually be expressed as Δ[N₂O] (μg N L^-1).

For the river network, the total N₂O emission from the water surface was calculated as:

\begin{array}{l} E f f l u x = S A \times N_{2} O_{E R} & (8) \end{array}

where Efflux is the annual N₂O emission (Gg N₂O-N yr^-1), SA is the surface area of river network (km²), and N₂O_ER is the N₂O emission rate (mg N m^-2 h^-1).

Relationship between observed N₂O and DIN concentrations in the CRN

We used a nonlinear Michaelis-Menten kinetics equation (MM) to determine the relationship between riverine DIN and dissolved N₂O concentrations from field sampling data (Boyer et al., 2006; Turner et al., 2016). The relationship was expressed as:

\begin{array}{l} [N_{2} O] = \frac{V_{\max} \times [D I N]}{K_{m} + [D I N]} & (9) \end{array}

where riverine dissolved N₂O concentration ([N₂O], mg N L^-1) was set as a function of DIN concentration ([DIN], mg N L^-1); V_max (mg N L^-1) represents the maximum reaction rate produced by microbial activates; K_m (mg N L^-1) is the affinity constant corresponding to [DIN] at half V_max.

Riverine DIN data sources

In order to analysis longitudinal variation of [DIN] in CRN, data from annual hydrological reports from 1970 to 1985 (China, 1970-2012) and public literature from 1999 to 2012 (Han and Liu, 2004; Chetelat et al., 2008; Li et al., 2008; Wu et al., 2008; Li et al., 2009; Wang and Wang, 2009; Han et al., 2010; Lu et al., 2011; Xu et al., 2011; Zheng et al., 2011; Ji and Jiang, 2012; Jiang et al., 2012; Yan et al., 2012; Liu et al., 2013; Sun et al., 2013; Wu et al., 2013; Han et al., 2014; Li et al., 2014a; Li et al., 2014b; Luo et al., 2014; An et al., 2015; Qu et al., 2015; Wang et al., 2015a) within the Changjiang River basin were collected (SI, Figure S4). For the [DIN] observed at hydrological stations, the sampling frequency ranged from 2 to 3 times per month. Descriptions of chemical analysis and data quality were similar to that reported by Yan et al., (Yan et al., 2003). A total of 102 stations in the Changjiang drainage basin were compiled and categorized by Strahler river order. Among these stations, 4 stations were located along the main stem of the river, including (from upstream downwards) the Cuntan, Yichang, Hankou, Datong hydrological stations (Figure 1). The other 98 stations were scattered along the Changjiang’s numerous tributaries. No stations in the delta region were included in this study because of the intersections with numerous streams and rivers originated from other drainage basins. In the CRN, nitrite concentrations were usually under test limit or much lower than nitrate and ammonium. The DIN in this study was only considered as the sum of nitrate and ammonium.

As nitrogen data were not routinely measured at hydrological stations since 1986, we modeled riverine [DIN] from the Global NEWS2-DIN (Mayorga et al., 2010). The summary of NEWS2-DIN model could be expressed as:

\begin{array}{l} D I N_{y i e l d} = (1 - Q_{r e m}) \times (1 - L_{d e n}) \times (1 - D_{d i n}) \times (D I N_{s e w} + F E_{w s} \times T N_{d i f f}) & (10) \end{array}

where DIN_yield is the modeled annual DIN yield per river basin area (kg N km^-2 yr^-1); Q_rem is the fraction (0-1) of consumption water from river; L_den is the fraction of nitrogen removed by denitrification in streams and rivers; D_din is the fraction (0-1) of nitrogen retained in the reservoir; TN_sew is the DIN exported from point source of industrial wastewater and domestic sewage (kg N km^-2 yr^-1); FE_ws is the fraction (0-1) of DIN exported from landscape to the rivers, and calibrated from observation data (Yan et al., 2010; Wang et al., 2014a); TN_diff is the total diffuse source of nitrogen input watershed (kg N km^-2 yr^-1), including nitrogen addition from chemical fertilizer, atmospheric deposition, animal manure nitrogen, and biological N-fixation. The model presents a sound prediction accuracy for [DIN] for the period of 1970 to 2003 (Yan et al., 2010), and extended to model [DIN] to 2013 (Wang et al., 2015b). Therefore, we updated the modeled [DIN] to the period of 2014 in the following analysis.

Calculation of the mass ratio of N₂O emission and DIN load based on nitrogen “spiraling” metrics

The mass ratio of emitted N₂O against DIN that entering the river network can be calculated by nitrogen “spiraling” theory (Workshop, 1990; Ensign and Doyle, 2006) (SI Table S3). Details of model structure and parameterization are presented in the SI Table S3.

First, areal rate of total DIN conversion velocity (V_f, m yr^-1) (Royer et al., 2004; Earl et al., 2006), defined as the speed of DIN converted through N₂O emissions from water per unit area of riverbed per unit time. V_f can be calculated as following:

\begin{array}{l} V_{f} = \frac{N_{2} O_{E R}}{C} & (11) \end{array}

where V_f is the N₂O conversion velocity (also DIN uptake velocity) (m yr^-1), N₂O_ER is N₂O emission rate (mg N m^-2 h^-1), and C is the [DIN] (mg N L^-1). We calculated V_f for each stream and river from 2^nd to 8^th river order in each subbasin, and estimated total values for 1^st order streams of each subbasin.

The N₂O loss rate (ζ), also called the first-order reaction rate, can be calculated based on equation (12):

\begin{array}{l} ζ = \frac{V_{f}}{D} & (12) \end{array}

where D is the mean river depth of stream or river (m).

Then, the mass ratio of N₂O emission and DIN load for each river order, also called the DIN removal fraction (R_RF) from the river water column. R_RF is related to N₂O loss rate and DIN retention time in rivers, so it can be expressed as:

\begin{array}{l} R_{R F} = 1 - e^{- ζ \cdot τ} & (13) \end{array}

where τ is the DIN retention time (d). The R_RF values of different order rivers are different.

Finally, the total mass ratio of N₂O emission against DIN load for the entire CRN can be calculated as:

\begin{array}{l} T R_{R F} = 1 - (1 - R_{R F 1}) \times (1 - R_{R F 2}) \times \dots \times (1 - R_{R F i}) & (14) \end{array}

Where, TR_RF is the total N removal fraction as N₂O emission, R_RF1, R_RF2…R_RFi is the removal fraction as N₂O emission for each river order, and i (i = 1,2,…,8) represents river order.

In this study, for each stream and river of each subbasin, the N₂O emission rate, V_f, and ζ were calculated; in each subbasin, the N₂O emission amount was estimated based on the N₂O emission rate and surface area of river network; for the whole basin, the total N₂O emission amount was the sum of all sub-basins involving all rivers from 1st to 8th orders.

Statistical analysis

In this paper, all the regression and nonparametric analysis and data figures were performed using MATLAB software (R2015b, MathWorks, Natick, USA). The regression parameters for hydraulic geometry scaling, discharge yield, water temperature, and channel slope were conducted using built-in fitting function regress under confidence level of 95%. The parameters in Michaelis-Menten kinetics equation was fitted using a nonlinear least square fitting function mimen under confidence level of 95% running in MATLAB software.

Results

Riverine surface area of the CRN

We used data of observed river discharge, velocity, width, and depth for multiple years (1970 - 2014) from the selected hydrological stations in the Changjiang River basin to analyze hydraulic geometry relationships for streams and rivers in CRN (Figure 1, SI Figure S1 - S2). The hydraulic relationships in this study are comparable with reported data (Harman et al., 2008; Raymond et al., 2012; Jung et al., 2013; Raymond et al., 2013). Table 1 gives the statistics of extracted streams and rivers for the CRN. The calculated whole river/stream surface area is about 18765.7 ± 686.7 km², or 1.07 ± 0.04% of the total land area of the Changjiang River basin.

TABLE 1

Table 1 Statistics for streams and rivers of Changjiang River Basin.

Our analysis demonstrated that 8^th-order rivers, or the main stem rivers, represents the highest proportion of total river surface area at 17.0%, and 7^th-order rivers have the smallest proportion of surface area at 6.8%. Surface areas for smaller 1st to 4th orders of the rivers are similar, accounting for about 12.5% of the total river surface area. The 5^th and 6^th order rivers together comprise the remaining 27.0% of the CRN surface area (SI Figure S5).

Riverine N₂O gas transfer velocity

N₂O transfer velocity (k_N2O) will vary due to differences in the physical controls on water-air turbulence across the different rivers. Our calculated k_N2O presents a longitudinal variation among river order, and tends to decrease with increasing river order (Table 1 and SI Figure S3). The value of k_N2O is the highest for 1^st order rivers, at 6.65 ± 0.11 m d^-1, whereas less than 2.10 ± 0.04 m d^-1 for main stem at 8th order. For the headwaters of 2^nd and 3^rd orders, values of k are 5.96 ± 0.10 and 5.33 ± 0.09 m d^-1, respectively. The values of k for 4^th, 5^th and 6^th order rivers are similar, at 4.14 ± 0.07, 3.98 ± 0.07 and 4.21 ± 0.07 m d^-1, respectively. The average value for k is about 4.61 ± 0.08 m d^-1 for the entire CRN. This value is close to results of recent regional global studies (Alin et al., 2011; Raymond et al., 2013), but is significantly higher than that used in a recent regional calculation for the Changjiang River (Yan et al., 2012; Wang et al., 2015a), and a calculation for the Amazon (Rasera et al., 2013).

River N₂O and its non-linear relationship with DIN

Our in situ measurements of DIN and N₂O covered seasonal, diurnal, and spatial variations in the Changjiang River basin (SI Table S1 and Figure 1). These rivers spanned widely in terms of discharge (ranging from about 6.0 to 28500 m³ s^-1; average: 11770 m³ s^-1), concentrations of ${NO}_{3}^{-}$ (0.163 to 3.09 mg N L^-1; average:1.283 mg N L^-1), ${NH}_{4}^{+}$ (0.04 to 1.619 mg N L^-1; average: 0.271 mg N L^-1), and DIN (0.308 to 3.839 mg N L^-1; average: 1.554 mg N L^-1), and other environmental factors (SI Table S1). Measured N₂O concentrations ranged from 0.104 to 1.11 µg N L^-1 (average 0.521 µg N L^-1), and the associated ΔN₂O concentrations (deviation from saturation) varied between -0.084 and 0.927 µg N L^-1 (average 0.305 µg N L^-1) at all sites.

In-river nitrification and denitrification account for N₂O formation in the water column (Beaulieu et al., 2010) and sediments (Barnes and Owens, 1999; Boyer et al., 2006; Marzadri et al., 2014). Therefore, DIN load could influence N₂O production, but synthesis of their relationship is lacking. Here, we obtained the relationship between riverine [DIN] and [N₂O] using the Michaelis-Menten kinetics equation (Boyer et al., 2006; Turner et al., 2016). We used field observations to calculate the V_max and K_m. Here [N₂O] can be derived from the [DIN] with a V_max of 0.00106 (± 0.000041) mg N L^-1 and a K_m of 1.42 (± 0.11) mg N L^-1 (Figure 2). V_max represents the maximum reaction rate produced by microbial activates. K_m is the affinity constant corresponding to the [DIN] at half V_max. We then used these parameters to regress riverine [N₂O] from [DIN] for the periods of 1970 to 2014 in the entire CRN.

FIGURE 2

Figure 2 Field measurement of N₂O concentrations as a function of DIN concentrations ( n = 55). Red line corresponds to Michaelis–Menten curves. (CH, Chaohu Lake water system, CJ, the Changjiang River, DHS, Datong Hydrological station at the Changjiang River mouth, PYH, Poyang Lake water system).

Temporal N₂O emissions in the CRN

Both compiled and collected [DIN] did not show significant longitudinal variation among different Strahler river orders (SI, Figure S4). Modeled [DIN] by the NEWS-DIN for the Changjiang River, presented significantly increasing temporal trend in the CRN. The average [DIN] was about 0.17 mg N L^-1 in 1970, and below 0.60 mg N L^-1 between 1970 and 1985, and then rapidly increased to about 2.10 mg N L^-1 in 2014. N₂O emission rates ranged from 0.0019 mg m^-2 h^-1 for 8^th river order to 0.006 mg m^-2 h^-1 for headwater streams (1^st and 2^nd river orders) in 1986, and the rate increased to 0.0298 and 0.0944 mg m^-2 h^-1 for 8^th river order and headwater streams, respectively, in 2014. The rate was significantly negatively related to Strahler river order in CRN (Figure 3). The total N₂O emission showed spatial differences for all 1^st to 8^th river orders. N₂O emission for headwater streams (0.234 ± 0.008 and 3.657 ± 0.131 Gg N₂O-N yr^-1 (1 Gg=10⁹g) in 1986 and 2014, respectively) were significantly higher than those for other river orders, and the value for 7^th order was the lowest owing to its less surface area and lower N₂O transfer velocity. N₂O emission for 6^th order is higher than that for 4^th, 5^th, 7^th, and 8^th orders between 1986 and 2014. The total emission was 0.66 Gg N₂O-N yr^-1 in 1986 and increased to 10.3 Gg N₂O-N yr^-1 in 2014 for the entire CRN.

FIGURE 3

Figure 3 The relationship between N₂O emission rate and the Strahler river order in the CRN at each given [DIN] between 1986 and 2014. (DIN = 0.6 mg L^-1 between 1986 and 1990, 1.3 mg L^-1 between 1991 and 1999, 1.6 mg L^-1 between 2000 and 2009, and 2.1 mg L^-1 for 2014).

Mass ratio of N₂O emission to DIN load

The mass ratio of DIN that emits as N₂O against that entering the river network can be considered as the EF (i.e., the ratio of N₂O emission to DIN load). Then, we calculated and characterized the ratio to scale up it within the observed concentration by combining with “spiraling” theory. This ratio can be compared to test the appropriateness of the default IPCC EF_5r (IPCC, 2006) with full considering N loads, biogeochemical and hydrological processes at watershed scale.

Areal rate of total DIN uptake through N₂O emissions (V_f), responds differently to river hydrological processes and [DIN]. We calculated V_f of rivers from 1^st to 8^th orders, and V_f ranged from 0.131 m yr^-1 to 0.436 m yr^-1. V_f was greater in headwater rivers, and declined linearly with river order (Figure 4). In river networks, river size or river hydrological conditions such as river depth and water residence times control the processes of nitrification and denitrification in several ways. DIN can convert efficiently as N₂O in rivers with lower order or small streams (due to their high ratios of riverbed area to water volume) and have a cumulative influence on entire network N₂O emissions because they account for most of the river length (about 74.3% of total river length for 1^st and 2^nd river orders) in the CRN (Table 1). Lower ratios of riverbed area to water volume for rivers with higher order or larger rivers make them less efficient in N₂O production.

FIGURE 4

Figure 4 Relationships between V_f and river order and [DIN] in the CRN. (A) Relationships between V_f and Strahler river order at six given [DIN] (0.6,1.3, 1.6,2.1,4.0, 10.0 mg L^-1); (B) V_f in response to [DIN]. V_f does not increase in parallel with [DIN] when it was higher than 1.6 mg L^-1.

V_f has different dynamic responses to [DIN] in the CRN. Thus, we further estimated the change of V_f values at six given [DIN] (0.6, 1.3, 1.6, 2.1, 4.0, and 10.0 mg L^-1) for each river order (Figure 4B). V_f increased linearly with [DIN] when it was less than 1.6 mg L^-1, but declined linearly as [DIN] continued to increase (Figure 4B). We hypothesized that the linear decline in N₂O V_f was the result of both increased [DIN] and lower efficiency in converting DIN into N₂O. Our Global-NEWS2-DIN model demonstrated that the average [DIN] in the CRN was higher than 1.6 mg L^-1 in 2000, and higher than 2.0 mg L^-1 in 2014. Thus, our study suggests that V_f did not increase in parallel with [DIN] when it was higher than 1.6 mg L^-1 in the CRN after 2000.

The N₂O loss rate (ζ), combining the conversion velocity and river depth, varied from 0.27×10^-4 d^-1 to 37.64×10^-4 d^-1 for the rivers with one to eight orders (Figure 5). The ζ values had significant spatial patterns. The highest ζ value (37.64×10^-4 d⁻¹) was observed in rivers with lower order (1^st to 3^rd) and the lowest (0.27×10^-4 d^-1) in rivers with higher order (8^th). We found that the ζ values declined exponentially with increased river order (R² = 0.99, p< 0.01).

FIGURE 5

Figure 5 Relationships between N₂O loss rate (ζ) and river order.

The mass ratio of DIN that emits as N₂O against that entering the river network, representing EF, has a dynamic pattern in response to both [DIN] and Strahler river orders (Figure 6A). EF was the highest for 6^th river order (average: 0.0556), and the lowest for 7^th order rivers (average: 0.0155) in the CRN during 1986-2014 when [DIN] was in the range of 0.55-2.10 mg L^-1. EF ranged from 0.0015 to 0.050 with an average of 0.0384 for the main stem of the Changjiang River (8^th order river) between 1986 and 2014. For all the Strahler orders, EF increased linearly with [DIN] when [DIN] was less than 1.6 mg L^-1, but declined linearly with increasing [DIN] when it was higher than 1.6 mg L^-1 (Figure 6A), indicating that EF does not increase in parallel with [DIN]. Thus, rivers became less efficient at conversion of DIN although excess DIN increased uptake rate per area of streambed. The threshold value is about 1.6 mg L^-1 in [DIN] for the CRN.

FIGURE 6

Figure 6 Dynamic response of the ratio of DIN emitting as N₂O against that entering the river network to [DIN] for (A) Strahler river orders and (B) the entire CRN.

Our study shows that the integrated EF for the entire CRN, was also not a constant (Figure 6B), depending on both N loads and hydrological conditions at river ecosystem scale. When [DIN] were lower than 0.5 mg L^-1, EF values were below zero, indicating no N₂O emission. EF value ranged from 0.0113% to 0.382%, significantly increased with [DIN], when [DIN] was in the range of 0.55-1.6 mg L^-1. EF value, however, declined linearly with increasing [DIN] when it was higher than 1.6 mg L^-1. The dynamic response of EF to DIN demonstrates the diminishing capacity of river ecosystems to convert excess DIN into N₂O when N load increased as a direct result of human activities.

Discussion

Uncertainties in N₂O emissions

Uncertainties in quantifying N₂O emissions from the CRN existed in two key variables: (1) water surface area of the river network, (2) N₂O transfer velocity.

Generally, the riverine surface area determines the magnitude of riverine N₂O emissions. First, our calculation of the riverine surface area is based on data extracted from a 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM data set in Arcgis 10.2 in combining with the Strahler river order method (Strahler, 1957).

In our study, the lower land of the CRN, such as Taihu Lake area and its adjacent area were excluded due to the limit of DEM (Figure 1). This region represents less than 5% area of the entire Changjiang River basin, and the ignorance of river surface area of this area could not result in significant difference. Thus, despite the uncertainty of DEM extraction result, its influence on water surface area were expected to be insignificant. Secondly, there are some kinds of ephemeral and intermittent rivers in the CRN. These rivers are dry up for certain time of one year and will not involve in air-water N₂O exchange. As there are no direct observation data in the CRN, we used the regression results of USGS data (Raymond et al., 2013) to determine the ephemeral surface area. As a result, the total calculated area ranged from 150.2 to 232.9 km² for the period between 1970 and 2014, accounting for ~1.03 ± 0.10% of the total water surface area. Therefore, the influence of the ephemeral and intermittent rivers to the total emissions could be negligible. Our calculated surface area of the CRN (18765.7 ± 686.7 km²) is much lower than the reported data used by Yan et al. (Yan et al., 2012), because our calculation does not include lakes, reservoirs and ponds within the basin. We could not find any comparable previous regional estimate for river surface water area in the CRN. Nevertheless, our estimated water surface area accounts for 1.07 ± 0.04% of the Changjiang River basin. The ratio is comparable, but a little higher than that estimated by Raymond (Raymond et al., 2013) with 0.83%. Therefore, our estimation of the water surface area for the CRN is considered to be defensible.

Uncertainty in quantifying of gas transfer (k₆₀₀) is significant, especially in the large river network. But considering the difficulties in measuring k₆₀₀ spatially and temporally cross a wide range of the network, empirical models can be adopted to estimate the value (Raymond and Cole, 2001; Striegl et al., 2012; Borges et al., 2015). Factors used to model this value include river channel slope, flow velocity, wind, etc. At one hand, river channel slope and flow velocity, as the major parameters to determine river energy dissipation and turbulence, are intrinsically connected with k₆₀₀ (Raymond et al., 2012). On the other hand, wind is considered to be the key factor controlling on k₆₀₀, particular in open waters such as lakes, estuaries. In this study, we adopted in-river flow velocity and channel slope to calculate k₆₀₀ that appears to vary across river orders. Our results are comparable to that measured by chambers measurements (Striegl et al., 2012), and fall in the range of model-estimated for the rivers of US (Butman and Raymond, 2011), Africa (Borges et al., 2015) and the world (Raymond et al., 2013), respectively, but higher than that from wind-based calculation (Borges et al., 2004; Yan et al., 2012; Rasera et al., 2013; Wang et al., 2015a), especially for the rivers with lower river orders. The inverse relationship between k₆₀₀ and river orders indicates the less gas exchange velocity for large river and more exchange for lower order rivers. Lower order rivers and middle-size rivers usually have a higher vegetation cover and shading (Bouwman et al., 2013b), river surfaces tend to be less affected by wind than larger rivers. This may explain gas transfer velocity in small and middle size rivers are more controlled by channel slope and flow velocity rather than wind speed. Therefore, it results in greater uncertainties in estimating k₆₀₀ based on wind. The surface area weighted average k₆₀₀ value for the entire CRN was 5.22 ± 0.04 m d^-1, much higher than those calculated with wind-based models (Yan et al., 2012; Wang et al., 2015a). Our revisiting calculation of k₆₀₀ values based on in-river flow velocity and channel slope, together with the estimation of the riverine surface area, can also contribute to understand regional riverine budgets of other greenhouse gases such as CO₂ and CH₄ in the CRN.

Previous estimations were assumed that N₂O emissions were linearly correlated with DIN loads, but riverine N₂O production was nonlinearly response to [DIN]. Therefore, it may overestimate N₂O emissions at low DIN loads such as rivers in Africa (Borges et al., 2015) and the Quebec (Soued et al., 2016), and underestimate N₂O emissions at high DIN loads in rivers such as in corn belt (Turner et al., 2015).

N₂O hotspots along the aquatic continuum from headwater streams to the estuary

The N₂O emission rates increased rapidly from 1986 to 2014, and reached to be about 2.332 ± 0.186 mg N m^-2 d^-1 during the period of 2010 - 2014 in headwater streams. The 1^st to 3^rd Strahler order rivers were defined as headwater streams (Gong et al., 2021). We used the surface areas of all 1^st to 3^rd order streams and rivers, and calculated the emission amount to be 1.998 ± 0.195 Gg N₂O-N yr^-1 in 2010-2014. The averaged N₂O emission rate in the mainstem was about 0.686 ± 0.054 mg N m^-2 d^-1 in 2010-2014. We used the surface areas of mainstem, and calculated the emission amount to be 0.786 ± 0.061 Gg N₂O-N yr^-1.

We compiled the data of N₂O emission rate and amount in the Changjiang river estuary reported by Wang et al. (2014b) compared with these for the rivers of the continuum (Figure 7). We found the N₂O emission rate and amount were higher in the headwater streams than these in the mainstem and the estuary. These suggested that headwater streams are hotspots of N₂O emission across the whole aquatic continuum.

FIGURE 7

Figure 7 (A) N₂O emission rate and (B) N₂O emission amoun across the aquatic continuum among the headwater streams (1^st to 3^rd order), mainstem, and the estuary in the same research period of 2010-2014. The data of N₂O emission in the estuary from Wang et al. (2014b).

Significance of the CRN N₂O emissions for the regional budget

As we discussed in the above section, N₂O emissions ranged from 0.66 to 10.3 Gg N₂O-N yr^-1 between 1986 and 2014 for the entire CRN. The estimations of the total N₂O emissions for the period between 1986 and 2002 are much higher than that of the previous report by Yan et al. (2012) for the same period owing to the revisiting higher values for both k₆₀₀ and EF (Yan et al., 2012; Wang et al., 2015a). This is because previous studies calculated N₂O emission from the entire CRN using wind speed method only in main stem, which neglected longitudinal variations of k600 and EF among different river orders. Lower order and middle size rivers usually with higher vegetation cover and shading, of which water surface is recognized having less affected by wind than larger rivers (Bouwman et al., 2013a), but more affected by river channel slope (Raymond and Cole, 2001; Raymond et al., 2013). As a result, mean k_N2O in 1st to 3rd order streams were about 3.55, 3.14, and 2.77 times that of 8^th order rivers (Figure S3b). We also compare our result with the both global and other regional scales in riverine N₂O emission budgets. Taking the N₂O emission for the period between 2010 and 2014, the regional budget in N₂O emissions from the CRN (averaged: 10.3 Gg N₂O-N yr^-1) accounts for about 0.82% - 5.31% of global riverine N₂O budget (Table 2) (Seitzinger and Kroeze, 1998; Seitzinger et al., 2000; Kroeze et al., 2005; Garnier et al., 2009; Beaulieu et al., 2011; Baron et al., 2012; Soued et al., 2016). Our results were lower than that estimated for US Rivers (Baron et al., 2012), but almost equal to the emission from Asian rivers (Hu et al., 2016). Our estimate may represent the minimum limit of N₂O emissions because our estimate only considers rivers and streams, not including other parts of the total inland water areas, such as lakes, reservoirs, and ponds in the Changjiang River basin. In recent decades, the construct of countless dams or reservoirs may alter the hydrological and biogeochemical processes, which have increased the travel time of water, reduced down-river suspended sediment content, and changed (promoted/weakened) denitrification, etc (Wisser et al., 2010). The collection of disturbance types, their legacy and influence on N₂O budgets is beyond this study, but deserve further study to accurately assessment N₂O emissions from the CRN.

TABLE 2

Table 2 N₂O emissions from regional and world rivers.

Conclusions

Our study demonstrates the significant nonlinear relation between DIN and N₂O in the CRN, and models evaluated in this study indicate that the CRN is supersaturated with N₂O since 1986. We find the N₂O emission rate is negatively related to Strahler river orders, and positively related to N loading. Our analysis showed that EF was dynamic rather than static in response to [DIN], it increases with N loading within a certain range, and then declines linearly with increasing N loading.

Our analysis of the N₂O budget in the CRN revealed a significant temporal increasing trend at the regional scale. On average, 0.66 Gg N₂O-N yr^-1 in 1986 and 10.3 Gg N₂O-N N yr^-1 in 2014 were emitted from the entire CRN, respectively.

A comparison of our regional-scale study with other studies from different scales showed a relatively wide range of variability in riverine N₂O emission. Our estimate of N₂O emission in the CRN during 2010-2014 is about 0.82% - 5.31% of global riverine N₂O budget. Our study shows that the headwater streams are hotspots of N₂O emission across the whole aquatic continuum.

The importance of the entire river network to N₂O emissions provides information on the origins of large river network N₂O budget, including artificial reservoirs, natural lakes and wetlands. Further study on the surface area and N₂O of these inland aquatic systems is needed to further evaluate the relative importance of different sources.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author contributions

FW, QY, and WY designed the research and wrote the manuscript. QY carried out the data analysis. ST, PZ, and JW provided the data support. All authors contributed to the article and approved the submitted version.

Funding

This work is jointly supported by National Natural Science Foundation of China (41877483, 42077304, and 42106214), National Key Research and Development Program of China (2019YFC0507805), the Strategic Leading Science and Technology Program of the Chinese Academy of Sciences (XDA20020202).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.1025912/full#supplementary-material

References

Alin S. R., de Fátima F. L., Rasera M., Salimon C. I., Richey J. E., Holtgrieve G. W., et al. (2011). Physical controls on carbon dioxide transfer velocity and flux in low-gradient river systems and implications for regional carbon budgets. J. Geophys. Res.: Biogeosci. 116 (G1), G01009. doi: 10.1029/2010JG001398