Background

Groundwater is a critical water resource around the globe ensuring food and water security [1]. Irrigation with groundwater for agricultural activities is a common practice in many arid and semi-arid regions [2]. However, over exploitation of groundwater has led to severe groundwater depletion in several regions of the world [3]. Areas most affected by groundwater depletion are California and Midwest in the US, Northern India, and the North China Plain (NCP) [4]. In the NCP approximately 70% of the irrigated area is currently groundwater-fed [5, 6] which causes the groundwater level to drop by more than 2 cm per year [7]. While several studies focus on the importance and influence of the dramatic groundwater level drop, the further fate of the pumped irrigation water is much less studied [8].

Groundwater is typically 10- to 100-fold over-saturated with CO2 [9], even up to 250 times the atmospheric equilibrium value [10]. Thus, if it is pumped to the surface, CO2 is released. First estimates show that CO2 emissions from groundwater pumping probably represent a globally significant source of CO2 [11], with both CO2 liberation from the water and CO2 production due to pump energy generation contributing to the negative climate impact of groundwater irrigation [12]. CO2 has been recognized as the dominant driver of climate change resulting in tremendous deteriorative impacts on the environment and society at a global scale [13]. As one of the major emission sources of CO2, agricultural land plays an important role in the global carbon cycle [14]. Direct CO2 emissions from agriculture have been well-documented in many regions [15,16,17]. However, indirect emission sources, such as irrigation drainages and rivers receiving drainage water, may account for a large part of the uncertainties in the carbon budgets of agricultural ecosystems [18,19,20,21].

If CO2 containing water equilibrates with the atmosphere, not only the dissolved CO2 but also other species of the carbonate system have to be considered [22]. Existing estimates of CO2 emissions from groundwater pumping assume complete evaporation of the pumped water, leaving solutes and solids at the surface. Under such conditions, each molecule of CO2 outgassing is theoretically produced from two molecule of carbonate hardness (refers to HCO3− + CO32−) to re-equilibrate with atmospheric CO2 accompanied by one molecule of carbonate precipitation simultaneously [11, 23]. Thus, maximum CO2 liberation from groundwater dissolved inorganic carbon could be calculated from the bicarbonate content of the groundwater [9]. When a large part of dissolved inorganic carbon (DIC) containing irrigation water drains into surface waters, DIC may have four different fates: evasion as CO2, carbonate precipitation, fixation into biomass by TOC production, and downstream transport [24].

Streams are recognized as important CO2 sources to the atmosphere [25, 26]. In streams the distribution among the four fates of groundwater inorganic carbon is controlled by both hydrological and biochemical processes and their complex interactions [10, 27]. Flow velocity, characterized by fluctuating hydraulic conditions due to flood irrigation or excessive utilization of water resources, usually results in enhanced carbon export from terrestrial and wetland habitats to fluvial networks that subsequently become a source of CO2 [25, 28, 29]. Turbulence at the air–water interface, which is dependent on the river geomorphology and the river flow, affects the gas transfer coefficient of CO2, which ultimately determines fluvial CO2 emissions [30,31,32]. On the other hand, increasing nutrient levels stimulate biological processes such as algal growth, which decrease DIC concentration by converting CO2 to organic carbon via photosynthesis [33]. Agricultural drainage water is typically high in inorganic nutrients. Thus, our hypothesis is that both nutrients and flow velocity mediate the relative importance of uptake versus outgassing of CO2 accompanied with calcite precipitation from the dense agricultural drainage channels transporting pumping groundwater.

However, quantifying these processes of DIC depletion remains challenging. It is difficult to disentangle the effects of both factors on the fate of inorganic carbon from field observations only, primarily due to uncontrolled environmental conditions. We, therefore, performed a mesocosm experiment in the field to explore how CO2 emissions from irrigation groundwater respond to nutrients and flow velocity treatments. The goals of this study are to (1) quantitatively testify the effects of nutrients and flow velocity on CO2 emissions and the carbon budget of irrigation drainage channels in a typical agricultural region in China, and (2) provide implications to CO2 emission potential from the irrigation groundwater in the NCP as well as other similar regions worldwide.

Methods

Study area

