Groundwater data improve modelling of headwater stream CO₂ outgassing with a stable DIC isotope approach

A large portion of terrestrially derived carbon outgasses as carbon dioxide (CO₂) from streams and rivers to the atmosphere. Particularly, the amount of CO₂ outgassing from small headwater streams is highly uncertain. Conservative estimates suggest that they contribute 36 % (i.e. 0.93 petagrams (Pg) C yr⁻¹) of total CO₂ outgassing from all fluvial ecosystems on the globe. In this study, stream pCO₂, dissolved inorganic carbon (DIC), and δ¹³C_DIC data were used to determine CO₂ outgassing from an acidic headwater stream in the Uhlířská catchment (Czech Republic). This stream drains a catchment with silicate bedrock. The applied stable isotope model is based on the principle that the ¹³C ∕ ¹²C ratio of its sources and the intensity of CO₂ outgassing control the isotope ratio of DIC in stream water. It avoids the use of the gas transfer velocity parameter (k), which is highly variable and mostly difficult to constrain. Model results indicate that CO₂ outgassing contributed more than 80 % to the annual stream inorganic carbon loss in the Uhlířská catchment. This translated to a CO₂ outgassing rate from the stream of 34.9 kg C m⁻² yr⁻¹ when normalised to the stream surface area. Large temporal variations with maximum values shortly before spring snowmelt and in summer emphasise the need for investigations at higher temporal resolution. We improved the model uncertainty by incorporating groundwater data to better constrain the isotope compositions of initial DIC. Due to the large global abundance of acidic, humic-rich headwaters, we underline the importance of this integral approach for global applications.

1 Introduction

Rivers and streams are the main carbon pathways from the continents to the oceans and thus constitute an important link in the global carbon cycle. In the process of transport, large amounts of carbon – mostly in the form of CO₂ – outgas from the water surface to the atmosphere (Cole et al., 2007; Aufdenkampe et al., 2011; Regnier et al., 2013; Wehrli, 2013). Globally, this form of CO₂ contributions to the atmosphere was estimated between 0.6 and 2.6 petagrams (Pg) of carbon per year (Raymond et al., 2013; Lauerwald et al., 2015; Sawakuchi et al., 2017; Marx et al., 2017a). The lower value of this range by Lauerwald et al. (2015) excluded streams of Strahler stream numbers below 3. This is because of sparse coverage of actual direct measurements of the partial pressure of CO₂ (pCO₂) in headwater streams. Note that the upper value of the range still lacks a representative contribution from headwater streams (Marx et al., 2017a).

The contributions of headwater streams are considered as a major unknown factor in these estimates of global carbon budgets for inland waters (Cole et al., 2007; Raymond et al., 2013). The main reasons are uncertainties in groundwater input as well as poorly defined surface areas and gas transfer velocities (Marx et al., 2017a; Schelker et al., 2016). In addition, pCO₂ and subsequent CO₂ outgassing fluxes typically decline rapidly from stream source areas to river sections further downstream (van Geldern et al., 2015; Stets et al., 2017). Poor definition of these gradients adds another uncertainty to the global carbon budget. The enormous number of small headwater streams making a significant contribution on a basin and thus on the continental scale combined with the scarcity of data led to the term aqua incognita (Bishop et al., 2008).

Various direct and indirect approaches to determine CO₂ fluxes exist. Most often fluxes are calculated from pCO₂ and gas transfer velocities (Teodoru et al., 2009; Raymond et al., 2012; Lauerwald et al., 2015; van Geldern et al., 2015). However, because of large variabilities of gas transfer velocities on different spatial and temporal scales, this type of determination remains controversial, especially for small-scale applications (Marx et al., 2017a; Regnier et al., 2013; Schelker et al., 2016). For applications in small streams and during changing flow conditions, direct methods such as floating chamber approaches also exhibit major drawbacks such as altered outgassing behaviour because of artificially created currents inside anchored chambers (Lorke et al., 2015; Bastviken et al., 2015). In addition, rapid downstream losses of CO₂ often imply that CO₂-rich groundwater inputs are lost before actual measurements can take place (Reichert et al., 2009).

Recent approaches have used stable carbon isotopes of dissolved inorganic carbon to reliably quantify CO₂ outgassing from streams and rivers (Polsenaere and Abril, 2012; Venkiteswaran et al., 2014), in which the model by Venkiteswaran et al. (2014) is a parsimonious, simpler version of the model by Polsenaere and Abril (2012). Both apply inverse modelling to calculate the amount of CO₂ lost upstream of a sampling point within a stream or at a catchment outlet. One clear advantage when compared to conventional methods is that these stable isotope approaches account for the potentially high CO₂ outgassing upstream of any sampling point. Moreover, they incorporate groundwater seeps in first-order headwaters, particularly at low discharge (Polsenaere and Abril, 2012). These factors are typically not covered by conventional methods.

