Gas transfer velocities of CO₂ in subtropical monsoonal climate streams and small rivers

CO₂ outgassing from rivers is a critical component for evaluating riverine carbon cycle, but it is poorly quantified largely due to limited measurements and modeling of gas transfer velocity in subtropical streams and rivers. We measured CO₂ flux rates and calculated k and partial pressure (pCO₂) in 60 river networks of the Three Gorges Reservoir (TGR) region, a typical area in the upper Yangtze River with monsoonal climate and mountainous terrain. The determined k₆₀₀ (gas transfer velocity normalized to a Schmidt number of 600 (k₆₀₀) at a temperature of 20 ^∘C) value (48.4±53.2 cm h⁻¹) showed large variability due to spatial variations in physical processes related to surface water turbulence. Our flux-derived k values using chambers were comparable with k values using the model derived from flow velocities based on a subset of data. Unlike in open waters, e.g., lakes, k₆₀₀ is more pertinent to flow velocity and water depth in the studied river systems. Our results show that TGR river networks emitted approx. 0.69 to 0.71 Tg CO₂ (1 Tg =10¹² g) during the monsoon period using varying approaches such as chambers, derived k₆₀₀ values and models. This study suggests that incorporating scale-appropriate k measurements into extensive pCO₂ investigations is required to refine basin-wide carbon budgets in subtropical streams and small rivers. We concluded that the simple parameterization of k₆₀₀ as a function of morphological characteristics is site specific for regions and watersheds and hence highly variable in rivers of the upper Yangtze. k₆₀₀ models should be developed for stream studies to evaluate the contribution of these regions to atmospheric CO₂.

1 Introduction

Rivers serve as a significant contributor of CO₂ to the atmosphere (Raymond et al., 2013; Cole et al., 2007; Li et al., 2012; Tranvik et al., 2009). As a consequence, accurate quantification of riverine CO₂ emissions is a key component to estimate net continental carbon (C) flux (Raymond et al., 2013). More detailed observational data and accurate measurement techniques are critical to refining riverine C budgets (Li and Bush, 2015; Raymond and Cole, 2001). Generally, two methods are used to estimate CO₂ areal fluxes from the river system: direct measurements using floating chambers (FCs) and the indirect calculation of a thin boundary layer (TBL) model, which is dependent on the gas concentration gradient at the air–water interface and gas transfer velocity, k (Guerin et al., 2007; Xiao et al., 2014). Direct measurements are normally laborious, while the latter method is simple and thus preferred (Butman and Raymond, 2011; Lauerwald et al., 2015; Li et al., 2012, 2013; Ran et al., 2015).

The areal flux of CO₂ (F, mmol m⁻² d⁻¹) via the water–air interface by TBL is described as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }F=k\times K_h\times p{\mathrm{CO}}_{\mathrm{2}},\text{(2)}& {\displaystyle }& {\displaystyle }{K}_h=\mathrm{1}{\mathrm{0}}^{-(\mathrm{1.11}+\mathrm{0.016}\times T-\mathrm{0.00007}\times T^{\mathrm{2}})},\end{array}$$

where k (m d⁻¹) is the gas transfer velocity of CO₂ (also referred to as piston velocity) at the in situ temperature (Li et al., 2016); pCO₂ (µatm) is the pCO₂ gradient at the air–water interface (Borges et al., 2004). K_h (mmol m⁻³ µatm⁻¹) is the aqueous-phase solubility coefficient of CO₂ corrected using in situ temperature (T in ^∘C) (Li et al., 2016).

pCO₂ can be measured well in various aquatic systems; however, the accuracy of the estimation of flux is dependent on the k value. Broad ranges of k for CO₂ (Raymond and Cole, 2001; Raymond et al., 2012; Borges et al., 2004) were reported due to variations in techniques, tracers used and governing processes. k is controlled by turbulence at the surface aqueous boundary layer; hence, k₆₀₀ (the standardized gas transfer velocity at a temperature of 20 ^∘C is valid for freshwater) is parameterized as a function of wind speed in open water systems of reservoirs, lakes and oceans (Borges et al., 2004; Guerin et al., 2007; Wanninkhof et al., 2009). While in streams and small rivers, turbulence at the water–air interface is generated by shear stresses at the streambed, and thus k is modeled using channel slope, water depth and water velocity in particular (Raymond et al., 2012; Alin et al., 2011). Variable formulations of k have been established by numerous theoretical, laboratory and field studies; nonetheless, a better constraint on k levels is still required as its levels are very significant and specific due to large heterogeneity in the hydrodynamics and physical characteristics of river networks. This highlights the importance of k measurements in a wide range of environments for accurate upscaling of CO₂ evasion and for parameterizing the physical controls on k₆₀₀. However, only a few studies provide information on k for riverine CO₂ flux in Asia (Alin et al., 2011; Ran et al., 2015), and those studies do not address the variability of k in China's small rivers and streams.

Limited studies have demonstrated higher levels of k in large Chinese rivers (Liu et al., 2017; Ran et al., 2015, 2017; Alin et al., 2011), which contributed to much higher CO₂ areal flux, particularly in China's monsoonal rivers that are impacted by hydrological seasonality. The monsoonal flow pattern and thus flow velocity is expected to be different than other rivers in the world; as a consequence, k levels should be different than others and potentially higher in subtropical monsoonal rivers.

Figure 1 Map of sampling locations for major rivers and streams in the Three Gorges Reservoir region, China.