Our study was conducted in the Yucheng Comprehensive Experiment Station of Dezhou Irrigation District (DID) between 115°45′–117°36′ E and 36°24′–38°00′ N, located in the North China Plain. The annual-averaged precipitation is 587 mm (ranges from 286 to 1034 mm). Precipitation occurs mostly from June to September, which accounts for 75% of the total annual amount. Annual evaporation is between 900 and 1400 mm and annual mean temperature is 12.8 °C [34]. The groundwater for agricultural irrigation is characterized by high alkalinity and enriched with agricultural nutrients when it enters the drainage channels [35]. Groundwater bicarbonate concentrations in the NCP can range between 1.55 and 7.65 mmol L−1 [36]. The ditches and canals are mainly artificial and lined with concrete.

Historically, DID has been characterized by anthropogenic activities, particularly intensive agriculture development, which has been causing a high demand of water resources in the NCP [37]. Irrigation drainages were built during the twentieth century to increase crop yield in the salinized soil by diverting water from the lower reaches of the Yellow River to DID. Given that the surface water cannot fulfill the demand from agricultural and industrial production, over-exploitation of groundwater had been frequently observed through the sharp decrease of the groundwater table and heavy pollution therein [38]. To eliminate the influence of soil salinity on crop growth, flooding irrigation is still commonly applied among local farmers. This enhances nutrient inputs from irrigated groundwater into the surface water network that subsequently become a source of CO2. Based on farmer’s practice when they irrigate the crops, our experiment simulate that stream drainage collected from groundwater-irrigated agricultural land. For better representation of local groundwater, we used groundwater from a local irrigation well as inlet water to our experimental system.

Mesocosm setup

Cubic mesocosms (volume 0.2 m3) were constructed to mimic drainage channels. The physical design of our mesocosms follows that from Petersen et al. [39]. The systems are composed of sixteen cuboid polyvinyl chloride plastic tanks with 0.9 m in length, 0.3 m in width, and 0.6 m in depth. Depth of water was set to 0.5 m. Mesocosms were placed outdoor and remained open and in contact with the field atmosphere to ensure natural light and weather conditions. The containers were sunk into the ground to buffer heating from ambient atmosphere temperature and solar radiation. The electrical wiring for power supplying was also buried. The upper edge of the tanks was 10 cm above the ground so that no surface runoff could flow into containers during rainstorms. The containers allow a factorial design combining two nutrient levels with two flow velocity levels in four replicates (Additional file 1: Figure S1). The tanks were equipped with two alternative submerged pumps to mimic river flowing. Each mesocosm had an individual pump (65 W or 15 W) to provide continuous flowing condition and to preclude water pressure differences of connected pipes when sharing the same pump under the same treatment. The pumps were placed in mesh cages to prevent snails, insects, and floating fragments from suction into the pump and to reduce the risk of clogging. The inlet of the pump was on the shorter side of the tank to simulate a unit of water channel (Fig. 1). We set the length to width ratio as 3:1 to establish a realistic flow regime in the central part of our setting. The middle section along the longitudinal direction of these tanks can be considered a homogeneous flow field with simple physical boundaries.

Fig. 1
figure 1

Conceptual diagram of a single container with a pump. The up-right insert panel shows a conceptual model of DIC transformation in the mesocosm. The fate of DIC was categorized into four pathways: CO2 evasion, TOC production, calcite precipitation, and downstream transport

Full size image

Experimental procedure

Groundwater was collected from a nearby agricultural irrigation well from a depth of 5 m and pumped into the pre-cleaned tanks without sediments. Nutrients were added into the tanks at the beginning of the experiment for the high-nutrient treatment, while low-nutrient tanks received no additional nutrients. Particularly, high-nutrient tanks were instantaneously fed with 40 mg P m−3 and 4000 mg N m−3 as K2HPO4 and Ca(NO3)2, respectively. We set the nutrient enrichment level by mimicking the two alternative states existing in the local drainage system. The nitrate and ammonia concentrations of the study groundwater well were approximately 0.02 mg N L−1 and 5.8 μg N L−1, respectively, while the nitrate concentration of drainages in a field survey in Yucheng (on 16th and 17th of September) varied from 0.02 to 0.42 mg N L−1 with a mean value of 0.17 mg N L−1 (Table 1). Other parameters in the experimental water were comparable to that in the surface water as well (Table 1). We designed the initial nitrate concentration for the high nutrient experimental group at a level above the average nitrate concentration (0.3 mg N L−1). Phosphorous was added with the N:P molar ratio of 10:1 following Lone et al. [40].