The integrative models exploit the fact that stable isotope ratios of dissolved inorganic carbon (expressed as δ¹³C_DIC) in stream water are controlled by ¹³C ∕ ¹²C ratios of its sources and the intensity of CO₂ outgassing. One important input is the stable isotope ratio of soil CO₂, which in turn depends on the plants' pathways used for photosynthesis and the organic matter sources fuelling plant and microbial respiration (Mook et al., 1983; Vogel, 1993). In general, this soil-internally produced CO₂ has a ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{CO}}_{\mathrm{2}}}$ value close to the initial substrate, which has a range from −30 to −24 ‰ for the most commonly occurring C₃ plants (Ehleringer and Cerling, 2002). After entering the stream, CO₂ outgassing to the atmosphere increases ¹³C ∕ ¹²C ratios in the remaining DIC pool because of the well-known equilibrium isotopic fractionation between CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (Myrttinen et al., 2015, 2012; Mook et al., 1974). This predictable and temperature-related process is calculated by the models of Polsenaere and Abril (2012) and Venkiteswaran et al. (2014). They are independent of the gas transfer velocity k and account for the upstream portion of a sampling site in headwater streams.

Figure 1 Location of the Uhlířská catchment and sampling sites modified from Sanda et al. (2014).

The aim of this work was to model stream CO₂ outgassing on the basis of the stable isotope approach by Polsenaere and Abril (2012) and to extend their method by including measured groundwater stable isotope composition in order to reduce modelling uncertainties. Our study utilises data from the well-studied Uhlířská catchment in the Jizera Mountains (Czech Republic) (Dusek et al., 2012; Sanda et al., 2014; Vitvar et al., 2016; Marx et al., 2017b). Since the background geology of the catchment consists of silicate rocks, carbonate weathering as a CO₂ source can be virtually excluded. The CO₂ saturation is then exclusively controlled by the mobilisation of terrestrial respired organic carbon and by the input of shallow groundwater (Humborg et al., 2010; Amiotte-Suchet et al., 1999).

2 Materials and methods

2.1 Study site

The Uhlířská catchment is situated in the northern Czech Republic, 9 km northeast of the city of Liberec (Fig. 1). The stream Černá Nisa flows in the catchment valley and is a tributary of the Lužická Nisa River that later merges with the Odra River and flows towards the Baltic Sea. The stream length in this experimental catchment is about 2100 m and water travel times are less than 1 h from the spring to the catchment outlet.

Table 1 lists the main characteristics of the Uhlířská catchment. The annual average precipitation exceeds 1200 mm yr⁻¹ and the annual average temperature is 5.5 ^∘C (1996–2009). Snow cover typically prevails during 6 months of the year, mostly between November and March (Hrncir et al., 2010). During our study period from September 2014 to April 2016 snow cover prevailed only between January and March. However, during snowmelt the monthly average discharge doubles (> 50 L s⁻¹) compared to the other months' discharges (Table 1).

The forest consists of a spruce monoculture (Picea abies) and isolated patches of larch, beech, and rowan trees (Sanda and Cislerova, 2009). Purely granitic bedrocks underlie this type of C₃ vegetation. The catchment has two basic types of soils. On the hillslopes, about 60–90 cm deep and highly heterogeneous soil profiles consist of dystric Cambisols, Podzols or Cryptopodzols that developed on weathered and fractured bedrocks. These soil types cover approximately 90 % of the catchment area. The valley bottom soils consist of a layer of peat of mostly Histosol types with depths up to 300 cm. The latter soils cover approximately 10 % of the catchment area, make up small wetlands along the stream and lie on top of fluvial material, which embodies the main perennial aquifer (Sanda et al., 2014).

Table 1 Characteristics of the Uhlířská catchment and Černá Nisa stream.

2.2 Water sampling and laboratory analyses

Between September 2014 and April 2016 shallow wetland groundwater and stream water were collected on a monthly basis and approximately weekly basis, respectively. Sampling points are shown in Fig. 1.

The stream discharge was determined at a V-notch weir at the catchment outlet (site UHL; see Fig. 1). Temperature (T), pH, and total alkalinity (TA) were determined in the field with portable HACH equipment, including a multi-parameter instrument and a digital titrator (all HACH Company, Loveland, CO, USA). The measured parameters exhibited a precision of 0.1 pH units, 0.1 ^∘C (2σ) for T and was better than 2 % (2σ) per 100 titration steps for TA (Marx et al., 2017b). Water samples were collected from the stream approximately 10 cm below the water surface and from boreholes with a peristaltic pump approximately 24 h after purging the boreholes.