Considerable efforts, such as the use of purposeful (Crusius and Wanninkhof, 2003; Jean-Baptiste and Poisson, 2000) and natural tracers (Wanninkhof, 1992) and FCs (Alin et al., 2011; Borges et al., 2004; Prytherch et al., 2017; Guerin et al., 2007), have been carried out to estimate accurate k values. The direct determination of k by FCs is more popular due to the simplicity of the technique for short-term CO₂ flux measurements (Prytherch et al., 2017; Raymond and Cole, 2001; Xiao et al., 2014). Prior reports, however, have demonstrated that k values and the parameterization of k as a function of wind and/or flow velocity (probably water depth) vary widely across rivers and streams (Raymond and Cole, 2001; Raymond et al., 2012). To contribute to this debate, extensive investigation was first accomplished for the determination of k in rivers and streams of the upper Yangtze using the FC method. Models of k were further developed using hydraulic properties (i.e., flow velocity, water depth) by flux measurements with chambers and a TBL model. Our recent study preliminarily investigated pCO₂ and air–water CO₂ areal flux as well as their controls from fluvial networks in the Three Gorges Reservoir (TGR) area (Li et al., 2018). The past study was based on two fieldwork studies, and diffusive models from other rivers and/or regions were used. Here, we derive k levels and develop the gas transfer model in this area (mountainous streams and small rivers) for more accurate quantification of CO₂ areal flux; it also serves for the fluvial networks in the Yangtze River or others with similar hydrology and geomorphology. Moreover, we conducted detailed field campaigns in two contrasting rivers, the Daning and Qijiang, for models (Fig. 1) and the rest were TGR streams and small rivers. The study thus clearly shows distinct differences from a previous study (Li et al., 2018) through the contribution of specific new objectives and data supplements with wider significance. Our new contributions to the literature thus include (1) the determination and controls of k levels for small rivers and streams in subtropical areas of China and (2) new models developed in subtropical mountainous river networks. The outcome of this study is expected to help in the accurate estimation of CO₂ evasion from subtropical rivers and streams and thus refine the riverine C budget over a regional and/or basin scale.

2 Materials and methods

2.1 Study areas

All field measurements were carried out in the rivers and streams of the Three Gorges Reservoir (TGR) region (28^∘44^′–31^∘40^′ N, 106^∘10^′–111^∘10^′ E) located in the upper Yangtze River, China (Fig. 1). This region is subject to a humid subtropical monsoon climate with an average annual temperature ranging between 15 and 19 ^∘C. Average annual precipitation is approx. 1250 mm with large intra-annual and interannual variability. About 75 % of the total annual rainfall is concentrated between April and September (Li et al., 2018).

The river sub-catchments include large-scale river networks covering the majority of the tributaries of the Yangtze in the TGR region, i.e., data for 48 tributaries were collected. These tributaries have drainage areas that vary widely from 100 to 4400 km² with width ranging from 1 m to less than 100 m. The annual discharges from these tributaries have a broad spectrum of 1.8–112 m³ s⁻¹. Detailed samplings were conducted in the two largest rivers of Daning (35 sampling sites) and Qijiang (32 sites) in the TGR region. These two river basins drain catchment areas of 4200 and 4400 km². The studied river systems had widths < 100 m, and we thus defined them as small rivers and streams. The Daning and Qijiang river systems are underlain by widely carbonate rock and located in a typical karst area. The location of the sampling sites is shown in Fig. 1. Detailed information on the sampling sites and primary data is presented in the Supplement (Table S5). The sampling sites are outside the Reservoirs and are not affected by dam operation.

2.2 Water sampling and analyses

Three fieldwork campaigns from the main river networks in the TGR region were undertaken during May through August in 2016 (i.e., 18–22 May for Daning, 21 June–2 July for all tributaries of TGR and 15–18 August for Qijiang). A total of 115 discrete grab samples were collected (each sample consisted of three replicates). Running waters were taken using pre-acid-washed 5 L high-density polyethylene (HDPE) plastic containers from depths of 10 cm below the surface. The samples were filtered through prebaked Whatman GF/F (0.7 µm pore size) filters on the sampling day and immediately stored in acid-washed HDPE bottles. The bottles were transported in an icebox to the laboratory and stored at 4 ^∘C for analysis. Concentrations of dissolved organic carbon (DOC) were determined within 7 days of water collection (Mao et al., 2017).

Water temperature (T), pH, DO saturation (DO %) and electrical conductivity (EC) were measured in situ by calibrated multiparameter sondes (HQ40d HACH, USA, and YSI 6600, YSI incorporated, USA). pH, the key parameter for pCO₂ calculation, was measured to a precision of ±0.01, and the pH sonde was calibrated by certified reference materials (CRMs) before measurements with an accuracy of better than ±0.2 %. Atmospheric CO₂ concentrations were determined in situ using an EGM-4 (environmental gas monitor; PP SYSTEMS Corporation, USA). Total alkalinity was measured using a fixed endpoint titration method with 0.0200 mol L⁻¹ hydrochloric acid (HCl) on the sampling day. DOC concentration was measured using a total organic carbon analyzer (TOC-5000, Shimadzu, Japan) with a precision better than 3 % (Mao et al., 2017). All the solvents and reagents used in the experiments were of analytical reagent grade.

Concomitant stream width, depth and flow velocity were determined along the cross section, and flow velocity was determined using a portable flow meter LS300-A (China); the meter shows an error of < 1.5 %. Wind speed at 1 m over the water surface (U₁) and air temperature (T_a) were measured with a Testo 410-1 handheld anemometer (Germany). Wind speed at 10 m of height (U₁₀, unit in m s⁻¹) was calculated using the following formula (Crusius and Wanninkhof, 2003):