Table 1 Characteristics of the groundwater and surface water in the study area
Full size table

We choose flow velocities in surrounding drainage channels as a reference. Flow velocity of local drainages varied from 0.001 to 0.340 m s−1 during our field sampling. Accordingly, two levels of flow velocity were set to 0.1 and 0.4 m s−1. Pumps were cleaned manually every 3 days. It is expected that there was little CO2 degassing through pumping itself, since the submerged pumps were tested for air-tightness before usage and it is expected that any artifact would be equal between low and high nutrient treatments, because the same pumps were used.

Overall, there were four (2 × 2) treatments with 4 replicates each in our experimental design, which are denoted as: high nutrient and high flow velocity (A + H); high nutrient and low flow velocity (A + L); low nutrient and high flow velocity (N + H); and low nutrient and low flow velocity (N + L).

The experiments was conducted between 31st August and 15th September 2018. Sampling and measurements of the physical and chemical parameters was performed at 0, 4, 16, 32, 64, 112, 160, 280 h. The experiment measurement was finished after 12 days when the CO2 concentration had reached equilibrium with the atmosphere. Due to the mesocosm malfunction, the N + H treatment has only one replicate after 50 h.

Sampling and analysis

Water temperature, dissolved oxygen (DO), and electrical conductivity (EC) were measured using multi-parameter probes (Hach H40d, USA). Water level and flow velocity were recorded using a flow meter (HR-2, China) in parallel to each sampling. Water samples were collected to quantify concentrations of nitrate, ammonium, DIC, alkalinity, and TOC. For determination of NO3−–N and NH4+–N concentration, water samples were filtered using 0.45 μm filters by spectrophotometry (TU-1810DSPC, China). Calibration curves were produced using reference samples according to quality control standards and were then applied to evaluate data from each set of samples. Reagents, procedural blanks, and samples were measured twice in parallel, with average values reported. The differences between two replicates were within 5% of the mean value for all samples. NO3− and NH4+ concentrations could be measured at precisions of 2.6% and 8.6%, respectively. DIC and TOC were analyzed with a Vario TOC Analyzer (Elementar Analysensysteme GmbH, Germany). The detection limits of DIC and TOC were 4 and 1 μg L−1 with precisions of 1.5% and 5%, respectively. Alkalinity was determined by titration with ~ 0.01 M H2SO4 at a precision of 6%. While the measured TOC data could not be directly used. We used two methods to quantify CO2 emissions during our experiments: a carbon budget approach and direct flux measurements using floating chambers (see Additional file 1 for method details).

Calculation of C budgets

Changes of DIC in the mesocosms can be separated into three parts: carbon dioxide evasion, calcite precipitation, and total organic carbon (TOC) production (Fig. 1):

$$ \Delta {\text{DIC}} = \Delta {\text{CO}}_{2} + \Delta {\text{CaCO}}_{3} + \Delta {\text{TOC}}{.} $$
(1)

To verify whether carbonate precipitation was likely to occur within the mesocosm, we calculated calcite saturation indices (SIs) using PHREEQC [41]. Because we did not measure calcium concentrations in the mesocosm, we used historical rainy season data from irrigation ditches in Dezhou District [42] and derived possible SIs with 1.28 ± 0.16 (mean value ± standard deviation). Therefore, we assumed calcite precipitation occurred in all mesocosms, with 1 mol each of CO2 outgassed and CaCO3 precipitated from 2 mol of HCO3− according to Wood and Hyndman [11]:

$$ 2{\text{{HCO}}_{3}}^{ - } + {\text{Ca}}^{2 + } = {\text{CO}}_{2} \uparrow + {\text{CaCO}}_{3} \downarrow + {\text{H}}_{2} {\text{O}}{.} $$
(2)

CO2 evasion depends on CO2 concentration and the physical gas transfer velocity [32]. The CO2 concentration is affected by the thermodynamic carbonate equilibrium. TOC formation is a biological process that converts inorganic carbon into organic matter. At any time, the variation for DIC is denoted as:

$$ {\text{DIC}}_{0} - {\text{DIC}}_{t} = \Delta {\text{CO}}_{2} + \Delta {\text{CaCO}}_{3} + \left( {{\text{TOC}}_{t} - {\text{TOC}}_{0} } \right). $$
(3)

Equation (3) is based on our hypothesis that CO2 in the system was over-saturated and TOC was generated during the experiment, \({\text{DIC}}_{0}\) and \({\text{TOC}}_{0}\) is the amount of DIC and TOC at the initial time, respectively, and \({\text{DIC}}_{t}\) and \({\text{TOC}}_{t}\) is the amount of DIC and TOC at time t, respectively. \(\Delta {\text{CO}}_{2}\) and \(\Delta {\text{CaCO}}_{3}\) is the amount of CO2 outgassing and calcite precipitating during time 0 to time t, respectively. Combined with Eq. (2),

$$ \Delta {\text{CO}}_{2} = \Delta {\text{CaCO}}_{3} . $$
(4)

This is a conservative estimate as we note that the ratio of \(\Delta {\text{CO}}_{2} :\Delta {\text{CaCO}}_{3}\) can range from 1.18 to 1.88 [9]. The CO2 evasion could be deduced by the difference between the variation amount of DIC and the variation amount of TOC, re-arranging Eqs. (3) with (4) gives:

$$ \Delta {\text{CO}}_{2} = \frac{{\left( {{\text{DIC}}_{0} - {\text{DIC}}_{t} } \right) - \left( {{\text{TOC}}_{t} - {\text{TOC}}_{0} } \right)}}{2}. $$
(5)

Because TOC data were not available, we estimated the TOC generation indirectly from total inorganic nitrogen (TIN) consumption. Assuming the molar C:N ratio in freshwater being higher than the Redfield ratio of 106:16 [43, 44], TOC production could be estimated from the consumption of nitrate and ammonium, which are the main species of nitrogen in the groundwater system with nitrate accounting for the majority, based on a C:N ratio of 10:1 [45]. We assume that denitrification as a process that removes N from the system can be neglected in our case, because our system was turbulent and well aerated.

Our mass balance was calculated by multiplying the concentrations with the varying water volume in the setup (mesocosm bottom area × decreasing water level with negligible volume of the water in pipes and pump), so the evaporation would not influence our result.

We performed a series of analysis of variance (ANOVA) to explore the patterns of DIC, TOC, and CO2 and their relationship to nutrient and flow velocity. The concentration of CO2 at given DIC and alkalinity was calculated using standard dissociation constants of the carbonate system [46]. Data were log-transformed to satisfy assumptions of residuals normal distribution. ANOVA was applied to test for overall statistical differences among treatments. The ratio to the initial amount of the above parameters for each of the mesocosm was carried out to test for the treatment effects using ANOVA. In the model, nutrient and flow velocity were treated as the fixed effect, and temporal pseudo-replication from repeated sampling over time was considered as nested within each mesocosm as random effects. The comparison of DIC and CO2 changes is given by the different slope of a first-order reaction of log10 DIC or log10 CO2 over time to test for statistical differences for each treatment and mesocosm. All calculations, statistical analysis, and data visualization were performed using R version 3.5.1 [47].

Results

Significant differences among four treatments

The differences of flow velocities and nutrients among the four treatments were significant at the beginning of the experiment (Table 2). The nitrate concentration was similar between high and low flow with a larger variation range in high flow treatments. Because of some clogging of the pumps, the flow velocity fluctuated at a larger range at high flow compared to the low flow velocity.

Table 2 Flow velocities and chemical parameters at start (t = 0 h) and the end (t = 280 h) of the experiment among four treatments
Full size table

At first, DO saturation increased exponentially over time in all mesocosms (Fig. 2). The initial DO concentration at high flow was higher than that at low flow at the starting time, probably because more turbulences at high flow mixed more oxygen into the water before the first sampling. DO saturation varied dissimilarly over time and treatments as well. Mesocosms with nutrient additions developed a clear over-saturation, showing that TOC production was stimulated by nutrient additions. The A + L treatment reached the highest over-saturation. On the contrary, the control mesocosms rapidly reached equilibrium with atmosphere. As expected, equilibration was faster at high flow.

Fig. 2
figure 2