Particulate organic carbon (POC) was sampled on 0.7 µm pore size glass fibre filter papers and analysed according to Barth et al. (2017). All other samples were filtered via syringe disk filters with 0.45 µm pore size (Minisart HighFlow PES, Sartorius AG, Germany) in the field. Before filtration, both syringe and membrane were pre-washed with sample water. Samples for the analysis of dissolved inorganic carbon (DIC) concentrations and isotopes were collected in 40 mL amber glass vials without headspace. The vials fulfil specifications of the US Environmental Protection Agency (EPA) and were closed with butyl rubber/PTFE septa and open-hole caps, with butyl rubber side showing towards the sample. Vials were poisoned with 20 µL of a saturated HgCl₂ solution to avoid biological activity after sampling. After collection, all samples were kept in the dark at 4 ^∘C until analyses.

Water samples were analysed for their carbon stable isotope ratio of dissolved inorganic carbon (δ¹³C_DIC) by an OI Analytical Aurora 1030W TIC-TOC analyser (OI Analytical, College Station, Texas) that was coupled in continuous flow mode to a Thermo Scientific Delta V plus isotope ratio mass spectrometer (IRMS) (Thermo Fisher Scientific, Bremen, Germany). The sample was reacted with 1 mL of 5 % phosphoric acid (H₃PO₄) at 70 ^∘C for 2 min to release the dissolved inorganic carbon (DIC) as CO₂. The evolved CO₂ was purged from the sample by helium. The Aurora 1030W system uses organic materials (see below) for the calibration to the VPDB scale and for concentration measurements. To oxidise these materials to CO₂, 2 mL of 10 % sodium persulfate (Na₂S₂O₈) was reacted for 5 min at 98 ^∘C. A trap and purge (T&P) system was installed for the analysis of low concentrations. Details of the coupling of the TIC/TOC analyser to IRMS are described in St-Jean (2003).

All values are reported in the standard δ notation in per mil (‰) vs. Vienna Pee Dee Belemnite (VPDB) according to

$$\begin{array}{}\text{(1)}& \mathit{\delta }={\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{reference}}}}-\mathrm{1},\end{array}$$

where R is the ratio of the numbers (n) of the heavy and light isotope of an element (e.g. n(¹³C) ∕ n(¹²C)) in the sample and the reference (Coplen, 2011). The data sets were corrected for instrumental drift during the run and linearity. Then they were normalised to the VPDB scale by two laboratory reference materials (C₄ sugar and KHP) measured in each run. The in-house reference materials were calibrated directly against USGS-40 (L-glutamic acid) and IAEA-CH-6 (sucrose) by using an elemental analyser (Costech ECS 4010). A value of −26.39 and −10.45 ‰ was assigned to USGS-40 and IAEA-CH-6, respectively. Concentration was determined from the signal of the OI Aurora 1030W internal non-dispersive infrared sensor (NDIR) and a set of calibration standard with known concentrations prepared from analytical (A.C.S.)-grade potassium hydrogen phthalate (KHP). Analytical precision based on the repeated analyses of a control standard (C₃ sugar) during all runs was better than ±0.3 ‰ for δ¹³C and better than 5 % for concentrations (1σ).

POC samples were analysed for δ¹³C_POC using a Costech Elemental Analyser (ECS 4010; Costech International, Pioltello, Italy; now NC Technologies, Bussero, Italy) coupled in continuous flow mode to a Thermo Scientific Delta V plus IRMS. The data sets were corrected for linearity and instrumental drift during the run. Values were normalised for carbon to the VPDB scale by analyses of internal reference materials (C₄ sugar and KHP) that were calibrated directly versus USGS-40 and USGS-41 (L-glutamic acid). Assigned values to USGS-40 and USGS-41 were −26.39 and +37.63 ‰ for δ¹³C, respectively. For precision and accuracy two laboratory standards (acetanilide and tartaric acid) were measured in each run. Precision, defined as the standard deviation of the control standard, was better than 0.1 ‰ (1σ) for δ¹³C_POC.

2.3 Calculations and model assumptions

At the Uhlířská headwater catchment the original streamCO₂-DEGAS model was applied to calculate stream CO₂ outgassing for different scenarios with varying values of river respiration (R) with 10, 19, 25, 50, and 75 % to test the model sensitivity to these values. In a second approach we modified the streamCO₂-DEGAS model as follows: instead of soil organic matter (δ¹³C_SOM) we used groundwater δ¹³C_DIC data to better constrain initial CO₂ values and to reduce the model uncertainty.

The Uhlířská catchment meets the assumption of the streamCO₂-DEGAS model with (i) stream waters being acidic with pH values between 3.9 and 6.7 (Table 2), and for our study we assumed that (ii) waters in the stream are unproductive. This means that secondary processes such as photosynthesis by algae or biofilms and organic matter (OM) degradation to CO₂ are neglected. This is a plausible assumption because high runoff and short residence times often leave insufficient time for substantial degradation of OM (Raymond et al., 2016; Catalan et al., 2016). However, the potential of temperature-dependent in-stream bio- and photodegradation (Demars et al., 2011; Moran and Zepp, 1997), particularly during summer, cannot be entirely excluded in the Černá Nisa stream.