$$\begin{array}{}\text{(3)}& U_{\mathrm{10}}=U_Z\left[\mathrm{1}+{{\displaystyle \frac{\left(C_{\mathrm{d}\mathrm{10}}\right)}{K}}}^{\mathrm{1}/\mathrm{2}}\times \ln \left({\displaystyle \frac{\mathrm{10}}{z}}\right)\right],\end{array}$$

where C_d10 is the drag coefficient at 10 m of height (0.0013 m s⁻¹), K is the von Karman constant (0.41) and z is the height (m) of wind speed measurement. $U_{\mathrm{10}}=\mathrm{1.208}\times U_{\mathrm{1}}$ as we measured the wind speed at a height of 1 m (U₁).

Aqueous pCO₂ was computed from the measurements of pH, total alkalinity and water temperature using a CO₂ system (k₁ and k₂ are from Millero, 1979) (Lewis et al., 1998). This program can yield high-quality data (Li et al., 2012, 2013; Borges et al., 2004).

2.3 Water-to-air CO₂ fluxes using the FC method

FCs (30 cm in diameter, 30 cm in height) were deployed to measure air–water CO₂ fluxes and transfer velocities. They were made of cylindrical polyvinyl chloride (PVC) pipe with a volume of 21.20 L and a surface area of 0.071 m². These nontransparent, thermally insulated vertical tubes covered by aluminum foil were connected via CO₂-impermeable rubber–polymer tubing (with outer and inner diameters of 0.5 and 0.35 cm, respectively) to a portable nondispersive infrared CO₂ analyzer EGM-4 (PPSystems). Air was circulated through the EGM-4 instrument via an air filter using an integral pump at a flow rate of 350 mL min⁻¹. The chamber method was widely used and more details of the advantages and limits of chambers were reviewed elsewhere (Alin et al., 2011; Borges et al., 2004; Xiao et al., 2014).

Chamber measurements were conducted by deploying two replicate chambers or one chamber two times at each site. In sampling sites with low and favorable flow conditions (Supplement Fig. S1), freely drifting chambers (DCs) were deployed, while sites in rivers and streams with higher flow velocity were sampled with anchored chambers (ACs) (Ran et al., 2017). DCs were used at sampling sites with a current velocity of < 0.1 m s⁻¹; this resulted in a limited number of sites (a total of six sites) using DCs. ACs would create an overestimation of CO₂ emissions by a factor of several-fold (i.e., > 2) in our study region (Lorke et al., 2015). Data were logged automatically and continuously at 1 min intervals over a given span of time (normally 5–10 min) after enclosure. The CO₂ area flux (mg m⁻² h⁻¹) was calculated using the following formula:

$$\begin{array}{}\text{(4)}& F=\mathrm{60}\times {\displaystyle \frac{\mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}\times M\times P\times T_{\mathrm{0}}}{\mathrm{d}t\times V_{\mathrm{0}}\times P_{\mathrm{0}}\times T}}H,\end{array}$$

where dpCO₂∕dt is the rate of concentration change in FCs (µL L⁻¹ min⁻¹); M is the molar mass of CO₂ (g mol⁻¹); P is the atmosphere pressure of the sampling site (Pa); T is the chamber absolute temperature of the sampling time (K); V₀ is the molar volume (22.4 L mol⁻¹); P₀ is atmosphere pressure (101 325 Pa); T₀ is absolute temperature (273.15 K) under the standard condition; and H is the chamber height above the water surface (m) (Alin et al., 2011). We accepted the flux data that had a good linear regression of flux against time (R²≥0.95, p < 0.01) following the manufacturer specifications. In our sampling points, all measured fluxes were retained since the floating chambers yielded linearly increasing CO₂ against time.

Water samples from a total of 115 sites were collected. Floating chambers with replicates were deployed in 101 sites (32 sampling sites in Daning, 37 sites in TGR river networks and 32 sites in Qijiang). The sampling period covered the spring and summer season, and our sampling points are reasonable considering a water area of 433 km². For example, 16 sites were collected for the Yangtze system to examine hydrological and geomorphological controls on pCO₂ (Liu et al., 2017), with 17 sites for dynamic biogeochemical controls on riverine pCO₂ in the Yangtze basin (Liu et al., 2016). Similar to other studies, sampling and flux measurements in the day would tend to underestimate the CO₂ evasion rate (Bodmer et al., 2016).

2.4 Calculations of the gas transfer velocity

The k was calculated by reorganizing Eq. (1). To make comparisons, k is normalized to a Schmidt (Sc) number of 600 (k₆₀₀) at a temperature of 20 ^∘C.

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{600}}=kT{\left({\displaystyle \frac{\mathrm{600}}{\mathit{Sc}T}}\right)}^{-\mathrm{0.5}}\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{Sc}T=\mathrm{1911.1}-\mathrm{118.11}T+\mathrm{3.4527}{T}^{\mathrm{2}}-\mathrm{0.04132}{T}^{\mathrm{3}}\end{array}$$

k_T is the measured value at the in situ temperature (T, unit in ^∘C), and SC_T is the Schmidt number of temperature T. A dependency of −0.5 was employed here as measurements were made in turbulent rivers and streams in this study (Alin et al., 2011; Borges et al., 2004; Wanninkhof, 1992).

2.5 Estimation of river water area