Variations of DO saturation (a–d) and EC (e–h) during the experiment. The black dashed lines show the exponential fits and colored areas indicate its 95% confidence intervals. The red dashed lines show 100% of the DO saturation

Full size image

EC in the high flow velocity treatments showed a similar trend with a sharp decrease from 2.24 to 1.95 mS cm−1 in the first 150 h after which it stayed constant (Fig. 2). The low flow velocity treatments showed a steadily decrease from 2.24 to 1.95 mS cm−1 from the start to the end of experiment.

DIC and alkalinity were used to calculate the CO2 concentrations at the beginning and in the end. At the end of the experiment in the high nutrient treatment, pCO2 in water was below that in the atmosphere, indicating CO2 uptake. Assuming alkalinity has a linear decrease over time accompanied with known DIC changes (Fig. 3a–d), we calculated the pCO2 at any time in each mesocosm using the CO2SYS program [48]. Thus, the time when pCO2 in each mesocosm was equal to pCO2 in the atmosphere (pCO2air = 395 μatm) could be estimated. After about 125 and 235 h, the systems switched to CO2 uptake under the A + H and A + L treatments, respectively (Additional file 1: Figure S2). The N + L treatment reached CO2 equilibration with atmosphere until the end of experiment.

Fig. 3
figure 3

DIC depletion (a–d) and TOC generation (e–h) in each container over time among four treatments, the dashed lines in a–d show the exponential fits (p < 0.001), and the dashed lines in e–h show the linear fits with p values and colored areas indicate 95% confidence intervals

Full size image

DIC and TOC

The DIC amount in all treatments decreased continuously by 31–48% during the experimental period (Fig. 3). Decrease of DIC was largest in the A + H treatment and lowest in the N + L treatment. The amount of DIC in the system was significantly different between the two nutrient treatments (p = 0.041), as well as between the two flow velocities (p = 0.020). DIC in the N + H treatment decreased rapidly at the beginning of the experiment (particularly at the first 50 h). The decreasing rate of DIC under the A + H treatment surpassed that under the N + H treatment, while the rate was relatively slow but steadily after 170 h. The rate of DIC decrease was also strongly related to flow velocity (p = 0.006) but had insignificant relationship with nutrients. Moreover, there was a significant interaction of flow velocity and nutrients in the DIC decrease slope showing a combined effect between the two factors (p = 0.043).

The TOC increase in the high nutrient treatments was evidently much higher than the one in the treatment without nutrients addition (Fig. 3) and independent of flow velocity. In the low nutrient treatments, TOC did not significantly change during the experiment, indicating little TOC production in these systems. This is consistent with the DO data which showed no O2 oversaturation in the low nutrient treatments (Fig. 2).

CO 2 evasions among four treatments

All mesocosms were net sources of CO2 during the experimental period (Fig. 4). Always positive fluxes (except at the end of the experiment) were confirmed by our floating chamber measurements (Fig. 4e–h). The high flow treatments reached CO2 equilibrium between water and atmosphere faster (Additional file 1: Figure S2). The rate of CO2 change also shows the significant difference between flow velocity treatments (p = 0.024). High flow velocity sped up CO2 evasion. More CO2 was emitted significantly under the high flow velocity (p = 0.005) and low nutrient condition (p = 0.096). CO2 emission from the low nutrient treatments was 35% higher than that in the high nutrient treatments in average.

Fig. 4
figure 4

Cumulative CO2 emissions in each container over time among four treatments calculated from DIC (a–d) and from floating chamber measurements (e–h), the dashed lines show the exponential fits (p < 0.001) and colored areas indicate 95% confidence level, the red dots indicate the negative CO2 fluxes, namely there was CO2 uptake occurring

Full size image

CO2 evasions of the N + H treatment increased rapidly during the first 80 h, and remained nearly constant after 150 h. The A + H treatment emitted markedly a lot of CO2 at the beginning but the flux became lower over time. CO2 evasion from the two low flow treatments exhibited a similar behavior with relatively gentle and steady CO2 emissions, and the N + L treatment released evidently more CO2 than the A + L treatment.

Combined with the result of pCO2 changes (Additional file 1: Figure S2), CO2 uptake could occur by the time when pCO2 concentration dropped below pCO2air. The floating chamber measurement also showed that some CO2 was taken up (CO2 flux < 0) under the high nutrient treatments at the time of 280 h. The rate of CO2 uptake can be calculated from pCO2 and the gas transfer velocity (k600):