The original streamCO₂-DEGAS model by Polsenaere and Abril (2012) demands the input variables of total alkalinity (TA), δ¹³C_DIC, pCO₂, and temperature (Table 2), as well as the proportion of river respiration (R) and the carbon isotope composition of δ¹³C_SOM. Note that TA and groundwater δ¹³C_DIC were only measured during monthly samplings and linear interpolated otherwise.

We used DIC concentrations to calculate ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (Dickson, 2007), and together with pH and T we were able to calculate pCO₂ values with the following equation (Plummer and Busenberg, 1982):

$$\begin{array}{}\text{(2)}& p{\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{{\mathrm{HCO}}_{\mathrm{3}}^{-}\times {\mathrm{H}}^{+}}{K_{\mathrm{H}}\times K_{\mathrm{1}}}},\end{array}$$

where ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ is the activity of bicarbonate, H⁺ is 10^−pH, K₁ is the temperature-dependent first dissociation constant for the dissociation of H₂CO₃ (all variables in mol L⁻¹), and K_H is the Henry's law constant in mol L⁻¹ atm⁻¹. The uncertainty of the calculation depends on the measurement uncertainties of DIC, pH, and T. The largest uncertainty is caused by pH measurements as they are on a logarithmic scale. The pH measurement uncertainty is typically smaller than ±0.1 pH units and causes a maximum uncertainty of ±21 % for pCO₂. We consider this a worst-case scenario, which is also indicated in Fig. 3.

R (0<R<1) is the proportion of CO₂ resulting from respiration in water along the entire stream course (waterlogged soils, stream waters, and sediments) and was approximated by a range between 0.1 and 0.19 (i.e. 10 to 19 %). This corresponds to the credible range of internal CO₂ production as a percentage of median stream CO₂ emissions from small streams (< 0.01 m³ s⁻¹) in the contiguous United States and was determined by Hotchkiss et al. (2015).