The water surface is an important parameter for CO₂ efflux estimation, as it depends on its climate, channel geometry and topography. River water area therefore largely fluctuates with a much higher areal extent of water surface, particularly in the monsoonal season. However, most studies do not consider this change, and a fraction of the drainage area is used in river water area calculation (Zhang et al., 2017). In our study, 90 m resolution SRTM DEM (Shuttle Radar Topography Mission digital elevation model) data and Landsat images in the dry season were used to delineate the river network and thus the water area (Zhang et al., 2018), while stream orders were not extracted. The water area of river systems is generally much higher in the monsoonal season in comparison to the dry season; for instance, the Yellow River showed a 1.4-fold higher water area in the wet season than in the dry season (Ran et al., 2015). The available dry-season image was likely to underestimate CO₂ estimation.

2.6 Data processing

Prior to statistical analysis, we excluded k₆₀₀ data for samples with an air–water pCO₂ gradient < 110 µatm, since the error in the k₆₀₀ calculations are drastically enhanced when pCO₂ approaches zero (Borges et al., 2004; Alin et al., 2011), and datasets with pCO₂ > 110 µatm provide an error of < 10 % on k₆₀₀ computation. Thus, we discarded the samples (36.7 % of sampling points with flux measurements) with pCO₂ < 110 µatm for k₆₀₀ model development, while all samples were included for flux estimations from the diffusive TBL model and floating chambers.

Spatial differences (Daning, Qijiang and all tributaries of TGR region) were tested using the nonparametric Mann–Whitney U test. Multivariate statistics, such as correlation and stepwise multiple linear regression, were performed for the models of k₆₀₀ using potential physical parameters of wind speed, water depth and current velocity as independent variables (Alin et al., 2011). Data analyses were conducted from both separated data and combined data for the river systems. k models were obtained by water depth using data from the TGR rivers and by flow velocity in the Qijiang, while models were not developed for Daning and the combined data. All statistical relationships were significant at p < 0.05. The statistical processes were conducted using SigmaPlot 11.0 and SPSS 16.0 for Windows (Li et al., 2009, 2016).

Table 1 Statistics of all the data from the three river systems (for separated statistics please refer to Figs. S2 and S3).

CI: confidence interval.

3 Results

3.1 CO₂ partial pressure and key water quality variables

Significant spatial variations in water temperature, pH, pCO₂ and DOC were observed among the Daning, TGR and Qijiang rivers, whereas alkalinity did not display such variability (Fig. S2). pH varied from 7.47 to 8.76 with the exception of two quite high values of 9.38 and 8.87 (total mean: 8.39±0.29). Significantly lower pH was observed in TGR rivers (8.21±0.33) (Table 1; p < 0.001; Fig. S2). pCO₂ varied between 50 and 4830 µatm with a mean of 846±819 µatm (Table 1). There were 28.7 % of samples that had pCO₂ levels lower than 410 µatm, while the studied rivers were overall supersaturated with reference to atmospheric CO₂ and act as a source of atmospheric CO₂. The pCO₂ levels were 2.1 to 2.6-fold higher in TGR rivers than Daning (483±294 µatm) and Qijiang River (614±316 µatm) (Fig. S2).

There was a significantly higher concentration of DOC in the TGR rivers (12.83±7.16 mg L⁻¹) than Daning and Qijiang River (3.76±5.79 vs. 1.07±0.33 mg L⁻¹ in Qijiang and Daning) (p < 0.001; Fig. S3). Moreover, Qijiang showed a significantly higher concentration of DOC than Daning (3.76±5.79 vs. 1.07±0.33 mg L⁻¹ in Qijiang and Daning) (p < 0.001 by Mann–Whitney rank sum test; Fig. S3).

3.2 CO₂ flux using floating chambers

The calculated CO₂ areal fluxes were higher in TGR rivers (217.7±334.7 mmol m⁻² d⁻¹, n=35), followed by Daning (122.0±239.4 mmol m⁻² d⁻¹, n=28) and Qijiang rivers (50.3±177.2 mmol m⁻² d⁻¹, n=32) (Fig. 2). The higher CO₂ evasion from the TGR rivers is consistent with high riverine pCO₂ levels. The mean CO₂ emission rate was 133.1±269.1 mmol m⁻² d⁻¹ (n=95) in all three rivers sampled. The mean CO₂ flux differed significantly between TGR rivers and Qijiang (Fig. 2).

Table 2 Comparison of different models for CO₂ areal flux estimation using combined data (unit is mmol m⁻² d⁻¹ for CO₂ areal flux and cm h⁻¹ for k₆₀₀).

^aCI: confidence interval. ^b The flow-velocity-based model is from a subset of the data (please refer to Fig. 4) ^c Mean value determined using floating chambers (FCs). ^d This figure is revised to be 49.6 cm h⁻¹ if the model (k₆₀₀=62.879 FV +6.8357, R²=0.52, p=0.019) is used (the model is obtained by taking out two extreme values; please refer to Fig. 4c), and the corresponding CO₂ areal flux is 203±190 mmol m⁻² d⁻¹.

Figure 2 Box plots of CO₂ emission rates from floating chambers in the three investigated river systems (different letters represent statistical differences at p < 0.05 by Mann–Whitney rank sum test). The black and red lines, lower and upper edges, and bars and dots in or outside the boxes demonstrate median and mean values, 25th and 75th, 5th and 95th, and < 5th and > 95th percentiles of all data, respectively. For an interpretation of the references to color in this figure legend, the reader is referred to the web version of this article. “Total” means combined data from the three river systems.

3.3 k levels