$$ f_{{{\text{CO}}_{2} }} = k_{600} \cdot k_{{\text{H}}} \cdot \left( {p{\text{CO}}_{2} - p{\text{CO}}_{{2{\text{air}}}} } \right). $$
(6)

\(k_{{\text{H}}}\) is Henry’s constant for CO2 at a given temperature and salinity [49]. Gas transfer velocity (\(k_{600}\)) could be calculated from the initial DIC depletion in the low nutrient treatments assuming that there was little TOC produced in these systems (Fig. 3g, h). Assuming k600 being constant and equivalent under the same flow velocity treatments, k600 was 0.93 and 0.23 m day−1 in the high and low flow velocity treatments, respectively. Using these k600 values and the final pCO2 we can calculate the CO2 flux at the end of the experiment using Eq. (6). Multiplying that flux with the duration of the under-saturated period results in an average total CO2 uptake of 42 and 9.5 μmol in the A + H and A + L treatments, respectively. This is much lower than the observed changes in DIC. Thus, uptake did not significantly contribute to the total C budget of the experiment (Fig. 4).

Discussion

Carbon dynamics

Our results show that flow velocity and nutrient level determine the rate and amount of CO2 evasion, respectively. CO2 evasion depends on both CO2 concentration and the gas transport coefficient (k) in Eq. (6) [50, 51]. With the CO2 evasions under the same flow treatments showing similar tendency over time, flow velocity altered k, which determined the rate of CO2 emission. Meanwhile, treatments with the same nutrient level ended with similar amount of CO2 evasion, suggesting that the nutrient level ultimately dictates the quantity of CO2 emission. Providing that more TOC was generated in the high nutrient treatment, our results confirm that TOC generation reallocate the depletion of inorganic carbon.

Our results also highlight the interaction of the flow and nutrient speeding up the depletion rate of DIC in drainage channels (p = 0.043). However, there was no significant synergistic effect of velocity and nutrients on the final DIC amount remaining in the system (p = 0.40). Flow velocity significantly altered the slope of DIC and CO2 over time, indicating that flow velocity dominated the rate of DIC consumption and CO2 emissions. Simultaneously, flow velocity influences the quantity of DIC consumption as well, possibly by means of regulating on ecosystem metabolism. High flow can intensify the capacity of an ecosystem to store and process carbon [52]. TOC production at high flow velocity could be enhanced by higher light and nutrient availability throughout the water column due to better mixing [53]. Low flow on the other hand may promote sedimentation, thereby improving water clarity which may increase TOC production [33]. Nutrients affect the quantity of both CO2 and DIC by altering the pattern of the inorganic carbon utilization through ecosystem processes (i.e., gross primary production, ecosystem respiration) [54, 55]. Nutrient enrichment may accelerate TOC processing, transformation, and export, potentially altering food-web dynamics and ecosystem stability in the long term [54]. Higher water temperature enhances photosynthesis enzyme activity that favors the biomass growth [56]. In our study, the synergy of both factors cannot be explained by light availability or temperature, because turbidity was low due to lack of sediment and temperature did not differ between treatments (mean water temperature of 25.0 ± 1.0 ℃ within each sampling). Most probable, better mixing at high flow prevented possible local nutrient limitation of photosynthesis in our experiment.

Uncertainty

Effects of C:N ratio on carbon budgets

Our TOC results are sensitive towards the chosen C:N ratio (10:1 in our case). We examined the responsiveness of the C:N ratio to the DIC allocation in our experiment by calculating carbon mass balances for different C:N ratios [C:N = 5, 106:16 (Redfield ratio), 10 (used in this study), and 15] (Additional file 1: Figure S3). Variation of the C:N ratio changed the absolute amount of CO2 emitted from the high nutrient treatments but did not change the general pattern.

Alkalinity dynamics

In a closed system changes in alkalinity only occur by processes including either a solid phase or conversion to a gas, which then escapes to the atmosphere [57]. Beside the definition of the alkalinity as acid neutralizing capacity of a solution titrated to the CO2 equivalence point, the alkalinity can be also described by the equivalent sum of conservative cations (those that do not affect alkalinity) minus the sum of conservative anions:

$$ {\text{Alkalinity}} = \left[ {{\text{Na}}^{ + } } \right] + \left[ {{\text{K}}^{ + } } \right] + 2\left[ {{\text{Mg}}^{2 + } } \right] + 2\left[ {{\text{Ca}}^{2 + } } \right] + 2\left[ {{\text{Mn}}^{2 + } } \right] + \left[ {{\text{NH}}_{4}^{ + } } \right] + 2\left[ {{\text{Fe}}^{2 + } } \right] - \left[ {{\text{Cl}}^{ - } } \right] - 2\left[ {{\text{SO}}_{4}^{2 - } } \right] - \left[ {{\text{NO}}_{3}^{ - } } \right]. $$
(7)

The changes of alkalinity caused by nitrate and ammonium consumption can account for 0.1–26.2% of the total alkalinity variation. The continuous EC decline during the experiment (Fig. 2) was also observed indicating that not only nitrogen transformation occurred, but there was also ion precipitation inducing alkalinity dropping down. The main ions in the groundwater were SO42−, HCO3− (affects the DIC change), Na+, and Ca2+ [58]. Chemical theory predicts that upon evaporation, half of the bicarbonate in the groundwater will re-equilibrate with the atmosphere releasing CO2 while the other half precipitates, mostly as calcite [11]. Thus, calcite precipitation must have affected our pCO2 calculations and led to alkalinity decrease through Ca2+ removal from the system. Assuming the measured alkalinity changes (153 μmol L−1 in average) were completely caused by calcite precipitation, then there would be 77 μmol L−1 Ca2+ being utilized for CaCO3 formed in each mesocosm. This is roughly consistent with our CO2 emission calculations, because assuming one calcite formed per emitted CO2 would result in 69–168 μmol L−1 Ca2+ being utilized for CaCO3 formed.

Method uncertainties in CO2 flux measurement

We used two methods to quantify CO2 emissions during our experiments: a carbon budget approach and direct flux measurements using floating chambers. Theoretically, the cumulative CO2 flux measured by the chamber should match the cumulative C loss from the water,

$$ S\mathop \smallint \limits_{0}^{t} f\left( t \right){\text{d}}t = \left( {c_{t} - c_{0} } \right)V, $$
(8)

where \(S\) is the surface area of the water; \(f\left( t \right)\) is the CO2 emission flux as a function of time; V is the volume of the water, \(c_{0}\) and \(c_{t}\) are the CO2 concentrations at time 0 and t, respectively.

We found that the results of the direct flux measurement and the DIC calculation were not consistent. The flux chamber method resulted in a twenty-time overestimation of the CO2 evasion especially at higher flow (Fig. 4). This is likely because of the artificial turbulence generated by the floating chamber [59, 60]. Deployment of floating chamber in the mesocosm system could cause a large overestimation of the CO2 evasion compared with the one in natural waters [61, 62], which would lead to an unrealistic carbon budget with all carbon removed from the system by outgassing in our carbon budget. Thus, floating chamber measurements cannot be used to quantify CO2 emissions in flume experiments.

Implications and upscaling

Drainage channels receive irrigation water from agriculture and route the nutrient containing groundwater directly into streams. High flow velocity accelerates the CO2 evasion to the atmosphere rapidly (Fig. 5), leaving less opportunity for channel metabolism and processing [32]. Previous research also suggests that CO2 concentrations are largely reduced by intensive primary production in rivers [63]. Thus, nutrient additions have the potential to attenuate CO2 emissions. The formed biomass is transported downstream and could be trapped in lentic parts of the drainage network, where particle sedimentation mainly occurs. Notably, collected particulate carbon needs to be properly buried or dredged. Otherwise methane production is likely to take place in the sediment accumulations. That would largely increase the GHG effect of the drainage network due to the larger global warming potential of methane compared to CO2 [21, 30]. Alternatively, in-stream plankton could potentially be mineralized by heterotrophic respiration [19, 64]. In that case, CO2 would only be temporarily bound and later liberated further downstream.

Fig. 5
figure 5

i Fitted curves of temporal changes of the DIC, TOC, and cumulative CO2 emission among the four treatments (a, c exponential regression, b linear regression) and ii the percentage fate of the dissolved inorganic carbon in the mesocosm

Full size image