According to Polsenaere and Abril (2012) POC was used as representative of catchment soil organic matter. Polsenaere and Abril (2012) chose a δ¹³C_SOM of −28 ‰ that was determined as average annual δ¹³C_POC at their study site. In our study we also used the average annual stream δ¹³C_POC. Thus the δ¹³C_SOM was confined with −29.5 ‰ (Table 2). In addition, the average of wetland groundwater δ¹³C_DIC with an assumed isotopic fractionation of +1 ‰ for movement of CO₂ in waterlogged soils, groundwaters, river waters and sediments (O'Leary, 1984), served as input of initial δ¹³C_DIC (δ¹³C-DIC_init). Measured δ¹³C_DIC values ranged between −25.2 and −10.1 ‰ (Supplement).

The model results are the partial pressure of initial CO₂ before gas exchanges (pCO_2init) with the atmosphere start and the fraction of stream DIC that has degassed into the atmosphere ([DIC]_ex). [DIC]_ex corresponds to the CO₂ loss upstream of the sampling point and is given as a concentration (Polsenaere and Abril, 2012). CO₂ fluxes were calculated by multiplication with average discharges that were established from daily values. To allow for inter-catchment comparisons the carbon losses were normalised to the stream surface area. In addition, to avoid often imprecise stream lengths and surface areas, the carbon losses were also normalised to the catchment area.

For comparison, the gas transfer velocity adjusted to the in situ temperature (k_T, in m d⁻¹) can be calculated from model results when assuming that the water (pCO_2,aq) is the average between the modelled soil pCO₂ and the in situ pCO₂ at the catchment outlet:

$$\begin{array}{}\text{(3)}& k_{\mathrm{T}}={\displaystyle \frac{F}{K_{\mathrm{H}}\times \left(p{\mathrm{CO}}_{\mathrm{2},\mathrm{aq}}-p{\mathrm{CO}}_{\mathrm{2},\mathrm{air}}\right)\times M_{\mathrm{C}}}},\end{array}$$

where F is the modelled CO₂ outgassing (in g m⁻² d⁻¹), pCO_2,air the partial pressure of CO₂ in the atmosphere considered with ∼ 400 ppmV (ESRL/GMD, 2017), K_H is the Henry's law constant (in mol L⁻¹ atm⁻¹) and M_C the molar mass of C (12.011 g mol⁻¹). k_T was then converted into the normalised gas transfer velocity of CO₂ at 20 ^∘C (k₆₀₀, in m d⁻¹) according to

$$\begin{array}{}\text{(4)}& k_{\mathrm{600}}={\displaystyle \frac{k_{\mathrm{T}}}{{\left(\frac{{\mathit{\text{Sc}}}_{\mathrm{T}}}{\mathrm{600}}\right)}^{-\mathrm{0.75}}}},\end{array}$$

where Sc_T is the Schmidt number at the measured in situ temperature (Raymond et al., 2012).

All modelling approaches were executed via MATLAB (MathWorks, Natick, MA, USA).

2.4 Model input variables from the Uhlířská catchment

Stream (UHL) and wetland groundwater (DST, HST, P84) pH values were always below 6.7 during the entire study period (Table 2). In agreement with generally low stream pH values, the capacity of waters to buffer acidic inputs indicated by TA was generally low. The TA ranged around a median of 130 µmol L⁻¹ with a standard deviation (1σ) of ±101 µmol L⁻¹ at the catchment outlet (Table 2). Note that in the following the standard deviations express the variation over the sampling period and not the measurement error. The δ¹³C_DIC values had a median of −18.4 ‰ with a range from −25.2 to −10.1 ‰ (Table 2). The most negative values were measured in wetland groundwater with median values of −24.5, −28.7 and −24.3 ‰ in HST, DST, and P84, respectively (Fig. 1). The range of these groundwaters between −23.6 and −29.6 ‰ fits the range of δ¹³C in C₃ plants and associated soil organic matter. In the different water compartments (wetland groundwater and stream, Table 2), the DIC concentrations decreased with increasing δ¹³C_DIC values. The lowest DIC concentrations were found at the catchment outlet (31 to 483 µmol L⁻¹) and the highest concentrations in wetland groundwater (284 to 540 µmol L⁻¹). The median stream pCO₂ determined at the catchment outlet was 1374 ±710 ppmV with a range from 450 to 3749 ppmV. Values in wetland groundwater were typically higher with median values of 5590, 4440, and 6210 ppmV at HST, DST, and P84, respectively. Stream water temperatures showed a characteristic annual trend with a range from 0.8 to 13.8 ^∘C and a median of 4.5 ^∘C (Table 2).

Table 2 pH, total alkalinity (TA), temperature (T), pCO₂ values, DIC concentrations and ¹³C ∕ ¹²C ratios expressed as median values with ±1σ standard deviation and ranges given in parentheses. Site names are according to Fig. 1.

3 Results and discussion

Our study shows that the uncertainty of the respiration parameter (R) can be circumvented by incorporating wetland groundwater δ¹³C_DIC into the streamCO₂-DEGAS model. In a first step the original CO₂-DEGAS model was run with an R value of 10, 19, 25, 50, and 75 %. In a second step fractionation-corrected groundwater δ¹³C_DIC was incorporated into the model. Thus the isotope composition of the initial CO₂ was better constrained and the uncertainty on the isotope fractionation in soils was reduced.

In addition, we were able to estimate DIC export with modelled CO₂ outgassing and calculated lateral export of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, CO${}_{\mathrm{2}}^{*}$, and total DIC (Table 3). These data cover a period of 20 months and measurements took place at the catchment outlet, whereas modelling results relate to CO₂ outgassing between the stream source and the catchment outlet (UHL). The Supplement for this article contains detailed analytical data used for model calculations in this study. It is accessible through the journal's website and is additionally archived in the World Data Center PANGAEA (https://doi.org/10.1594/PANGAEA.889689) for long-term storage and free access.

3.1 Model sensitivity

The parameter of proportion river respiration (R) has a large impact on the results by the model. It is attributed to the average percentage of respiration occurring along the entire stream course (waterlogged riparian soils, stream waters, sediments, and the hyporheic zone). Thus, it is important to evaluate the model sensitivity to the assumed respiration. Modelled CO₂ outgassing and modelled initial pCO₂ values relate to an assumed R of 10 and 19 %, which corresponds to an average in-stream respiration as determined by Hotchkiss et al. (2015). However, this value does not account for respiration in groundwater and soil water. The contributions of these compartments are typically also considerable in headwaters. Although this type of respiration was not measured directly, we can assume a large potential of respiration in waters of the organic-rich wetland and of riparian soils with peak values during late summer and early autumn (Pacific et al., 2008). In addition, respiration in gravel bar waters along the stream (Boodoo et al., 2017) can lead to an exceedance of 19 % for R along the Černá Nisa stream.

Figure 2 Data of measured discharge, in situ pCO₂, δ¹³C_DIC, DIC concentrations, and CO₂ fluxes calculated via model equations in Raymond et al. (2012) for the Uhlířská catchment outlet (UHL) (a) and modelled soil pCO₂ and CO₂ loss, as well as k₆₀₀ calculated from model results and CO₂ outgassing normalised to the stream surface area (b).

In some cases the model was not able to produce reliable results and the convergence criteria were not fulfilled. The initial soil pCO₂ had to be extremely large (> 150 000 ppmV) for selected sampling events to reach the convergence of both δ¹³C_DIC and pCO₂ at the same iteration (Polsenaere and Abril, 2012). This was during few cases in February and March 2015 and during summer between May and September 2015. Moreover, for an R of 10 to 19 % during 22 to 23 cases, for an R of 25 and 50 % during 26 cases, and for an R of 75 % during 29 cases, the convergence criteria were not attained. Consequently, with increasing R a necessary convergence was more often not possible due to exceeding reasonable boundaries. These findings suggest that 10 to 19 % are reasonable R values for most months, except for June to August in 2015. Low respiration values are plausible for springtime; however, during summer an increased respiration can be assumed due to warmer temperatures and increased biological activity. One possible explanation for higher model uncertainty during summer is that the modelled CO₂ loss shows a strong non-linear dependence on in situ pCO₂ and δ¹³C_DIC (Polsenaere and Abril, 2012). Thus, the relative error of modelled CO₂ loss increases with the modelled CO₂ loss itself. This is particularly the case for low TA (e.g. 0.1 mmol L⁻¹; Polsenaere and Abril, 2012), where the δ¹³C_DIC increase is not buffered by the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ pool. According to Polsenaere and Abril (2012), large losses of CO₂ only occur with high in situ pCO₂ and with high δ¹³C_DIC. TA was low with 218 and 196 µmol L⁻¹ in the Uhlířská catchment during June and July 2015. In contrast, δ¹³C_DIC values and in situ pCO₂ were increased with on average −15.8 and −16.2 ‰ as well as 2225 and 2411 ppmV, respectively. Thus, increased δ¹³C_DIC values together with the model's non-linear dependence on in situ pCO₂ and δ¹³C_DIC may have caused an overestimation of modelled CO₂ losses and – as a consequence – led to failed convergence criteria during summer. Note that processes such as photosynthesis, which are not part of the model, may also lead to increased δ¹³C_DIC and lower in situ pCO₂ values. This may further increase the uncertainty of model results.

3.2 Gas transfer velocities

A huge benefit of the applied isotope approach is that the difficult to estimate gas transfer velocity parameter (k) is not needed for the calculation of CO₂ outgassing fluxes. On the other hand, k₆₀₀ (gas transfer velocity of CO₂ at 20 ^∘C, Eq. 4) values could be estimated from model results when assuming the water pCO₂ being the average between the modelled soil pCO₂ and the in situ pCO₂ at the catchment outlet (Table 3; Fig. 2).

For the original model calculated k₆₀₀ values varied from 0.5 to 7.0 and from 0.6 to 7.9 m d⁻¹ for R= 10 and 19 % (Table 3), with higher values mostly occurring during springtime and lower values during summer (Fig. 2). These values were higher than k₆₀₀ values that were determined in a lowland boreal landscape of Québec (Campeau et al., 2014) (Table 4). However, values were slightly higher than in the silicate Renet catchment in France (1.2–7.2 m d⁻¹) and were comparable to results from temperate silicate catchments in Germany (Halbedel and Koschorreck, 2013) (Table 4). They also equalled median values for boreal and Arctic streams (latitude 50–90^∘) determined by Aufdenkampe et al. (2011), whereas k₆₀₀ values for Sweden, for the United States, and for this study – all determined with model equations (i.e. Raymond et al., 2012) – showed higher values (Table 4). For modelling with groundwater data, we obtained lower k₆₀₀ values than with the original model (Table 3). Values ranged from 0.1 to 5.4 m d⁻¹ (Table 3) with higher values during springtime and lower values during summer (Fig. 2). The mean k₆₀₀ value was similar to the value of 2.2 m d⁻¹ calculated from gas transfer coefficients for propane in the temperate Zillierbach stream in Germany (Table 4).

Table 3 DIC partitioning according to the streamCO₂-DEGAS model (Polsenaere and Abril, 2012) and calculated from DIC measurements from September 2014 to April 2016 in the Uhlířská catchment expressed as median/average values ±1σ standard deviation and ranges given in parentheses. The total DIC export corresponds to the sum of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export, CO${}_{\mathrm{2}}^{*}$ export, and CO₂ outgassing. Initial soil pCO₂ was calculated from DIC_init., pH and temperature. k₆₀₀ was calculated from model results.

^a Modelling with proportion river respiration (R)=10 and 19 %. ^b Modelling with groundwater data.

Table 4 Gas transfer velocities (k₆₀₀) normalised to a stream temperature of 20 ^∘C of low order streams. Values calculated from model results are displayed in Table 3.

^a Mean values. ^b Median values. ^c Running waters in Aufdenkampe et al. (2011) have < 60–100 m width.

3.3 Soil pCO₂

The streamCO₂-DEGAS model assumes that the modelled initial pCO₂ represents soil pCO₂ (Polsenaere and Abril, 2012). For the Uhlířská catchment, this would mean that soil pCO₂ values ranged between 460 and 106 770 ppmV for original model and between 460 and 131 050 ppmV for modelling with groundwater input (Table 3). These were mostly within the range of values reported by Amundson and Davidson (1990), who determined pCO₂ from 400 to 130 000 ppmV for the upper metre of soils around the world. For three acidic catchments in France, modelled soil pCO₂ varied between 2120 and 77 860 ppmV (Polsenaere and Abril, 2012). In a temperate hardwood-forested catchment Jones and Mulholland (1998) modelled soil pCO₂ values between 907 in winter and 35 313 ppmV in summer. Higher soil pCO₂ values were also modelled during summer in the Uhlířská catchment. The main reasons for elevated soil pCO₂ in warmer seasons are higher temperatures with enhanced soil respiration in summer (Jones and Mulholland, 1998).

A clear positive correlation between soil pCO₂ and CO₂ outgassing for modelling with groundwater inputs (r² > 0.96, Supplement) stresses the importance of soil pCO₂ values and their dynamics for stream CO₂ loss.

3.4 CO₂ outgassing

Figure 3a displays modelled CO₂ outgassing for the proportion of river respiration (R) inputs that were selected with 10, 19, 25, 50, and 75 %. Figure 3b displays results for the modified modelling with groundwater input. Shaded areas show pCO₂ measurement uncertainties. Points that plot outside the measurement uncertainty did not fulfil the convergence criteria (Polsenaere and Abril, 2012). They were thus replaced by interpolated values. Note that, for modelling with groundwater, data the fluxes where the convergence criteria were not fulfilled could become reduced by ∼ 50 % (Fig. 3b).

For modelling with R=10 % and 19 % the modelled CO₂ outgassing varied from 0.01 to 72.52 mg C L⁻¹ and showed a mean of 6.35 mg C L⁻¹ (Table 3). When normalised to the catchment area, fluxes showed a mean of 5.10 mg C m⁻² d⁻¹ with values from 0.03 to 57.07 mg C m⁻² d⁻¹.

The highest fluxes > 30 mg C m⁻² d⁻¹ were modelled in May 2015, September 2015, November 2015, and March 2016 (Fig. 2). Translating these results to annual carbon outgassing fluxes yielded a mean of 23.9 to 34.5 g C m⁻² yr⁻¹ for 12 consecutive months and R= 10 and 19 %. Normalised to the stream surface area, an average of 28.6 to 41.3 kg C m⁻² yr⁻¹ was outgassed annually.

When including shallow wetland groundwater data in the model, the mean CO₂ outgassing was 1.34 mg C L⁻¹ and varied from 0.003 to 85.40 mg C L⁻¹, with minimum and maximum values during April 2016 and June 2015 (Table 3). Corresponding fluxes with a mean of 8.8 kg C d⁻¹ and a maximum of 88.5 kg C d⁻¹ were smaller when compared to mean CO₂ losses of 952, 258, and 10 671 kg C d⁻¹ for similar organic-rich catchments in France calculated with the same model approach (Polsenaere and Abril, 2012). However, stream discharge in the French catchments was larger by at least 1 order of magnitude, thus yielding higher outgassing rates. When normalised to the catchment area, modelled CO₂ losses had a mean of 4.9 mg C m⁻² d⁻¹ and varied from 0.02 to 49.7 mg C m⁻² d⁻¹ in the Uhlířská catchment. These values translate to a mean of 5.9 g C m⁻² d⁻¹ and a range between 0.02 and 59.5 g C m⁻² d⁻¹ when normalised to the stream surface area. Those groundwater-improved CO₂ outgassing estimates yielded the same outgassing trend with slightly decreased values (Fig. 3).

The most used method to calculated CO₂ fluxes is via k values, which are predicted via slope, flow velocity, stream depth, and discharge (Raymond et al., 2012). Corresponding Eqs. (1) to (6) yielded a mean flux of 8.1 g C m⁻² d⁻¹ with a range between 0.6 and 31.1 g C m⁻² d⁻¹ for the Uhlířská catchment. This mean value was larger than model results, whereas the range of fluxes was smaller (Table 3; Fig. 2). Huotari et al. (2013) and Wallin et al. (2011) observed large uncertainties of calculated k values and corresponding outgassing fluxes on the temporal scale. For instance, in headwater streams Wallin et al. (2011) determined errors of up to 100 % in median outgassing rates compared to measured k values. Annual carbon outgassing for modelling with groundwater data yielded 34.9 kg C m⁻² yr⁻¹ when normalised to the stream surface area and was in the range between 32.7 and 42.9 kg C m⁻² yr⁻¹ of annual values obtained by equations for k in Raymond et al. (2012).

Results indicate that fluxes according to Raymond et al. (2012) yielded reasonable estimates for annual fluxes, but they showed deficiencies in reproducing temporal variability (Fig. 2). The discrepancies are likely due to uncertainties in the selection of an appropriate k value in temporal highly variable headwater streams where variables such as slope, flow velocity, stream depth, and discharge are insufficient to predict the variability of k.

Figure 3 Modelled carbon dioxide loss via outgassing from the stream to the atmosphere upstream of the Uhlířská catchment outlet (UHL) based on the model by Polsenaere and Abril (2012) for proportion river respiration (R) between 10 and 75 % (a) and modified with measured groundwater δ¹³C_DIC (a, b). The areas shaded in blue (a) and red (b) indicate uncertainties in calculation of pCO₂ with ±21 %. Convergence criteria were not fulfilled for data points that lie outside the shaded area (see Sect. 3.1).

Figure 4 Modelled CO₂ loss via outgassing upstream of the measurement point (UHL) versus daily average discharge. Modelling results correspond to modelling with groundwater input.

Moreover, for CO₂ loss ([DIC]_ex) in mg C L⁻¹ and daily average discharge in L s⁻¹, a negative concentration–discharge relationship was observed (Fig. 4). Higher pCO₂, modelled soil pCO₂ and modelled CO₂ loss occurred during low flow, when relative contributions of CO₂-enriched groundwaters to stream waters were high. A similar concentration–discharge relationship was observed in a peatland stream, where deep soil and groundwater were considered as major CO₂ sources (Dinsmore and Billett, 2008).

3.5 Export of DIC proportions

Measured DIC together with modelled CO₂ outgassing from the stream surface to the atmosphere allowed the calculation of exported DIC species distributions. On an annual basis, the DIC proportions of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export, CO${}_{\mathrm{2}}^{*}$ (i.e. the sum of CO_2(aq) and H₂CO₃) export, and CO₂ outgassing had averages of 0.6, 1.2, and 7.1 to 10.3 mg C L⁻¹ for CO₂ outgassing with R=10 and 19 % (Table 3). This corresponds to the relative amounts of approximately $\mathrm{6}:\mathrm{12}:\mathrm{82}$ % ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export, CO${}_{\mathrm{2}}^{*}$ export, and CO₂ outgassing with respect to total inorganic carbon loss from the Uhlířská catchment outlet.

For modelling with groundwater, the relative proportion of CO₂ outgassing to annual DIC export was even higher. With concentrations of 0.6, 1.2, and 10.3 mg C L⁻¹ of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ export, CO${}_{\mathrm{2}}^{*}$ export, and CO₂ outgassing they were approximately $\mathrm{5}:\mathrm{10}:\mathrm{85}$ %. These values are comparable to relations found by other studies (Billett et al., 2004; Johnson et al., 2008; Davidson et al., 2010; Polsenaere and Abril, 2012) and largely differ from findings from a headwater stream in karstic bedrock, where < 30 % of DIC was outgassed as CO₂ (Lee et al., 2017). The determined DIC proportions suggest that a very large proportion (> 80 %) of DIC entering the headwater stream in the silicate Uhlířská catchment rapidly outgasses as CO₂.

4 Conclusions

Numerous studies have demonstrated the importance of CO₂ outgassing from rivers and streams to the atmosphere on global carbon budgets and pointed out the restricted data on headwater catchments (Raymond et al., 2013; Lauerwald et al., 2015; Aufdenkampe et al., 2011). Here we present a new study that successfully applied, validated, and modified the CO₂ degassing model by Polsenaere and Abril (2012) to a carbonate-free headwater catchment.

The modified isotope model was able to reproduce logical seasonal patterns of soil pCO₂ with a high variability of soil pCO₂ and CO₂ outgassing. It showed increased fluxes in summer and during snowmelt in 2015. Modelled CO₂ losses also negatively correlated with stream discharge. Results indicate maximum values of CO₂ outgassing from the stream to the atmosphere shortly before snowmelt. Modelled annual CO₂ outgassing was comparable to results obtained with model equations to calculate gas transfer velocities according to Raymond et al. (2012). However, results of the modified streamCO₂-DEGAS model showed larger variability, which indicates its potential to assess temporal CO₂ dynamics.

The model sensitivity to changing parameters of streamCO₂-DEGAS model, such as the proportion of river respiration (R), in situ pCO₂ and δ¹³C_DIC, was high. This indicates that the potential to assess temporal variations becomes compromised by a larger potential for errors. This was particularly the case during summer.

Because of its decreased uncertainty, future CO₂ modelling would further benefit from direct pCO₂ measurement instead of its calculation. Such direct in situ methods include submerged infrared gas analysis (Johnson et al., 2010), equilibrator systems (Polsenaere et al., 2013), an off-axis integrated cavity output spectrometer combined with a gas analyser (Gonzalez-Valencia et al., 2014), and non-dispersive infrared sensors inside floating chambers (Bastviken et al., 2015; Lorke et al., 2015).

We circumvented the uncertainty of the respiration parameter (R) by incorporating wetland groundwater δ¹³C_DIC into the model. By using groundwater data, the modelling results have shown substantial improvements when compared to modelling without groundwater data. We therefore stress the importance of adding more accurate groundwater measurements to such studies. This implies installation of more test wells in headwater catchments or the alternative use of local spring water as a proxy for groundwater. In addition, modelling at higher temporal resolution, particularly at the beginning of snowmelt, is needed to better reproduce dynamics and quantities of CO₂ outgassing.