A total of 64 data points were used (10 for Daning River, 33 for TGR rivers and 21 for Qijiang River) to develop the k model after removal of samples with pCO₂ less than 110 µatm (Table 2). No significant variability in k₆₀₀ values was observed among the three rivers sampled (Fig. 3). The mean k₆₀₀ was relatively higher in Qijiang (60.2±78.9 cm h⁻¹), followed by Daning (50.2±20.1 cm h⁻¹) and TGR rivers (40.4±37.6 cm h⁻¹), while the median k₆₀₀ was higher in Daning (50.5 cm h⁻¹), followed by TGR rivers (30.0 cm h⁻¹) and Qijiang (25.8 cm h⁻¹) (Fig. 3; Supplement Table S1). Combined k₆₀₀ data were averaged to 48.4±53.2 cm h⁻¹ (95 % CI: 35.1–61.7), and this is 1.5-fold higher than the median value (32.2 cm h⁻¹) (Fig. 3).

Contrary to our expectations, no significant relationship was observed between k₆₀₀, water depth and current velocity using all data in the three river systems (TGR streams and small rivers, Daning and Qjiang) (Fig. S4). There were no statistically significant relationships between k₆₀₀ and wind speed using separated data or combined data. Flow velocity showed a slightly linear relation with k₆₀₀, and an extreme high value of k₆₀₀ was observed during periods of higher flow velocity (Fig. S4a) using combined data. A similar trend was also observed between water depth and k₆₀₀ values (Fig. S4b). k₆₀₀ as a function of water depth was obtained in the TGR rivers, but it explained only 30 % of the variance in k₆₀₀. However, the model using data from Qijiang could explain 68 % of the variance in k₆₀₀ (Fig. 4b), and it was in line with general theory.

Figure 3 Box plots of k₆₀₀ levels in the three investigated river systems (there is not a statistically significant difference in k among sites by Mann–Whitney rank sum test). The black and red lines, lower and upper edges, and bars and dots in or outside the boxes demonstrate median and mean values, 25th and 75th, 5th and 95th, and < 5th and > 95th percentiles of all data, respectively. For an interpretation of the references to color in this figure legend, the reader is referred to the web version of this article. “Total” means combined data from the three river systems.

4 Discussion

4.1 Uncertainty assessment of pCO₂ and flux-derived k₆₀₀ values

The uncertainty of flux-derived k values mainly stems from ΔpCO₂ and flux measurements (Bodmer et al., 2016; Golub et al., 2017; Lorke et al., 2015). Thus, we provided uncertainty assessments for dominant sources of uncertainty from measurements of aquatic pCO₂ and CO₂ areal flux since the uncertainty of atmospheric CO₂ measurement could be neglected.

In our study, aquatic pCO₂ was computed based on pH, alkalinity and water temperature rather than directly measured. Recent studies highlighted pCO₂ uncertainty caused by systematic errors over empiric random errors (Golub et al., 2017). Systematic errors are mainly attributed to instrument limitations, i.e., sondes of pH and water temperature. The relative accuracy of the temperature meters was ±0.1 ^∘C according to the manufacturer specifications, and thus the uncertainty of water T propagated on uncertainty in pCO₂ was minor (Golub et al., 2017). Systematic errors therefore stem from pH, which has been proved to be a key parameter for biased pCO₂ estimation calculated from the aquatic carbon system (Li et al., 2013; Abril et al., 2015). We used a high-accuracy pH electrode and the pH meters were carefully calibrated using CRMs; in situ measurements showed an uncertainty of ±0.01. We then run an uncertainty of ±0.01 pH to quantify the pCO₂ uncertainty, and an uncertainty of ±3 % was observed. Systematic errors thus seemed to show little effect on pCO₂ errors in our study.

Figure 4 The k₆₀₀ as a function of water depth (WD) using data from TGR rivers (a), flow velocity (FV) using data from Qijiang (b) and comparison of the developed model with other models (c) (others without significant relationships between k and physical factors are not shown). The solid lines show regression, the dashed lines represent the 95 % confidence band, and the red dash-dotted line represents the model developed by Alin et al. (2011). Extreme values of 260 and 274 cm h⁻¹ are removed in panel (b); the revised model would be k₆₀₀=62.879 FV +6.8357, R²=0.52, p=0.019 (in panel c, 1 – the revised model, 2 – model from Ran et al., 2017, 3 – model from Ran et al., 2015, 4 – model from Alin et al., 2011, 5 – model from Liu et al., 2017; 1 – k₆₀₀=62.879 FV +6.8357; 2 – k₆₀₀=58.47 FV +7.99; 3 – k₆₀₀=13.677 exp (1.1 FV); 4 – k₆₀₀=35 FV +13.82; 5 – $k_{\mathrm{600}}=\mathrm{6.5}{\mathrm{FV}}^{\mathrm{2}}+\mathrm{12.9}$ FV +0.3). The unit of k in models 1–4 is cm h⁻¹, and the unit m d⁻¹ for model 5 is transferred to cm h⁻¹.

Random errors are from the repeatability of carbonate measurements. Two replicates for each sample showed an uncertainty within ±5 %, indicating that uncertainty in pCO₂ calculation from alkalinity measurements could be minor.

The measured pH ranges also exhibited great effects on pCO₂ uncertainty (Hunt et al., 2011; Abril et al., 2015). At low pH, pCO₂ can be overestimated when calculated from pH and alkalinity (Abril et al., 2015). Samples for CO₂ fluxes estimated from pH and alkalinity showed pH average of 8.39±0.29 (median 8.46 with quartiles of 8.24–8.56) (n=115). Thus, the overestimation of calculated CO₂ areal flux from pH and alkalinity is likely to be minor. Further, the contribution of organic matter to non-carbonate alkalinity is likely to be neglected because of low DOC (mean 6.67 mg L⁻¹; median 2.51 mg L⁻¹) (Hunt et al., 2011; Li et al., 2013).