At a given rate of CO2 emission, the flow velocity determines the required distance until the CO2 concentration in a drainage ditch equilibrates with the atmosphere. In a hypothetical drainage system of infinite length with no biological processes involved, it could be anticipated that the only fate of inorganic carbon is to outgas into atmosphere. Our study suggest that up to 20% of DIC would be converted to CO2 and 58% of the initial DIC would be transported downstream considering no other processes were involved until the pCO2 in the water equilibrates with the atmosphere. Previous study reported large variations of the release rate of atmospheric carbon from rivers to the atmosphere (2–30%) [65]. From field measurements of GHG emissions, Ran et al. [66] concluded that 35% of the carbon exported into the Yellow River network was degassed from the entire watershed during fluvial transport, which is comparable with our estimation.

Irrigation in the study area mainly occurs in March–April (to ensure the water requirement of winter-wheat during the growth period), and October (before wheat sowing). In the case of extreme drought in summer, irrigation is also needed for maize growth in July or August. During other periods the water residence time in the agricultural drainage is quite long, because runoff in irrigated agriculture is mainly related to the intensity of irrigation and rainfall events in comparison to the infiltration capacity of the soil [67]. In NCP, the streamflow has dramatically decreased because of the human activities. Low level of groundwater and large evaporation potential in this semi-arid region often lead to lentic water in drainages [68]. High primary production in the stagnant drainage channels will make them a carbon sink.

The median bicarbonate concentration of aquifer systems in the NCP is approximately 241 mg L−1 based on representative samples from three plain types [36]. Using the conservative assumption that when groundwater reaches the surface, half of the bicarbonate (121 mg L−1) is converted to CO2 [11] which is equivalent to 87 mg L−1 CO2. Annual groundwater net depletion in NCP was estimated as a volume of 8.3 ± 1.1 km3 year−1 [7]. The estimated annual CO2 emissions due to groundwater depletion in the NCP are thus approximately 0.52 ± 0.07 Mt (1 Mt = 1012 g). The irrigation efficiency in this region is about 0.5 as reported earlier [69]. Thus, half of the pumped groundwater reaches the drainage channels. However, we do not know how much of the CO2 is set free already on the fields and which part of it is entering the drainage channels. Furthermore, no biological carbon fixation is considered. This makes our estimate a worst case scenario. Indeed, high CO2 emissions during irrigation the fields have been observed [70]. Agriculture is always an integrative practice with soil management, fertilization, and irrigation. Further studies are needed to understand how much groundwater CO2 is emitted during the agricultural management and how much ends up in the drainage channels.

Groundwater depletion in the NCP contributed a great number of CO2 which is quite close to the groundwater CO2 efflux in US and India (Table 3), and accounted for ~ 5% of the global groundwater extraction. CO2 released from groundwater depletion has not yet been included in the Chinese carbon inventory [71]. A previous study indicates that groundwater loss in the NCP was approximately 50 km3 within 8 years, which is greater than the capacity of China's Three Gorges Dam (39.3 km3), the world’s largest power station [7]. Thus, CO2 emissions from groundwater pumping should be considered in global carbon budgets [72, 73]. This is especially relevant, because groundwater depletion is expected to acceleration in the future [74].

Table 3 Comparison of CO2 emissions from groundwater pumping in the NCP with other regions of the world as well as with other CO2 sources in the NCP
Full size table

Conclusions

Up to now, carbon dioxide emissions from irrigation drainage networks remain poorly understood. Mesocosm experiments help to unravel the effect of different drivers by experimental manipulation and replications. Here, we applied controlled mesocosm experiments to investigate the impact of flow velocity and nutrients on carbon dynamics in irrigation groundwater entering drainage channels in an intensive agricultural landscape. High nutrients contributed to TOC generation and in turn reduce CO2 evasion on a short term. High flow velocity, on the other hand, promoted rapid CO2 evasion at high nutrient level. Overall, TOC production counteracted CO2 evasion, whereas both of the flow velocity and nutrients stimulate DIC depletion. For a proper quantification of the real CO2 emissions from the entire network, information about flow velocity and residence time of the water in the ditch network as well as information about the fate of organic matter formed in the ditches is necessary. Controlled mesocosm experiments are a useful tool to disentangle the interaction of the various processes involved.