Efforts have been devoted to measurement techniques (comparison of FC, eddy covariance (EC) and boundary layer model; BLM) for improving CO₂ quantification from rivers because of a notable contribution of inland waters to the global C budget, which could have a large effect on the magnitude of the terrestrial C sink. Prior studies have reported inconsistent trends of CO₂ area flux with these methods. For instance, CO₂ areal flux from FC was much lower than EC (Podgrajsek et al., 2014), while areal flux from FC was higher than both EC and BLM elsewhere (Erkkila et al., 2018); however, Schilder et al. (2013) demonstrated that areal flux from BLM was 33 %–320 % of in situ FC measurements. Despite unresolved errors in various techniques and additional perturbations from FC, the FC method is currently a simple and preferred technique for CO₂ flux because choosing the right k value remains a major challenge and others require high workloads (Martinsen et al., 2018).

A recent study further reported fundamental differences in CO₂ emission rates between ACs and free DFs (Lorke et al., 2015); i.e., ACs biased the gas areal flux higher by a factor of 2.0–5.5. However, some studies observed that ACs showed reasonable agreement with other flux measurement techniques (Galfalk et al., 2013), and this method is straightforward, inexpensive and relatively simple; hence, it is widely used (Ran et al., 2017). Water–air interface CO₂ flux measurements were primarily made using ACs in our studied streams and small rivers because of relatively high current velocity; otherwise, floating chambers will travel far during the measurement period. In addition, inflatable rings were used for sealing the chamber headspace and the submergence of ACs was minimal; therefore, our measurements were potentially overestimated, but reasonable. We could not test for the overestimation of ACs in this study; studies with modified FCs, i.e., DCs and the integration of ACs and DCs, and multi-method comparison studies including FCs, ECs and BLM should be conducted for a reliable chamber method.

Our model was from a subset of the data (i.e., Qijiang); CO₂ flux from our model was in good agreement with the fluxes from FC, which determined k and other models when the developed model was applied for the whole dataset (please refer to Tables 2 and 3). The comparison of the fluxes from variable methods suggested that the model can be used for riverine CO₂ flux at catchment scale or regional scale, though it cannot be used at individual sites. Recent studies, however, did not test the applicability of models when k₆₀₀ models from other regions were employed. Our k₆₀₀ values were close to the average of Ran et al. (2015) (measured with drifting chambers) and Liu et al. (2017) (measured with static chambers in canoe shape); this indicated that our potential overestimation was limited. However, since we had very limited drifting chamber measurements because of high current velocity, the relationships with chamber-derived k₆₀₀ values and flow velocity–depth only with the drifting chamber data could not be tested. We acknowledge that k₆₀₀ could be overestimated using AFs.

Table 3 CO₂ emission during the monsoonal period (May through October) from all rivers sampled in the study. (a) Upscaling using CO₂ areal flux (mean ±SD) by FC during the monsoonal period. (b) Upscaling using determined k₆₀₀ average and models (whole datasets are used here).

A total water area of approx. 430 km² for all tributaries (water area is from Landsat ETM+ in 2015); CO₂ emission upscaling (Tg CO₂ during May through October) was conducted during the monsoonal period because of the sampling in this period.

The extremely high values (two values of 260 and 274 cm h⁻¹) are outside of the global ranges and also considerably higher than k₆₀₀ values in Asian rivers. Furthermore, the revised model was comparable to published models (Fig. 4), i.e., models of Ran et al. (2015) (measured with drifting chambers) and Liu et al. (2017) (measured with static chambers in canoe shape), which suggested that the exclusion of the two extreme values was reasonable, and this was further supported by the CO₂ flux using different approaches (Tables 2 and 3).

Sampling seasonality considerably regulated riverine pCO₂ and gas transfer velocity and thus water–air interface CO₂ evasion rate (Ran et al., 2015; Li et al., 2012). We sampled waters in the wet season (monsoonal period) because it showed a wider range of flow velocity and thus it covered the k₆₀₀ levels in the whole hydrological season. The wet season generally had a higher current velocity and thus a higher gas transfer velocity (Ran et al., 2015), while aquatic pCO₂ was variable with seasonality. We recently reported that riverine pCO₂ in the wet season was 81 % of the level in the dry season (Li et al., 2018), and a prior study on the Yellow River reported that the k level in the wet season was 1.8-fold higher than in the dry season (Ran et al., 2015). Another study on the Wuding River demonstrated that the k level in the wet season was 83 %–130 % of that in the dry season (Ran et al., 2017). Thus, we acknowledge a certain amount of errors on the annual flux estimation from sampling campaigns during the wet season in the TGR area, but this uncertainty is not significant because the diluted pCO₂ could alleviate the overestimated emission with an increased k level in the wet season (for a more detailed discussion, please refer to SOM).

4.2 Determined k values relative to world rivers

We derived k values for the first time in subtropical streams and small rivers. Our determined k₆₀₀ levels with a 95 % CI of 35.1 to 61.7 (mean: 48.4) cm h⁻¹ compared well with a compilation of data for streams and small rivers (e.g., 3–70 cm h⁻¹) (Raymond et al., 2012). Our determined k₆₀₀ values are greater than the global river average (8–33 cm h⁻¹) (Raymond et al., 2013; Butman and Raymond, 2011) and much higher than the mean for tropical and temperate large rivers (5–31 cm h⁻¹) (Alin et al., 2011). These studies evidence the fact that k₆₀₀ values are highly variable in streams and small rivers (Alin et al., 2011; Ran et al., 2015). Though the mean k₆₀₀ in the TGR, Daning and Qijiang is higher than the global mean, it is consistent with k₆₀₀ values in the main stream and river networks of the turbulent Yellow River (42±17 cm h⁻¹) (Ran et al., 2015) and Yangtze (38±40 cm h⁻¹) (Liu et al., 2017) (Table S2).

The calculated pCO₂ levels were within the published range, but towards the lower end of published concentrations compiled elsewhere (Cole and Caraco, 2001; Li et al., 2013). The total mean pCO₂ (846±819 µatm) in the TGR, Daning and Qijiang rivers was one-third lower than the global river average (3220 µatm) (Cole and Caraco, 2001). The pCO₂ lower than most of the world's river systems, particularly the under-saturated values, demonstrated that heterotrophic respiration of terrestrially derived DOC was not significant. Compared with high alkalinity, the limited delivery of DOC, particularly in the Daning and Qijiang river systems (Figs. S2 and S3), also indicated that in-stream respiration was limited. These two river systems are characterized by karst terrain and underlain by carbonate rock, for which photosynthetic uptake of dissolved CO₂ and carbonate mineral dissolution considerably regulated aquatic pCO₂ (Zhang et al., 2017).

Higher pH levels were observed in the Daning and Qijiang river systems (p < 0.05 by Mann–Whitney rank sum test), where more carbonate rock exists that is characterized by karst terrain. Our pH range was comparable to the recent study on karst rivers in China (Zhang et al., 2017). Quite high values (8.39±0.29, ranging between 7.47 and 9.38; 95 % confidence interval: 8.33–8.44) could increase the importance of the chemical enhancement; nonetheless, few studies have taken chemical enhancement into account (Wanninkhof and Knox, 1996; Alshboul and Lorke, 2015). Chemical enhancement can increase the CO₂ areal flux by a factor of several-fold in lentic systems with low gas transfer velocity, with the enhancement factor decreasing quickly as k₆₀₀ increased (Alshboul and Lorke, 2015). Our studied rivers are located in mountainous area with high k₆₀₀, which could cause minor chemical enhancement. This chemical enhancement of CO₂ flux was also reported to be limited in high-pH and turbulent rivers (Zhang et al., 2017).

4.3 Hydraulic controls of k₆₀₀

It has been well established that k₆₀₀ is governed by a multitude of physical factors, particularly current velocity, wind speed, stream slope and water depth, among which wind speed is the dominant factor of k in open waters such as large rivers and estuaries (Alin et al., 2011; Borges et al., 2004; Crusius and Wanninkhof, 2003; Raymond and Cole, 2001). In contrast, k₆₀₀ in small rivers and streams is closely linked to flow velocity, water depth and channel slope (Alin et al., 2011; Raymond et al., 2012). Several studies reported that the combined contribution of flow velocity and wind speed to k is significant in large rivers (Beaulieu et al., 2012; Ran et al., 2015). Thus, k₆₀₀ values are higher in the Yellow River (ca. 0–120 cm h⁻¹) compared to the low-gradient Mekong River (0–60 cm h⁻¹) (Alin et al., 2011; Ran et al., 2015) due to higher flow velocity in the Yellow River (1.8 m s⁻¹) than the Mekong River (0.9±0.4 m s⁻¹), resulting in greater surface turbulence and a higher k₆₀₀ level in the Yellow (42±17 cm h⁻¹) than Mekong River (15±9 cm h⁻¹). This could substantiate the higher k₆₀₀ levels and spatial changes in k₆₀₀ values of our three river systems. For instance, similar to other turbulent rivers in China (Ran et al., 2015, 2017), high k₆₀₀ values in the TGR, Daning and Qjiang rivers were due to mountainous terrain catchment, high current velocity (10–150 cm s⁻¹) (Fig. 4b), bottom roughness and shallow water depth (10–150 cm) (Fig. 4a). It has been suggested that shallow water enhances bottom shear, and the resultant turbulence increases k values (Alin et al., 2011; Raymond et al., 2012). These physical controls are highly variable across environmental types (Fig. 4a and b); hence, k values are expected to vary widely (Fig. 3). The k₆₀₀ values in the TGR rivers showed a wider range (1–177 cm h⁻¹; Fig. 3; Table S1), spanning more than 2 orders of magnitude across the region, and this is consistent with the considerable variability in physical processes related to water turbulence across environmental settings. A similarly broad range of k₆₀₀ levels was also observed in China's Yellow basin (ca. 0–123 cm h⁻¹) (Ran et al., 2015, 2017).

Insignificant relationships between riverine k₆₀₀ and wind speed were consistent with earlier studies (Alin et al., 2011; Raymond et al., 2012). The lack of strong correlation between k₆₀₀ and physical factors using the combined data was probably due to the combined effect of both flow velocity and water depth, as well as large diversity of channel morphology, both across and within river networks in the entire catchment (60 000 km²). This is further collaborated by weak correlations between k₆₀₀ and flow velocity in the TGR rivers (Fig. 4), where one or two samples were taken for a large-scale examination. We provided new insights into k₆₀₀ parameterized using current velocity. Nonetheless, k₆₀₀ from our flow-velocity-based model (Fig. 4b) was potentially largely overestimated with consideration of other measurements (Alin et al., 2011; Ran et al., 2015, 2017). When several extreme values were removed, k₆₀₀ (cm h⁻¹) was parameterized as follows: k₆₀₀=62.879 FV +6.8357, R²=0.52, p=0.019, flow velocity with a unit of m s⁻¹. This revised model was in good agreement with the model in the river networks of the Yellow River (Ran et al., 2017), but much lower than the model developed in the Yangtze system (Liu et al., 2017) (Fig. 4c). This was reasonable because k₆₀₀ values in the Yangtze system were from large rivers with higher turbulence than the Yellow and our studied rivers. Furthermore, the determined k₆₀₀ using FCs was, on average, consistent with the revised model (Table 2). These differences in the relationship between spatial changes in k₆₀₀ values and physical characteristics further corroborated the heterogeneity of channel geomorphology and hydraulic conditions across the investigated rivers.

Subtropical streams and small rivers are biologically more active and recognized to exert higher CO₂ areal flux to the atmosphere; however, their contribution to riverine carbon cycling is still poorly quantified because of data paucity and the absence of k in particular. Larger uncertainty of riverine CO₂ emission in China was anticipated by the use of k₆₀₀ from other continents or climate zones. For instance, k₆₀₀ for CO₂ emission from tributaries in the Yellow River and karst rivers originated from the model in the Mekong (Zhang et al., 2017), Pearl (Yao et al., 2007), Longchuan (Li et al., 2012) and Metropolitan rivers (Wang et al., 2017), which are mostly from temperate regions. Our k₆₀₀ values will therefore largely improve the estimation of CO₂ evasion from subtropical streams and small rivers and improve the refinement of the riverine carbon budget. More studies, however, are clearly needed to build the model based on flow velocity and slope–water depth given the difficulty in k quantification on a large scale.

4.4 Implications for large-scale estimation

We compared CO₂ areal flux from FCs, the models developed here (Fig. 4) and other studies (Alin et al., 2011) (Tables 2 and 3). CO₂ evasion was estimated for rivers in China with k values ranging between 8 and 15 cm h⁻¹ (Li et al., 2012; Yao et al., 2007; Wang et al., 2011) (Table S2). These estimates of CO₂ evasion rate were considerably lower than using present k₆₀₀ values (48.4±53.2 cm h⁻¹). For instance, CO₂ emission rates in the Longchuan River (e.g., k=8 cm h⁻¹) and Pearl River tributaries (e.g., k=8–15 cm h⁻¹) were 3 to 6 times higher using present k values compared to earlier estimates. We found that the determined k₆₀₀ average was marginally beyond the levels from the water-depth-based model and the model developed by Alin et al. (2011), while it was equivalent to the flow-velocity-based revised model, resulting in similar patterns of CO₂ emission rates (Table 2). Hence, the selection of k values would significantly hamper the accuracy of the flux estimation. Therefore, k must be estimated along with pCO₂ measurements to make accurate flux estimations.

We used our measured CO₂ emission rate from FCs for upscaling flux estimates during the monsoonal period given the sampling in this period and it was found to be 0.70 Tg CO₂ (1 Tg =10¹² g) for all rivers sampled in our study (Table 3a). The estimated emission in the monsoonal period was close to that of the revised model (0.71±0.66 (95 % confidence interval: 0.46–0.94) Tg CO₂) and using the determined k average, i.e., 0.69±0.65 (95 % confidence interval: 0.45–0.93) Tg CO₂, but slightly higher than the estimation using the water-depth-based model (0.54±0.51 Tg CO₂) and Alin's model (0.53±0.50 Tg CO₂) (Table 3b). This comparable CO₂ flux further substantiated the exclusion of extreme k₆₀₀ values for developing the model (Fig. 4). The CO₂ evasion comparison for variable approaches also implied that the original flow-velocity-based model (two extreme k₆₀₀ values were included; Fig. 4b) largely overestimated the CO₂ fluxes, i.e., 1.66±1.55 (1.08–2.23) Tg CO₂ was 2.3–3-fold higher than other estimations (Table 3b) and our earlier evasion using TBL on the TGR river networks (Li et al., 2018). Moreover, our estimated CO₂ emission during the monsoonal period also suggests that CO₂ annual emissions from rivers and streams in this area were previously underestimated, i.e., 0.03 Tg CO₂ yr⁻¹ (Li et al., 2017) and 0.37–0.44 Tg CO₂ yr⁻¹ (Yang et al., 2013), as the former used a TBL model with a lower k level; the latter employed floating chambers, but they both sampled very limited tributaries (i.e., two to three rivers). Therefore, measurements of k must be made mandatory along with pCO₂ measurement in river and stream studies.

5 Conclusions

We provided the first determination of gas transfer velocity (k) in subtropical streams and small rivers in the upper Yangtze. High variability in k values (mean 48.4±53.2 cm h⁻¹) was observed, reflecting the variability of the morphological characteristics of water turbulence both within and across river networks. We highlighted the fact that k estimates from empirical models should be pursued with caution and the significance of incorporating k measurements along with extensive pCO₂ investigation is highly essential for upscaling to watershed- and/or regional-scale carbon (C) budget.

Riverine pCO₂ and CO₂ areal flux showed pronounced spatial variability with much higher levels in the TGR rivers. The CO₂ areal flux was averaged at 133.1±269.1 mmol m⁻² d⁻¹ using FCs, and the resulting emission during the monsoonal period was around 0.7 Tg CO₂, similar to the scaling-up emission with the determined k and the revised flow-velocity-based model, while also marginally above the water-depth-based model. More work is clearly needed to refine k modeling in the river systems of the upper Yangtze River for evaluating regional C budgets.