Spatio-temporal variations in lateral and atmospheric carbon fluxes from the Danube Delta

River deltas, with their mosaic of ponds, channels and seasonally inundated areas, act as the last continental hot spots of carbon turnover along the land–ocean aquatic continuum. There is increasing evidence for the important role of riparian wetlands in the transformation and emission of terrestrial carbon to the atmosphere. The considerable spatial heterogeneity of river deltas, however, forms a major obstacle for quantifying carbon emissions and their seasonality. The water chemistry in the river reaches is defined by the upstream catchment, whereas delta lakes and channels are dominated by local processes such as aquatic primary production, respiration or lateral exchange with the wetlands. In order to quantify carbon turnover and emissions in the complex mosaic of the Danube Delta, we conducted monthly field campaigns over 2 years at 19 sites spanning river reaches, channels and lakes. Here we report on the greenhouse gas fluxes (CO₂ and CH₄) from the freshwater systems of the Danube Delta and present the first seasonally resolved estimates of its freshwater carbon emissions to the atmosphere. Furthermore, we quantify the lateral carbon transport of the Danube River to the Black Sea.

We estimate the delta's CO₂ and CH₄ emissions to be 65 GgC yr⁻¹ (30–120 GgC yr⁻¹, a range calculated using 25 to 75 percentiles of observed fluxes), of which about 8 % are released as CH₄. The median CO₂ fluxes from river branches, channels and lakes are 25, 93 and 5.8 mmol m⁻² d⁻¹, respectively. Median total CH₄ fluxes amount to 0.42, 2.0 and 1.5 mmol m⁻² d⁻¹. While lakes do have the potential to act as CO₂ sinks in summer, they are generally the largest emitters of CH₄. Small channels showed the largest range in emissions, including a CO₂ and CH₄ hot spot sustained by adjacent wetlands. Thereby, the channels contribute disproportionately to the delta's emissions, considering their limited surface area. In terms of lateral export, we estimate the net total export (the sum of dissolved inorganic carbon, DIC, dissolved organic carbon, DOC, and particulate organic carbon, POC) from the Danube Delta to the Black Sea to be about 160 ± 280 GgC yr⁻¹, which only marginally increases the carbon load from the upstream river catchment (8490 ± 240 GgC yr⁻¹) by about 2 %. While this contribution from the delta seems small, deltaic carbon yield (45.6 gC m⁻² yr⁻¹; net export load/surface area) is about 4 times higher than the riverine carbon yield from the catchment (10.6 gC m⁻² yr⁻¹).

1 Introduction

In an attempt to improve global climate models, the role of rivers and their deltas and estuaries in the carbon cycle has received increased attention for more than a decade (IPCC, 2007). Back then, the perception shifted from rivers as being mere lateral conduits of particulate and dissolved carbon species to a so-called active pipe concept, where rivers are considered as being efficient biogeochemical reactors with the potential to release significant amounts of carbon as CO₂ and CH₄ directly to the atmosphere (Cole et al., 2007; IPCC, 2013). A multitude of global upscaling studies (e.g. Tranvik et al., 2009; Regnier et al., 2013; Raymond et al., 2013) estimated the riverine and lacustrine fluxes of CO₂ and CH₄ to the atmosphere on a persistently fragmentary database, considering spatial and temporal coverage – especially of headwater streams and large lowland rivers (Hartmann et al., 2019; Drake et al., 2018).

Along the land–ocean aquatic continuum, about 0.9–0.95 PgC yr⁻¹ are estimated to be transferred laterally by rivers to the ocean (Regnier et al., 2013; Kirschbaum et al., 2019). Half of the carbon exported to the ocean is in the form of dissolved inorganic carbon (DIC), while the other half consists of particulate and dissolved organic carbon (POC and DOC, respectively) in about equal shares (Li et al., 2017; Kirschbaum et al., 2019). Recent estimates suggest that about 50 % to > 70 % of the carbon inputs from terrestrial ecosystems degas as CO₂ and CH₄ along the way to the ocean (Drake et al., 2018; Stumm and Morgan, 1981; Kirschbaum et al., 2019; Cole et al., 2007), making this the most important export flux of terrestrial carbon from inland waters. While rivers could emit 0.65–1.8 PgC yr⁻¹ (Lauerwald et al., 2015; Raymond et al., 2013), lakes and reservoirs could add another 0.3–0.58 PgC yr⁻¹ (Raymond et al., 2013; Holgerson and Raymond, 2016). Earlier works on inner estuaries, salt marshes and mangroves estimate their contribution to be another 0.39–0.52 PgC yr⁻¹ (Borges, 2005; Borges et al., 2005). So, river deltas and estuaries seem to contribute almost equally to CO₂ and CH₄ emissions as lakes and reservoirs, despite representing only about one-sixth of their global surface area (Cai et al., 2013; Holgerson and Raymond, 2016).

Deltas and estuaries represent hot spots of carbon turnover and CO₂ and CH₄ emissions due to the high nutrient load, large productivity and seasonal flooding. However, differences in geomorphology, anthropogenic alterations, complex hydrology and the influence of tides are just a few of the factors which make it very difficult to compare different deltaic and estuarine systems amongst each other (Galloway, 1975; Postma, 1990). Dürr et al. (2011) attempted to classify this diverse group of coastal habitats, which led to lower global emission estimates of 0.27 ± 0.23 PgC yr⁻¹ for CO₂ and 0.0018 PgC yr⁻¹ for CH₄ (Laruelle et al., 2010; Borges and Abril, 2011). These studies, however, did not explicitly consider the deltas and inner estuaries of large rivers such as the Amazon, Changjiang, Congo, Zambezi, Nile, Mississippi, Ganges or Danube.

The close connection of river deltas to adjacent wetlands has the potential to fuel CO₂ and CH₄ emissions. Almeida et al. (2017) show that peak concentrations of CO₂ in the Madeira River, a tributary of the Amazon, are linked to extreme flood events, and riparian wetlands in the Amazon basin have been identified as significant sources for the outgassing of terrestrial carbon in the form of CO₂ (Richey et al., 2002; Mayorga et al., 2005; Abril et al., 2014). Global wetlands were estimated to contribute 1.1 PgC yr⁻¹ (Aufdenkampe et al., 2011) to the carbon emissions in the land–ocean aquatic continuum. The uncertainty of these estimates is large, due to the difficulty in delineating global wetland areas (Tootchi et al., 2019) and the complex interaction between potential emissions and carbon uptake by vegetation and soils (Hastie et al., 2019). While the lower river basins of the Amazon, Mississippi and Zambezi have been subject to CO₂ and CH₄ evasion studies (Sawakuchi et al., 2014; Dubois et al., 2010; Teodoru et al., 2015), others, such as the Nile and Danube, remained unchartered territory in that respect. Both the Nile and Danube rivers represent one end of the river delta spectrum since they show little exposure to tidal action. Therefore, these deltas experience seasonal flooding, instead of (semi-)diurnal flooding determined by tidal action. Flooding can, in addition to groundwater drainage and surface runoff, transport substantial amounts of terrestrial carbon to aquatic systems (Abril and Borges, 2019). We thus anticipate seasonal variability in CO₂ and CH₄ emissions and in lateral carbon transport from the Danube Delta to the ocean.

In this study, we estimate delta-scale atmospheric CO₂ and CH₄ emissions for the Danube Delta and the lateral carbon transport of the Danube River to the Black Sea. We hypothesized that the hydromorphology of the different waterscapes would influence the outgassing behaviour of greenhouse gases by governing gas exchange and biogeochemical processes. The resulting differences in atmospheric fluxes would require treating the waterscapes separately in the upscaling process. Furthermore, we anticipated that the seasonality of the flooding affects both atmospheric and lateral fluxes.

To capture this spatial and temporal variability, we conducted a systematic study covering 19 sites in the Danube Delta over 2 years, with monthly sampling intervals. Based on this time series, we address the systematic differences between the delta's main waterscapes (river branches, channels and lakes) to classify different open-water sources for greenhouse gas emissions and dominating biogeochemical processes. Furthermore, we estimate lateral and atmospheric carbon fluxes, considering the spatio-temporal variability, discuss uncertainties linked to the upscaling process and compare the estimates to other major river systems.

Figure 1 Sampling stations in the Danube Delta, Romania. Near Tulcea, the Danube River splits into three branches, namely Chilia, Sulina and St George. Station 16 was removed from the study because of limited access during lower water level (clogged access channel). Shape files for map creation in QGIS are adapted from https://mapcruzin.com/ (last access: 13 December 2016). Contains information from https://www.openstreetmap.org/, which is made available under the Open Database License (ODbL) at https://opendatacommons.org/licenses/odbl/1.0/.

2 Methods

2.1 The Danube Delta

The Danube Delta is the second-largest river delta in Europe after the Volga Delta. It is located on the Black Sea coast in eastern Romania and southern Ukraine (Fig. 1). Close to the city of Tulcea, the Danube River splits and forms the Chilia, Sulina and St George branch (or Sfantu Gheorghe in Romanian). In the vast wetland area between the main river sections, the seasonal floods maintain an aquatic mosaic of reed stands and more than 300 shallow flow-through lakes of different sizes, which are hydrologically connected to the Danube via natural and artificial channels (Oosterberg et al., 2000). Since 1998, the Danube Delta has been a UNESCO Biosphere Reserve, with nearly 10 % of the area being strictly protected and another 40 % of the total surface area being declared as buffer zones (UNESCO, 2019). While five of the larger lakes of the Danube Delta have been subject to CO₂ and CH₄ evasion studies in the past (Durisch-Kaiser et al., 2008; Pavel et al., 2009), the main branches of the river and the small channels are considered unchartered territory with respect to CO₂ and CH₄ concentrations and fluxes.

2.1.1 Hydrology

The hydrology of the Danube River, which drives water exchange with the delta, has a pronounced seasonality. Receiving meltwater from the Alps and Carpathians, the Danube shows peak discharge in spring from April to June (Fig. 2), whereas the discharge minimum occurs in autumn from September through November. December and January often show a small peak in discharge. The discharge provided by the Danube River drives the seasonal and annual hydrological changes in the delta. From 2000 to 2014, the Danube's average annual discharge was 6760 m³ s⁻¹ (ICPDR, 2018), which is a 3 % increase compared to the period from 1930 to 2000 (Oosterberg et al., 2000). In the delta region, the discharge splits into the different main branches as follows: Chilia – 53 %; Sulina – 27 %; St George – 20 % (ICPDR, 2018). Approximately 10 % of the Danube's total discharge (620 m³ s⁻¹; averaged over 1981–1990) flows through the delta, of which about 20 % (120 m³ s⁻¹) is lost via evapotranspiration (Oosterberg et al., 2000).

To assess the hydrological conditions during the time of observation with respect to the long-term average, we compared water level observations from Isaccea, Romania (INHGA; Feodorov, 2017), to the discharge data set from Reni, Ukraine (ICPDR, 2018). Reni is located about 30 km upstream of Isaccea, without any major tributary joining in between. Water level data from Isaccea were converted to discharge using rating curves created from paired water level and discharge data from the National Institute of Hydrology and Water Management (INHGA). The comparison shows that 2016 was quite an average year in terms of discharge (Fig. 2), while, contrastingly, the Danube had very low discharge in 2017, especially during the period between March and October. Average discharge in 2017 was 5237 m³ s⁻¹ or 23 % below the average flow calculated from the International Commission for the Protection of the Danube River (ICPDR) data set; hence, we refer to it as a dry year. Water temperature and conductivity of our sampling period were also, in general, comparable with data from the ICPDR's long-term monitoring (see the Supplement). Although water temperature measured during summer months in both 2016 and 2017 was up to 3 ^∘C warmer than the long-term mean, these values did not exceed the maximum temperatures measured in the last 20 years.

2.1.2 Categorization into river branches, channels and lakes

We categorized our sampling stations into three groups based on geomorphological characteristics, namely main river branches, lakes and channels. River branch stations are all located along the three main branches of the Danube River, exhibiting velocities of about 0.75 m s⁻¹ (Danube Commission, 2018), large hydraulic cross sections and frequent embankments. The category of lake refers to shallow (2–3.5 m) open-water bodies within reed bed areas, and five out of six sampling stations showed abundant macrophytes in summer. Natural and artificial channels represent the third category. They provide a surface water connection between the lakes and the river branches. We included old meanders of the Danube and small channels within the delta. Both of these features show a low flow velocity of up to 0.3 m s⁻¹, yet span quite a range in terms of surface area and depth. Accessibility by motor boat determined the sampling stations in lakes and channels and restricted our monitoring to deeper lakes and larger channels. Both lakes and channels are connected to adjacent reed beds and marsh areas. Very shallow or isolated lakes, which are not represented in our data set, may receive a significant part of their water from adjacent reed beds (Coops et al., 2008) and have a higher residence time of up to 300 d compared to the investigated lakes, which have an estimated residence time of 10–30 d (Oosterberg et al., 2000).

2.2 Sampling

Our research area was located in the southern part of the delta enclosed by the Sulina and St George branches, which we studied intensively in 2016 and 2017. We focused on the southern part of the delta, since it is less impacted by agriculture compared to the area north of the Sulina branch (Niculescu et al., 2017). Samples and in situ measurements were taken once per month at 19 stations (Fig. 1), representing river main branches (n=7), channels (n=6) and the larger delta lakes (n=6). The sampling stations in the channels and lakes cover both the fluvial (west of station 18; Fig. 1) and the fluvio-marine parts of the delta. In situ measurements and sampling with a Niskin bottle was carried out 50 cm below the water surface. Sample analyses were conducted at the Eawag laboratories in Switzerland.

2.3 Dissolved and particulate carbon species

For DIC measurements, filtered (0.2 µm) and bubble-free water samples were stored in 12 mL Labco Exetainers under cool and dark conditions until analysis with a Shimadzu TOC-L analyser. For the analysis of POC and DOC, water was filtered through 7 µm pre-combusted and pre-weighed Hahnemühle glass fibre (GF) 55 filters. The filters were stored at −20 ^∘C until analysis, when they were dried and weighed for total suspended matter, subsequently fumigated with HCl for 24 h to remove the inorganic fraction and analysed by EA-IRMS (elemental analyser) for organic carbon content, which we used to calculate POC. The filtered water was acidified using 100 µL 10M HCl and stored in the dark at 4 ^∘C until the analysis of DOC with a Shimadzu TOC-L analyser. Due to potential contamination during sampling, DOC data prior to May 2016 was discarded.

2.4 Dissolved gases

2.4.1 Concentration measurements

We used mostly field-based methods for the analysis of dissolved CH₄, CO₂ and O₂. In 2016, samples for CH₄ analysis were taken for laboratory-based analysis by gas chromatography. Bubble-free water was filled into 120 mL septa vials by allowing an overflow of approximately three times the sample volume before preserving the sample by adding CuCl₂. Depending on the expected concentrations, a headspace of 15–25 mL was created in the lab using pure N₂. Samples were equilibrated overnight at 23 ^∘C on a shaker, and the headspace was analysed using gas chromatography with a flame ionization detector (GC-FID; Agilent Technologies, USA). In 2017, we used 1 L Schott bottles to prepare headspace equilibration directly in the field, using air. Samples were transferred to gasbags and analysed in the field for CH₄ using an Ultraportable CH₄ $/$ N₂O analyser (Los Gatos Research – LGR). We corrected for atmospheric contamination during the processing by subtracting the amount of CH₄ introduced with the air during equilibration. As tests showed that there was no significant difference between the lab- and field-based methods (see the Supplement), we pooled the data in our analysis.

CO₂ concentrations were measured in the field using a syringe headspace equilibration of 30 mL sampling water with 30 mL air. The syringes were shaken for 2 min and allowed to equilibrate before the transfer of the headspace into a dry syringe and analysis in an infrared gas analyser (EGM-4; PP Systems). The method is explained in more detail in Teodoru et al. (2015).

Dissolved O₂ concentration was measured in situ using a YSI ODO probe. The sensor was calibrated daily using water-saturated air and cross-checked with oxygen readings from a YSI Pro Plus multimeter sensor. We measured local in-stream respiration rates to evaluate if community respiration could sustain our measured CO₂ fluxes. The respiration rate was measured as O₂ drawdown over a 24 h period. For the measurement, six biological oxygen demand (BOD) bottles were filled with water samples, and three were measured immediately afterwards at t=0. The other three bottles were stored in the dark at approximately in situ temperatures, and the O₂ concentration was measured after 24 h. The O₂ consumption rate was derived from the time and concentration difference, assuming a linear decrease over time. We used this respiration rate to estimate the local CO₂ production rate by assuming a 1:1 aerobic respiration relation of O₂:CO₂. Ward et al. (2018) argue that respiration rate measurements in BOD bottles underestimate the respiration rate because microbial processes are limited by both the bottle size and the lack of turbulence, and they suggest a correction factor of 2.7 to correct BOD-derived respiration rates for size effects only or a factor of 3.7 for size and low turbulence effects. Applying these correction factors did not change the main point of our comparison between fluxes and CO₂ production rates.

2.4.2 CO₂ and CH₄ flux measurements

CO₂ and CH₄ fluxes were measured using a floating chamber. The chamber had an internal area of 829.6 cm² and an internal volume of 10 080 cm³, leading to a volume/area ratio of 12.15 cm. An aluminium foil coating minimized heating during deployment. CO₂ was routinely measured in the field over a 30 min period by coupling an infrared gas analyser (EGM-4; PP Systems) to the chamber in a closed loop. In 2016, CH₄ was sampled from the chamber by syringe and transferred overhead into 60 mL septa vials that had been pre-filled with a saturated NaCl solution until the liquid was replaced by gaseous sample. These discrete samples for lab analysis were taken at time t=0, 10, 20 and 30 min and analysed by GC-FID. In 2017, this laborious procedure was replaced by attaching the LGR analyser directly to the floating chamber.

Flux chamber measurements were conducted, unless conditions were too windy or boat traffic was too frequent in the main channel. In total, we took 265 flux measurements for CO₂ and 122 for CH₄. Of the latter, 91 measurements seemed to be without any significant influence of ebullition (i.e. R² of linear regression > 0.96; for more detail, see the Supplement) and are henceforth referred to as diffusive CH₄ fluxes. In the high-resolution LGR analyser time series, the influence of gas bubbles could easily be identified. We calculated the diffusive flux by fitting a linear regression to periods where data showed no influence of ebullition. In this case, the flux is calculated from the slope and the height of the gas volume in the chamber. In the discrete time series, it was hard to distinguish between diffusive flux and ebullition. When the linear regression of the discretely measured samples had an R²<0.96, we considered the flux measurement to be influenced by bubbles. In this case, we calculated the total flux by dividing the total concentration increase by the observation time, as we did to calculate the total flux of the LGR analyser measurements. A total of three cases with R²>0.96 showed fluxes >20 mmol m⁻² d⁻¹ and were thus also classified as total flux. Discrete time series showing a non-monotonous course (n=12) were excluded from further processing. Missing monotony can have several explanations, including sampling captured a bubble or a sample mix up.

2.4.3 Calculation of k₆₀₀

We used our CO₂ flux measurements to calculate the gas transfer coefficient k₆₀₀ as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{k}_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{F_{{\mathrm{CO}}_{\mathrm{2}}}}{(p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{water}}-p_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{air}})\cdot K_{H,{\mathrm{CO}}_{\mathrm{2}}}}}\text{(2)}& {\displaystyle }{k}_{\mathrm{600}}={\displaystyle \frac{k_{{\mathrm{CO}}_{\mathrm{2}}}}{(S{c}_{{\mathrm{CO}}_{\mathrm{2}}}/\mathrm{600}{)}^{-\frac{\mathrm{1}}{\mathrm{2}}}}},\end{array}$$

where $F_{{\mathrm{CO}}_{\mathrm{2}}}$ is the flux of CO₂, $p_{{\mathrm{CO}}_{\mathrm{2}}}$ is the measured partial pressure of CO₂ in water and air, respectively, and $K_{H,{\mathrm{CO}}_{\mathrm{2}}}$ is the solubility coefficient for CO₂ according to Weiss (1974). $S{c}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the Schmidt number for CO₂, calculated based on temperature (Wanninkhof, 1992). We estimated missing flux measurements using the median k₆₀₀ of the respective water type and the measured CO₂ concentrations.

Analogously, diffusive CH₄ fluxes were estimated from the individually calculated k₆₀₀, using the solubility coefficient from Wiesenburg and Guinasso Jr. (1979), the mean global atmospheric CH₄ mole fraction of 1.84 parts per million (ppm; Nisbet et al., 2019) and the Schmidt number for CH₄ from Wanninkhof (1992). We attributed the difference between this estimate and the total measured flux to ebullition.

Table 1 Surface area of the Danube Delta features. Assuming a 19 m channel width means the estimation of the surface area of the channels is on the lower end. The surface areas of freshwater and wetland do not add up to the total area since parts of the delta are covered by forest and agricultural polders.

^∗ This is based on shape files adapted from https://mapcruzin.com/ (last access: 13 December 2016) and contains information from https://www.openstreetmap.org, which is made available under the Open Database License (ODbL) at https://opendatacommons.org/licenses/odbl/1.0/.

2.5 Upscaling atmospheric fluxes to delta scale

Spatial upscaling of heterogeneous and scarce data is very difficult and handled in various ways in the literature. Like other authors in a global context (Aufdenkampe et al., 2011; Raymond et al., 2013), we believe that median fluxes give a more reliable representation of the fluxes in systems with large gradients. Based on the different characteristics of the three waterscapes, we estimated the delta-scale atmospheric CO₂ and CH₄ fluxes by multiplying the median flux of each waterscape with its respective area (Table 1). We did this separately for each month and summed up the results, considering the respective number of days per month. For example, the median annual flux from the rivers, ${\stackrel{\mathrm{‾}}{F}}_{\text{R}}$, was calculated as follows:

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{‾}}{F}}_{\text{R}}=\sum _{m=\mathrm{1}}^{\mathrm{12}}{F}_{\text{R},m}\times A_{\text{R}}\times n_{\text{days},m}\times {\mathrm{10}}^{\mathrm{3}},\end{array}$$

where F_R,m is the median flux in mmol m⁻² d⁻¹ measured in the river stations in month m, A_R is the area of the river branches in square kilometres (see Table 1) and n_days,m represents the respective number of days per month m. The factor 10³ is used to convert to the units of mol yr⁻¹. To obtain the annual flux from the channels, ${\stackrel{\mathrm{‾}}{F}}_{\text{C}}$, and the lakes, ${\stackrel{\mathrm{‾}}{F}}_{\text{L}}$, we proceeded in the same way. We converted the resulting annual fluxes of the different waterscapes from mol yr⁻¹ to GgC yr⁻¹ and GgCO₂ eq yr⁻¹, with the latter assuming a global warming potential for CH₄ of 28 over 100 years, i.e. neglecting climate feedback (IPCC, 2013). The total annual water–air flux, ${\stackrel{\mathrm{‾}}{F}}_{\text{tot}}$, from the delta was the sum of the following three fluxes:

$$\begin{array}{}\text{(4)}& {\stackrel{\mathrm{‾}}{F}}_{\text{tot}}={\stackrel{\mathrm{‾}}{F}}_{\text{R}}+{\stackrel{\mathrm{‾}}{F}}_{\text{C}}+{\stackrel{\mathrm{‾}}{F}}_{\text{L}}.\end{array}$$

We also performed this calculation using 25 and 75 percentiles instead of the median to assess the upper and lower boundaries of our estimate.

For a reliable upscaling of fluxes, we determined the surface area of each waterscape as precisely as possible (Table 1). We estimated the area covered by the Danube's branches by refining publicly available shape files for Romania and Ukraine (https://mapcruzin.com/, last access: 13 December 2016), using the OpenLayers Plugin in QGIS, which allowed a comparison of the shape file with satellite images. We used the same procedure for the lakes and arrived at the surface area reported by Oosterberg et al. (2000). Assessment of the surface area of the delta channels was more difficult as many of the small channels are hard to identify on satellite images. Generally, estimating the width of the channels is challenging due to emergent macrophyte coverage, which, depending on the image quality, blends in with the adjacent reed. Instead of mapping the channels, we therefore used the overall channel length reported by Oosterberg et al. (2000) and assumed an average channel width of 19 m, which means the resulting surface area is on the lower end. Especially the old, cut-off meanders of the Danube River (Dunarea Veche), which we also consider as belonging to the channel category, do have a much larger width ranging on the order of 100–200 m.

2.6 Import by the Danube River and export to Black Sea

To compare the delta's CO₂ and CH₄ emissions to the lateral transfer of carbon from the catchment to the Black Sea and the influence of the delta region, we also calculated the loads of dissolved and particulate carbon species transported by the Danube at the delta apex, F_D, and close to the Black Sea, F_BS. As a first step, we calculated the daily average load of each month, F_m, for the different carbon species, as follows:

$$\begin{array}{}\text{(5)}& F_m=C_m\cdot \stackrel{\mathrm{‾}}{Q_m},\end{array}$$

where C_m is the concentration of DIC, DOC or POC measured in month m, and $\stackrel{\mathrm{‾}}{Q_m}$ is the respective averaged daily discharge of month m. Since CH₄ showed much smaller concentrations (∼ factor 100–1000 with respect to DOC and DIC), we did not include it into the calculation. In a second step, we weighed F_m by the number of days per month, n_days,m, and took the sum over all the months of the year. The load transported by the Danube River upstream of the delta, F_D, was calculated based on the concentrations measured at station 1 (Fig. 1), which is located in the Tulcea branch close to the apex of the delta and represents the water signature from the catchment, as follows:

$$\begin{array}{}\text{(6)}& F_{\text{D}}=\sum _{m=\mathrm{1}}^{\mathrm{12}}{F}_m\left(\mathrm{st}.1\right)\times n_{\text{days},m}.\end{array}$$

Data from the stations in the three main branches close to the Black Sea (stations 3, 4 and 5; Fig. 1) were used to estimate the amount of carbon exported to the Black Sea, F_BS, as follows:

$$\begin{array}{}\text{(7)}& F_{\text{BS}}=\sum _{m=\mathrm{1}}^{\mathrm{12}}\left[F_m\left(\mathrm{st}.3\right)+F_m\left(\mathrm{st}.4\right)+F_m\left(\mathrm{st}.5\right)\right]\times n_{\text{days},m}.\end{array}$$

Stations 4 and 5 are located slightly upstream of the settlements of Sulina and St George to avoid measuring the effect of these two settlements. Station 3 is located in a small side arm of the Chilia branch, marking the border between Romania and Ukraine, which, during comparison measurements, showed the same water composition as the main branch.

In our data processing, we decided to exclude one unusually high POC value in April at the Sulina branch (station 4) from our load calculation as we assume it is caused by a high-discharge, high-turbidity event that does not represent the monthly mean well. Instead, we interpolated between March and May. For DOC, we replaced missing data from January to April 2016 with the measurements at the same stations in 2017, assuming that they are also good estimates for the previous year. This way, we arrived at DOC estimates that cover the same period as DIC and POC.

We calculated the lateral transfer of carbon between the Danube Delta and its river, F_lateral, by subtracting the load exported to the Black Sea via the three main branches, F_BS, from the load imported to the delta from the catchment, F_D, as follows:

$$F_{\text{lateral}}=F_{\mathrm{D}}-F_{\text{BS}}.$$

The resulting lateral flux, in our case, is comparably small, and we used Gaussian error propagation to estimate its range. The basis for the error propagation was the measurement uncertainties in the concentrations (0.5 % DIC; 4 % DOC; 10 % POC) and discharge (3 %, assumed), which were used to calculate the loads.

2.7 Statistical analysis

We used MATLAB R2016a and R2017b for the statistical analysis of the data set. The data were evaluated for normal distribution, using histograms and quantile-quantile plots. In case of O_2,sat, CH₄ and POC, data distribution improved towards normality using log transformation; however, the results were not fully satisfying. Levene's test revealed, furthermore, the heteroscedastic nature of our data. Results for tests of significant difference between the three aquatic categories, from the non-parametric Kruskal–Wallis test (De Muth, 2014) followed by a multiple comparison test after Dunn–Sidak, were therefore taken very cautiously. Given the non-normality of the data, we report median instead of mean values and give ranges as 25 to 75 percentiles or minimum to maximum measured values, as indicated.

Box plots shown in this paper indicate the 25 and 75 percentiles and the median. Outliers are detected using the interquartile range (> 1.5 × IQR). The whiskers indicate the minimum and maximum values that are not detected as outliers by this procedure.

Figure 2 Daily average discharge close to the apex of the Danube Delta. The dotted line and the shaded area show mean and minimum to maximum daily discharge, respectively, for the period from January 1997 to October 2009 at Reni, Ukraine (ICPDR, 2018). Blue and red lines show daily average discharge at Isaccea, Romania, in 2016 and 2017 (INHGA; Feodorov). Horizontal bars indicate the timing of sampling campaigns. The x axis ticks indicate the 15th day of the respective month.

Figure 3 Measured DIC, DOC and POC concentrations in the different waterscapes (river, channel and lake). (a, c, e) The 2-year observation period. (b, d, f) Seasonality of the data. Dotted lines connect median values. The x axis ticks indicate the 15th day of the respective month. Box plots indicate 25 and 75 percentiles and median; whiskers indicate maximum and minimum, with data > 1.5 × IQR shown as outliers.

3 Results

3.1 Dissolved and particulate carbon species

DIC concentrations measured during our study ranged from 1.6 to 4.2 mM (Fig. 3a). Median DIC concentrations were around 3.0 mM over the whole observation period, with channels showing 10 % higher and lakes showing 3 % lower median concentrations than the main river. In 2016, concentrations were lowest in August and highest in December in all three groups (Fig. 3b). In 2017, median concentrations were 10 % (rivers) to 20 % (channels and lakes) lower than in 2016.

DOC levels in the delta were about 1.8 times the concentrations observed in the river (Fig. 3c, d). Channels and lakes had very similar concentrations, and both showed a general increasing trend from May to October 2016, but in the river, concentrations already peaked in July 2016 and were lowest in October. Median concentrations were quite comparable for 2017, with a tendency towards lower values. DOC in the main river in August 2017 was nearly 30 % lower than in the previous year. Most of the year, DOC concentrations were nearly a factor 10 smaller than measured DIC concentrations.

In 2016, we observed the lowest median POC concentration in the channels (Fig. 3e, f). Median concentrations in both rivers and lakes were nearly twice as high compared to channels but showed a distinctly different seasonality. POC was highest in the main river from March to June, while it peaked in lakes during August to October, suggesting different carbon sources.

Figure 4 Here, k₆₀₀ (a), daily average discharge close to the apex (b) and measured concentrations of dissolved gases in the different waterscapes, i.e. river, channel and lake (c–h), are shown. (a, c, e, g)  Pooled data from 2 years. (d, f, h) Seasonal dynamics, with dotted lines connecting median values. The x axis ticks indicate the 15th day of the respective month. (c, d) CH₄ in 2016, with four channel values (ranging from 22.2 to 58.0 µM) and one lake station (12.5 µM) exceeding 10 µM cut off. (e, f) Dotted black line represents the equilibrium concentration of CO₂ at 15 ^∘C (18.2 µM). Box plots indicate the 25 and 75 percentiles and median; whiskers indicate maximum and minimum, with data > 1.5 × IQR shown as outliers.

3.2 Dissolved gases

3.2.1 Concentrations

During the entire monitoring period, CH₄ in water samples of the delta was always oversaturated with respect to atmospheric equilibrium concentrations of 0.0046 to 0.0023 µM at T=0 to 30 ^∘C (Fig. 4c, d). Median concentrations in the river samples were thus ∼ 100 times oversaturated (0.33 µM). The channels exhibited a more than 3 times higher median concentration than the main river (1.1 µM), with the highest concentrations in July to September 2016 (up to 59 µM). In contrast, the median concentration in the lakes exceeded the value of the main river only slightly (0.43 µM) yet with a much larger range. In all three subsystems, concentrations increased from February 2016 to maximum values in July to October 2016. In 2017, concentrations were lower in the channels compared to 2016.

Similarly to CH₄, we found CO₂ concentrations to be constantly supersaturated with respect to the atmosphere in the main branches of the Danube, ranging from 26 to 140 µM (Fig. 4e, f). The median concentration of 59 µM was more than 3 times as high as the equilibrium concentration of CO₂ at 15 ^∘C (18.2 µM). Channels showed a much higher range (2.4 to 790 µM), with a significantly higher median of 140 µM. During the entire monitored period, we encountered undersaturated conditions in this class at only two stations (17 and 18) in August 2017. Lakes, however, were undersaturated on 11 occasions in 2016 and 32 occasions in 2017. Dissolved concentrations in this category ranged from 0 to 95 µM, with a median of 28 µM.

In 2016, CO₂ concentration showed a pronounced seasonality in all three subsystems. In the main river, median CO₂ nearly doubled from January 2016 to April 2016 and subsequently decreased to reach levels around 60 µM. In 2017, no clear seasonal pattern emerged. That year, median values mostly ranged around 60 µM, with the lowest median concentration recorded in June (44 µM) followed by the maximum in July (81 µM).

Channels showed the largest increase in CO₂ during the warm season. Median concentrations increased more than 4 times, from 66 µM in February 2016 to 290 µM in July 2016. In terms of inter-annual CO₂ variability, 2017 showed a later and less pronounced increase in concentration (72 µM in March to 187 µM in May), followed by an earlier decline than 2016. From August 2017 to November 2017, median monthly concentrations ranged around 50 µM and were lower than the concentrations in the main river during this period. In general, CO₂ concentrations in the channels in 2017 were 18 % to 75 % below the values observed in 2016. We found the highest concentrations in the eastern part of the delta (station 10; Fig. 1), where concentrations reached around 360 µM in winter and up to 785 µM in summer 2016.

Compared to rivers and channels, lakes generally had the lowest CO₂ concentrations and showed a distinctly different seasonal pattern. Most of the observed lakes (stations 7, 8, 13 and 14) were undersaturated in the period from May to November 2016. CO₂ undersaturation in these lakes (including station 20) occurred 3 times more often and over a longer period, from March to December, in the drier year of 2017. In 2016, lakes showed highest median CO₂ concentrations in April (74.4 µM) and the lowest concentrations in July and August (20.5 and 14.6 µM, respectively). With the concentration increase in early spring, the decrease in summer and the following increase in autumn, the seasonal signal in 2016 recalls a sinusoidal curve. The pattern in the drier year, 2017, however, showed less variation with lower concentrations which were ranging from 0 to 71 µM.

O₂ saturation, as one might expect, often showed a mirror image to the CO₂ time series in all three systems (Fig. 4g, h). The main river was generally slightly undersaturated, with a median O₂ saturation of 93 %. O₂ saturation in river water ranged between 75 % and 109 % during the whole observation period. Median saturation in the channels was 14 % lower (79.5 %) and – as for CO₂ – covered a much broader range than in the main river. The lowest values observed were as low as 5 % O₂ saturation (0.4 mg L⁻¹) in July 2016, while maximum saturation reached nearly 150 % in August 2017. In winter, O₂ saturation in the channels was comparable with the river stations. Station 10 showed an exceptional behaviour and never exceeded a saturation of 72 % or 9 mg L⁻¹. O₂ saturation in the channels strongly decreased in the spring and summer months, resulting in concentrations of less than 2 mg L⁻¹ at station 9 in July 2016 and at station 10 from July to September 2016 and in June, July and October 2017. Contrastingly, most lakes showed a strong oversaturation of up to 180 % from April to October, resulting in a median saturation that slightly exceeded 100 %.

Table 2 Median and range of measured CO₂ and CH₄ fluxes (mmol m⁻² d⁻¹) and calculated k₆₀₀ values (m d⁻¹). Additionally, n states the number of measurements. The range indicates the minimum and maximum observations.

^a The data in this table rely only on measured $F_{{\mathrm{CH}}_{\mathrm{4}},\mathrm{dif}}$. Missing diffusive CH₄ fluxes for the upscaling were calculated from k₆₀₀. ^b Measurement uncertainty led to negative k₆₀₀ values in nine cases (i.e., twice in the channels and seven times in the lakes). These values were deleted manually; thus, n(k₆₀₀) is < n($F_{{\mathrm{CO}}_{\mathrm{2}}})$ for channels and lakes.

3.2.2 Measured atmospheric CO₂ and CH₄ fluxes

Median CO₂ fluxes were largest in channels (93 mmol m⁻² d⁻¹; see Table 2) where we also observed the highest overall flux of 880 mmol m⁻² d⁻¹. Lakes were the only locations that showed significant negative fluxes, i.e. CO₂ uptake during summer, when O₂ was strongly oversaturated.

The highest median diffusive fluxes of CH₄ were observed in the channels with 1.1 mmol m⁻² d⁻¹. Diffusive efflux from the river was generally lowest, while the lakes showed the largest variability, with a minimum of 0.03 and a maximum of 6.7 mmol m⁻² d⁻¹. Considerable ebullition occurred only in the delta lakes and channels, which accounted for ∼ 70 % of the total CH₄ flux.

The gas transfer coefficient, k₆₀₀, was calculated from the measured CO₂ fluxes. Median k₆₀₀ was lowest in the river branches and in the channels at 0.69 and 0.74 m d⁻¹, respectively (see Table S1). As lakes were more exposed to wind, median k₆₀₀ was considerably higher (1.2 m d⁻¹), and we observed the maximum k₆₀₀ of 8.6 m d⁻¹ in this category.

Figure 5 Flux rate and production rate of CO₂, as calculated from O₂ community respiration incubations, for selected river, channel and lake stations. Fluxes marked with asterisks were calculated from median k₆₀₀ of our observations in the respective waterscape. Dark purple bars represent measured respiration rates; light purple bars indicate the effect of a correction for measurement limitations using BOD bottles (factor 2.7; see Ward et al., 2018).

3.2.3 CO₂ production rate vs. CO₂ flux

We find respiration rates ranging between 0.8–390 mmol m⁻² d⁻¹ for rivers, while in the channels and lakes they ranged from 2.3–560 mmol m⁻² d⁻¹ and 1.0–350 mmol m⁻² d⁻¹, respectively (Figs. 5 and S7–S9). Median respiration rate is highest in rivers (54 mmol m⁻² d⁻¹), followed by lakes (48 mmol m⁻² d⁻¹) and channels (45 mmol m⁻² d⁻¹). Many stations showed a pronounced seasonality, with the highest respiration rates occurring mostly between July and October. Respiration rates, i.e. CO₂ production rates, generally exceed CO₂ fluxes in river and lake stations throughout the year (Fig. 5), which implies that local instream CO₂ production sustained the observed fluxes. At the channel stations, we frequently observed fluxes exceeding the local production, even if we account for the potential underestimation of the CO₂ production, which implies the presence of other CO₂ sources. This was most striking at station 10, the CO₂ hot spot, where CO₂ outgassing exceeded local respiration on average by a factor of 40. At the other channel stations (also see Fig. S8), there seems to be a seasonally occurring pattern. CO₂ fluxes exceed local production in the first half of the year, while they fall below local production for the remainder of the year. While this pattern is very distinct in 2016, it is less pronounced in the drier year of 2017, which suggests that the additional CO₂ source is linked to hydrology.

Figure 6 Annual greenhouse gas fluxes to the atmosphere obtained by upscaling the monthly median flux to the total area of each waterscape and taking the sum over all months (see text for details). Black vertical lines indicate the uncertainty and were calculated using the 25 percentile and 75 percentile, respectively, instead of median values for the calculation. CO₂ flux in panels (a) 2016 and (b) 2017. (c) Diffusive and (d) total CH₄ flux in 2016. Due to large data gaps, this calculation was not done for CH₄ in 2017. All fluxes are in GgC yr⁻¹. For tabulated values, see Table S2.

3.3 Upscaling atmospheric fluxes to delta scale

The upscaling of the freshwater CO₂ and CH₄ fluxes to the freshwater surface of the delta according to Eq. (4) led to a net CO₂ flux of 60 GgC in 2016 and less than half (23 GgC) in the drier year of 2017 (Figs. 6 and 7a; case “c”) when the overall contribution of the three compartments was lower and lakes turned into a net sink. The diffusive CH₄ flux (Fig. 7c) was 1 order of magnitude smaller than the CO₂ flux (Fig. 7a), but it increased three-fold when ebullition was considered (Fig. 7d).

Figure 7 Comparison of greenhouse gas fluxes from the delta's freshwaters to the atmosphere obtained by the different upscaling approaches, i.e. pooled, and cases “a” (discrimination by year), “b” (discrimination by year and waterscape) and “c” (discrimination by year, waterscape and month). Black vertical lines indicate the uncertainty when performing calculations using the 25 and 75 percentiles instead of median values. (a) CO₂ flux. (b) Diffusive and (c) total CH₄ flux. All fluxes are in GgC yr⁻¹. Bold y axis labels indicate the calculation approach (case “c”), which is shown in more detail in Fig. 6 for the individual contributions from rivers, channels and lakes.

Especially the CO₂ fluxes seem to be subject to considerable inter-annual variability (Fig. 7a, b), which highlights the need to discriminate between different years during the upscaling process. It is likely that the different hydrological conditions triggered different amounts of lateral inflow from the reed-covered wetlands and contributed to the large variability in CO₂ fluxes. For CH₄, this effect appears to be much smaller.

Considering the contributions from the different waterscapes shows that the river branches were the main source of CO₂ to the atmosphere in both years (Fig. 6a, b). Despite their small surface area (7 %), channels contributed 32 %–37 % to the total CO₂ flux. Lakes, on the other hand, switched from a net CO₂ source of 19 GgC in 2016 to a small net CO₂ sink of −3.3 GgC in the drier year of 2017. In 2016, the lakes emitted the largest share of CH₄, with 66 %, considering only diffusive fluxes (Fig. 6c), and 86 %, considering total CH₄ fluxes (Fig. 6d). Considering the global warming potential of CH₄ (IPCC, 2013), CH₄ was responsible for 17 % of the total 260 GgCO₂ eq yr⁻¹ emitted in 2016.

3.4 Lateral carbon transport

The annual import of carbon to the apex of the delta amounts to 8490 ± 240 GgC yr⁻¹ (Fig. 8). This flux consists mostly of inorganic carbon (DIC; 91 %), while DOC and POC comprise only small fractions of 6 % and 3 %, respectively. Lateral fluxes are highest in spring when discharge is highest. About 10 % of the Danube's water is channelled into the delta before reaching the Black Sea (Oosterberg et al., 2000); thus, we assume that 10 % of the annual carbon load of the Danube reaches the delta (i.e. 849 GgC yr⁻¹).

The water export from the delta, however, is poorly constrained. The balance between precipitation minus evaporation is negative, poorly quantified and quite variable. We therefore rely on the flux balance of the three branches to estimate carbon export from the delta. The resulting export to the Black Sea via the Danube's main branches amounts to 8650 ± 150 GgC yr⁻¹ and is less than 2 % higher than the inflow load reaching the apex of delta. The slightly higher load mainly relates to increased DOC levels reaching the main branches from the delta, especially during the spring flood in March and April. The relatively small fraction of water that passes through the delta changes the relative fraction of DOC and POC only marginally to 7 % and 4 %, respectively, while the largest fraction in the water reaching the Black Sea remains DIC (89 %; Fig. 8). DIC import and export is fairly comparable throughout the year, while POC export to the Black Sea strongly exceeded the imports from the catchment in April. DOC exports are highest in the first half of the year (see Fig. S5).

4 Discussion

4.1 The main waterscapes of the Danube Delta

As we had hypothesized, carbon dynamics differed significantly across the three different waterscapes. The non-parametric Kruskal–Wallis test, followed by the Dunn–Sidak test, showed that the median of the three classes is significantly different for concentrations of CH₄, CO₂, O₂ and DIC (see the Supplement). In the case of DOC, only the rivers differ significantly from the other two groups, while in the case of POC, only channels are significantly different. Rivers and lakes, however, may differ significantly in the quality of their POC, as observed by the seasonality of the signal, which shows that high POC in the river actually occurs during high discharge in spring, while high POC in the lakes occurs during algal blooms in late summer. The non-parametric Kruskal–Wallis test does not require normal distribution of the data, but it requires equal variance of the data groups investigated for the difference in median (Hedderich and Sachs, 2016). Our observations in the seasonal plots (Figs. 3, 4) support the results of the test. In most cases, the box plots do not overlap, indicating that the three groups are significantly different. For example, DOC is significantly higher in the delta lakes and channels due to the strong primary productivity of these systems. O₂ is significantly lower in the channels than in the other two categories due to the lateral inflow of oxygen-depleted waters from the wetland (Zuijdgeest et al., 2016; Zurbrügg et al., 2012). The large difference between the waterscapes, with respect to CO₂ and CH₄ fluxes, supports our approach to treating the waterscapes independently when upscaling the flux measurements to the total water surface of the delta.

4.2 Dominating processes

4.2.1 River branches

The main river branches of the Danube are mostly influenced by the hydrology and chemistry of the catchment, as shown by the comparison between the concentrations at the delta apex with concentrations in the three main branches close to the Black Sea. There is comparably little variation between the stations with respect to DIC, DOC and POC. At all sites, O₂ is slightly undersaturated most of the time, but we do not see a strong influence of the delta close to the Black Sea.

4.2.2 Channels

Carbon dynamics in the channels are strongly affected by the water source. The channels are connecting the river branches to the delta lakes. The direction of this connection depends primarily on hydrologic gradients between the delta and the main branches, which means that flow direction can reverse in individual channels and, thus, alter their chemical signature due to a change in the main inflow. Seasonally, the channels transport dissolved carbon into the delta and provide nutrients to the reed stands during the high-water season. During times of receding water levels in the main branches, the channels act as the delta's drainage pipes. The comparison between CO₂ fluxes and local CO₂ production rates (Fig. 5) shows that the high CO₂ fluxes in the channels are often not sustained by in-stream respiration alone, in contrast to what we observed in the river and lakes. While this discrepancy mainly occurs during high discharge in spring, it is most evident at station 10, where it occurs throughout the year of 2016. Station 10 is located in Canalul Vatafu-Împutita, Romania, at the border of a core protection zone of the biosphere reserve. During this study, it stood out as a CO₂ hot spot, responsible for the highest CO₂ concentrations (Fig. 4f). Additional CO₂-rich water inflows from adjacent wetlands could explain the large CO₂ fluxes, which exceed CO₂ production. The water at station 10 was always exceptionally clean, low in oxygen content and had a low pH, supporting the hypothesis of a pronounced input from the reed beds. During times of unusually low water levels, such as in August and September 2017, the lateral influx from the reed seems to cease (Fig. 5). The, at first glance, contradictory timing of the increased lateral inflow during increasing water levels at the other channel stations could be explained by a pressure wave. Water flooding the vegetated area in the west will push out old water with a long residence time in the vegetated area at the other edges further east. In general, channel water in the Danube Delta is therefore a mixture of three main sources, namely Danube river water, lake water and water infiltrating from the wetland. The importance of the individual water source depends on the location of the channel sampling sites and on the water levels, which trigger flooding or draining conditions.

Figure 8 Overview of carbon flux estimates in GgC yr⁻¹. The total area between the main branches is 3510 km² (see Table 1). Black and grey numbers refer to fluxes estimated during this study based on data from 2016. The following italicized values refer to estimates based on the data in the literature from different study periods (studies for carbon burial and primary production do not explicitly consider seasonality): ^∗ carbon burial in lakes, based on average sedimentation rate measured in seven lakes in the Danube Delta, with an organic carbon content range of 3 %–30 % (Begy et al., 2018); ${}^{\ast \ast }$ net CO₂ uptake of Phragmites australis upscaled to the area covered by the Scripo-Phragmitetum plant community (Zhou et al., 2009); ${}^{\ast \ast \ast }$ upscaled primary productivity of the Scripo-Phragmitetum plant community (Sarbu, 2006). The green area in the plot symbolizes the reed area without indicating all the locations of its occurrence.

4.2.3 Lakes

In the lakes, residence times of 10–30 d allow primary production and local decomposition of organic matter to become important factors driving carbon cycling. We observed abundant macrophytes like Ceratophyllum demersum and Elodea canadensis growing in spring and early summer, which, depending on lake depth, even reached the lake water surface. A change in the abundance of submerged vegetation to vegetation with floating leaves might be linked to changes in the CO₂ and CH₄ fluxes (Grasset et al., 2016). Around July, algal blooms coincided with a significant reduction in macrophyte abundance. This pattern seems to be reoccurring due to the eutrophic state of the delta lakes (Tudorancea and Tudorancea, 2006; Coops et al., 2008, 1999). During our observations, both macrophytes and algal blooms caused a drawdown of CO₂ and supersaturation in O₂ (Fig. 4f, h). The algal blooms also partly explain the peak in measured POC from July to November, which extended to most of the delta's channels (Fig. 3f). The degradation of the macrophyte biomass coincided with locally elevated CH₄ concentrations from July to October (Fig. 4d).

In constructed wetlands, macrophytes were found to influence the composition of methanogenic communities by affecting dissolved O₂ and nitrogen in the rhizosphere, which had a direct impact on the amount of CH₄ released to the atmosphere (Zhang et al., 2018). Potamogeton crispus, for example, which is also found in the delta lakes and channels, seasonally sustained CH₄ fluxes that were up to 3 times higher than CH₄ fluxes from Ceratophyllum demersum (Zhang et al., 2018). Studies showed that the plant community composition in the delta lakes shifted since the 1980s due to increasing eutrophication, which also led to an increase in Potamogeton species recorded in the delta (Sarbu, 2006). It remains unresolved whether this change in vegetation also affected the CH₄ release in the Danube Delta.

4.3 Uncertainties linked to the upscaling procedure

4.3.1 Spatial heterogeneity

In a hydrologically complex system like the Danube Delta, upscaling CO₂ and CH₄ is prone to several sources of uncertainties, most of them linked to the delta's small channels and lakes. First, the channel category showed a large range, not only in DOC and POC concentrations but also in dissolved gases and their fluxes. We attribute this primarily to the varying contribution from the three different water sources, with lateral influx from the reed stands drastically increasing the local CO₂ concentration and fluxes. One could thus argue that this group is too broad and should be refined. However, in a complex system like the Danube Delta, this is a laborious task since individual channels are known to reverse the flow direction (Irimus, 2006), and potentially, the amount of lateral inflow also depends on the hydrologic conditions in the main branches. The existing 1D hydrological model, SOBEK (DANUBS, 2005), could assist in delineating periods of reversed flow, but a detailed model for the exchange with the wetlands would have to be developed.

Second, the surface area of the channels is estimated based on the channel length, given in Oosterberg et al. (2000), and an assumed channel width of 19 m, which leads to an estimated surface area that we consider quite conservative. More exact mapping or better spatial data, which might exist with local authorities but was not at our disposal, could improve this estimate. A larger or smaller surface area attributed to the channel would influence the flux estimates from this category accordingly.

Third, we identified station 10 as a CO₂ hot spot, with concentrations reaching up to 22 000 ppm, during our study. The hot spot channel had an east to west orientation and was draining a core protection zone. Considering channels with these two criteria indicates that potential hot spots could account for up to 2 % of the channel length (see the Supplement) and contribute up to 20 % of the CO₂ and CH₄ fluxes of the channel category. The overall emissions from the channels (including hot spot channels) was decreased by 10 % to 30 % in this scenario, since considering the high fluxes separately lowered the median value used for the calculation of the channel fluxes. A first step to improve the upscaling would thus be to map the spatial distribution of dissolved gases in the delta. This would give insight on important questions linked to the hot spots – how many hot spots did we miss with our discrete sampling approach? What is their lateral extent? And how steep are the concentration gradients between hot spots and nearby sites?

Fourth, our study neglected small, hardly accessible and remote lakes. A study of lakes of various sizes in northern Quebec, Canada, revealed a strong, negative relation between lake CO₂ concentration (and fluxes to the atmosphere) and lake area, suggesting a higher CO₂ emission potential of smaller lakes compared to lakes with a large area (Marchand et al., 2009). Previous studies of lakes with a small area in the Danube Delta characterize them as very clear-water lakes (Coops et al., 1999) that have little or no surface water connection to the main branches (Coops et al., 2008), with increased water residence times and O₂ concentrations below 5 mg L⁻¹ during midday (Oosterberg et al., 2000). This indicated that these lakes, like the hot spot channel in this study, receive the majority of their water from the reed stands (Oosterberg et al., 2000). In contrast to the channels, which are wind sheltered by 2–4 m high reed stands, these small lakes provide a larger surface area and, thus, a larger wind fetch. Depending on the primary productivity in these lakes, better wind fetch, in combination with water contributions from the reeds, could result in higher fluxes to the atmosphere – at least compared to the larger lakes measured in this study. Based on the research in the literature, we estimate the area of potentially isolated lakes to be 99 km². Attributing these isolated lakes with channel-like flux properties would raise the total CO₂ and CH₄ emissions of the lakes several times and turn them from a potential CO₂ sink into a CO₂ source in 2017 (see the Supplement). The scenario, as such, represents an extreme case, but it highlights the potentially large contribution from small, thus far overlooked, lakes in the delta.

4.3.2 Seasonality

Seasonal data coverage is often not sufficient to address the seasonality of the fluxes, which might bias the estimates towards either higher or lower emissions. However, not only the under-representation of certain seasons or events but also the pooling of the data during the upscaling process influences the resulting estimates. In the following, we look at the effects of data pooling for our 2 year data set by comparing different upscaling approaches. In addition to the approach presented in Eq. (4), where we discriminate by year, month and waterscape (case “c”), we also calculated the yearly fluxes in more simple ways by either pooling all data (pooled), discriminating between years only (case “a”), and by discriminating according to year and waterscape (case “b”). In case “c”, where we considered individual months, data coverage of CH₄ did not allow the calculation for 2017. In all approaches, we treated the reed stands in the wetlands as a terrestrial part of the system, i.e. excluding them from the analysis.

For the Danube Delta, CO₂ flux estimates decreased when considering spatial heterogeneity and seasonality because the channel data, which showed the most pronounced seasonality and the highest fluxes, are treated independently and assigned to a comparably small area. Independent consideration of data from different years allows the exploration of the inter-annual variability, which is quite pronounced for CO₂ (Fig. 7a). CH₄ emissions tend to be higher in 2017, but the trend is not as clear, especially considering total fluxes (Fig. 7c; case “b”). The lower CO₂ flux in 2017 can be explained by the weaker connection of the wetland to the freshwater system of the Danube. We expect that, in 2017, most of the water exchange, especially during low-discharge conditions, between the river and the inner delta was along the channels as surface water connections, with comparably little water bypassing laterally through the wetland. While the CO₂ fluxes from the river were only marginally smaller than in 2016, channels emitted less than 50 %, and the lakes even turned into a net CO₂ sink in 2017 (Fig. 6a, b). The importance of the flooded vegetated area on CO₂ concentration in rivers was also found in the Congo and Amazon river basins (Borges et al., 2015, 2019; Amaral et al., 2019), where larger inundated areas correlated with higher pCO₂ values. In the case of the lakes, reduced lateral inputs from adjacent wetlands reveal their large CO₂ uptake potential. However, this might result in higher CH₄ emissions, as calculations according to case “b” indicate. Neglecting seasonality, diffusive CH₄ fluxes from the lakes were 0.3 GgC yr⁻¹ higher in 2017 (1.0 GgC yr⁻¹; data not shown).

Durisch-Kaiser et al. (2008) found Danube lakes to be sources of CO₂ and CH₄ to the atmosphere in both May and September 2006. Their measured fluxes fall well within the range of our observations. The comparison of data for corresponding months shows, however, that their CO₂ concentrations in May are, on average, twice as high as the ones we measured in 2016, while September concentrations are, on average, 18 % smaller. The higher fluxes in May could have been due to the aftermath of the severe flood, which reached Romania in the second half of April 2006 and inundated large parts of the delta, thereby promoting lateral exchanges.

4.4 Lateral and atmospheric carbon fluxes

The freshwaters of the Danube Delta export, in total, about 225 GgC yr⁻¹ (Fig. 8). About 40 % of this carbon is directly released to the atmosphere, while 60 % of the carbon is transported laterally to the Danube and subsequently to the Black Sea. However, the majority of the carbon reaching the Black Sea originates from the catchment (8490 ± 240 GgC yr⁻¹). The contribution from the delta is therefore comparably small, and the fraction of dissolved and particulate carbon species is only marginally changed by the delta. The anthropogenic alterations to the river's main branches, like the straightening and deepening to allow for commercial navigation, might be an explanation for this. The Sulina and the St George branch were especially strongly altered in that respect, which has increased the discharge along these branches and decreased the lateral exchange with the delta. Excavation of the channels furthermore increased the surface water connection between different features of the delta.

Considering the area between the three main branches (A_Delta=3510 km²; Table 1) and the catchment area (A_catchment=817 000 km²; Table 1), the deltaic carbon yield amounts to 46 gC m⁻² yr⁻¹, while the riverine carbon yield to the Black Sea is 11 gC m⁻² yr⁻¹. So, although the Danube Delta contributes only about 2 % to the total carbon load reaching the Black Sea, its role as a carbon source should not be underrated, as the carbon yield (net export/surface area) of the delta is about 4 times higher than the yield of the overall catchment.

In total, the Danube River and its delta supplied the Black Sea with 8650 ± 150 GgC yr⁻¹ in 2016, which fuels carbon emissions in the river plume. Based on concentration measurements in July 1995, Amouroux et al. (2002) estimated the CH₄ flux from the Danube River plume close to the St George branch to be 0.47 mmol m⁻² d⁻¹, which compares very well with the CH₄ flux we measured in the Danube River branches. As CH₄ concentrations in the river plume were 5 to 10 times higher than in the rest of the water column, the authors expect this flux to be fuelled by the carbon reaching the Black Sea from the delta. They estimate the total CH₄ emissions from river plumes in the Black Sea to be 28–52 GgC yr⁻¹, based on the total surface area of the plumes. Since the Danube River is providing more than 50 % of the total discharge and is thus the largest freshwater contributor to the Black Sea (BSC, 2008), the majority of this emission might be released from the Danube River plume. Assuming a share of 50 % of the total river plume emissions would mean that 8 %–16 % of the carbon laterally transported to the Black Sea might reach the atmosphere in the form of CH₄. This corresponds approximately to the share of DOC and POC transported to the Black Sea.

Table 3 Selected major rivers and their carbon fluxes to the ocean and atmosphere.

^a Coynel et al. (2005); ^b Teodoru et al. (2015); ^c Meybeck and Ragu (1997); ^d Li et al. (2017); ^e Sawakuchi et al. (2017); ^f flux from Sawakuchi et al. (2014) and area for upscaling from Sawakuchi et al. (2017); ^g Laruelle et al. (2010); ^h Borges and Abril (2011); ⁱ Jiang et al. (2019) – total DIC flux estimated using the same discharge as ^c; ^j Soltan and Awadallah (1995) – total flux estimated using the same discharge as ^c; ^k Badr (2016) total flux estimated; ^l Bianchi et al. (2007); ^m Bianchi et al. (2004); ⁿ Kirschbaum et al. (2019); ^o Moquet et al. (2016) – estimated from HCO${}_{\mathrm{3}}^{-}$ flux; ^p Druffel et al. (2005) – total flux estimated using the same discharge as ^a; ^q taken from this study.

The comparison of our lateral DOC and POC fluxes (see Table 3) to available estimates of lateral carbon transport of European rivers to the ocean (Ludwig et al., 1996; Dai et al., 2012) indicates that about 3 % and 4 % of the POC and DOC could be exported by the Danube River alone. On a global scale, the lateral export of POC compares to the amount exported by the Zambezi River (Teodoru et al., 2015) but is about 20 % lower than the export from the Nile, despite the much higher discharge (Meybeck and Ragu, 1997). Absolute DOC export, on the other hand, is about twice as high in the Danube compared to Zambezi and Nile (Teodoru et al., 2015; Badr, 2016). Differences in DOC and POC export are strongly correlated to catchment area or river discharge, while factors such as climate, forest cover, population density or seasonality also affect the respective export fluxes (Alvarez-Cobelas et al., 2012; Hope et al., 1994). Looking at the organic carbon export yields (see Table 4), we observe that this general trend also prevails for the selected rivers, yet the DOC yield of the Danube's catchment surpasses the one of the Mississippi. This might be due to the lower population pressure and lesser agricultural usage of the Danube Delta, potentially resulting in a better connection of the floodable land to the river. DIC yield, however, is strongly influenced by the lithology of the catchment via silica and carbonate weathering (Gaillardet et al., 1999). The DIC yields of the Mississippi and the Danube catchment, where siliciclastic and carbonate rocks are abundant, are also highest, especially in comparison to the Amazon, where a Precambrian basement covers a large part of the heavily weathered catchment. This might explain why the Danube is transporting as much as one-third of the Amazon's DIC load, while only having 3 % of its discharge (Moquet et al., 2016; Druffel et al., 2005).

Table 4 Annual discharge, catchment area and carbon yields of selected major rivers.

^a Yield calculated based on catchment area and lateral carbon flux to the ocean (see Table 3); ^b Dai et al. (2009); ^c Meybeck and Ragu (1997); ^d ICPDR (2018); ^e Tudorancea and Tudorancea (2006); ^f the “average literature value” as cited by Teodoru et al. (2015); ^g Badr (2016).

CO₂ concentration in large rivers positively correlates with DOC concentration (Borges and Abril, 2011), which can be explained both by simultaneous lateral inputs and by terrestrial organic matter degradation in these net heterotrophic systems. For the selected rivers, the positive correlation also roughly holds for the CO₂ fluxes. The CO₂ fluxes per unit area from the Danube are much smaller than the ones from the Amazon, but they are closer to those observed in the Mississippi, the Zambezi and the average deduced for estuarine systems (Jiang et al., 2019; Borges and Abril, 2011). Based on this correlation, we would expect the CO₂ fluxes per unit area for the Nile to be somewhere between the ones from the Amazon and the Zambezi (see Table 3). Sites with high CO₂ concentrations are also likely to have a high CH₄ content. However, the relation is more complex and not always straightforward (Borges and Abril, 2011). The CH₄ fluxes per unit area in the Danube Delta were comparable with those of the Zambezi River but exceeded the fluxes of the Amazon's large inner estuary reported by Sawakuchi et al. (2014).

4.5 The role of the wetland

Based on a literature review, Cai (2011) suggested that estuarine CO₂ degassing is strongly supported by the microbial decomposition of organic matter produced in adjacent coastal wetlands. While CO₂ produced in marsh areas and transported to the estuaries was lost to the atmosphere, riverine DIC and DOC content were not greatly altered. Also, several other studies highlight the impact of the lateral input of wetlands or floodplain-derived water on river water O₂ content (Zurbrügg et al., 2012) and in-stream CO₂ levels (D'Amario and Xenopoulos, 2015). Abril and Borges (2019) recently suggested that the active pipe concept of carbon transport in the aquatic continuum indeed needs to be extended to consider floodable and non-floodable land as separate carbon sources. This is in agreement with the present study, highlighting how an exchange with the wetland can raise CO₂ fluxes well above locally sustained in-stream respiration. In the following, we therefore assess the potential role of the wetland in this complex hydrological system.

The Danube Delta is dominated by the plant association of Scripo-Phragmitetum, which covers nearly 89 % of the total marsh area (1600 km²). Its net primary productivity ranges between 1500–1800 g m⁻² yr⁻¹ (Sarbu, 2006). Assuming a carbon content of 0.42 gC gBiomass⁻¹, determined by Greenway and Woolley (1999) for Phragmites australis, primary production in the reed amounts to 1000–1210 GgC yr⁻¹ (Fig. 8), which is about 8 times less than the carbon load transported by the river. A large fraction of the net carbon assimilation by the Phragmites stands is decomposed and released back to the atmosphere. In a Danish wetland, more than 50 % of the carbon was respired and released back to the atmosphere, with 48 % being released as CO₂ and 4 % as CH₄ (Brix et al., 2001). In the Danube Delta, the 50 % accretion rate would correspond to about 500 gC m⁻² yr⁻¹. However, net primary production and carbon accretion change seasonally with environmental factors such as temperature and irradiation. Accordingly, net CO₂ assimilation in the Danish study was limited to the warm season, from April to September, whereas CO₂ and CH₄ emission occurred during the whole year but with maxima of 0.2 mol m⁻² d⁻¹ during July–August. Qualitatively, we observed the same seasonality in CO₂ oversaturation in the channels that drain water from the Phragmites stands (Fig. 4f). For a wetland dominated by Phragmites australis in China, at a latitude comparable to the Danube Delta, Zhou et al. (2009) estimated the annual net uptake of CO₂ to be 62 gC m⁻² yr⁻¹. Scaled to the area of the Danube Delta, this would result in 99 GgC yr⁻¹ remaining in the delta, which is in the same order of magnitude as the total annual input of organic C from the catchment (79 GgC yr⁻¹). Similar to the Danish study, Zhou et al. (2009) also did not account for the potential lateral transport of carbon to adjacent water bodies. Our results show that channels in the Danube Delta are receiving carbon from the wetland, with peaks in CO₂ and CH₄ concentrations that match the maxima in the gross ecosystem production in China. Comparing the estimated carbon fluxes from the channels with the yearly carbon accumulation estimates of the wetland suggests that up to 20 % of the latter could be released to the atmosphere via lateral transport, assuming the carbon fluxes from the channels were exclusively sustained by the wetland. With a lag phase of about 3 months, the Danube Delta reed beds release peak concentrations of DOC and POC during October to November when the biomass in the reed stands start degrading (Fig. 3d, f).

Assessing the amount of carbon input needed to sustain the observed carbon fluxes in the delta by a simple mass balance approach shows that inputs need to be even higher (Eq. 8). For the mass balance, we consider the net export to the Danube River (F_Danube=160 GgC yr⁻¹) and to the atmosphere (F_atm=65 GgC yr⁻¹) and assume that sedimentation is predominantly occurring in the lakes of the delta (F_sedi). Begy et al. (2018) found averaged sedimentation rates in the lakes of the delta in the range of 0.84 g cm⁻² yr⁻¹. Carbon content in the sediment cores ranged between 3 % and 30 %, translating into a carbon burial rate of 65–650 GgC yr⁻¹ across all delta lakes. For the purpose of this simple balance, we neglect anthropogenic effects, e.g. removal of fish biomass or burning of the harvested reed areas during winter, and potentially associated carbon inputs.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& F_{\mathrm{In}}\approx -F_{\mathrm{Danube}}-F_{\mathrm{atm}}-F_{\mathrm{sedi}}\\ & F_{\mathrm{In}}\approx \mathrm{290}\mathrm{to}\mathrm{875}\mathrm{GgC}{\mathrm{yr}}^{-\mathrm{1}}.\end{array}\end{array}$$

Assuming that freshwaters are a net balanced system and these three fluxes represent all major export fluxes suggests that inputs of 290–875 GgC yr⁻¹ are required to sustain the export to the Danube, the atmosphere and the sediment. Since long-term carbon burial is most likely an order of magnitude smaller (DeLaune et al., 2018) than the decadal sedimentation rate, we expect the required input to be rather at the lower end of the determined range. Nevertheless, it still surpasses the potential contribution from the wetland, as estimated above, by a factor of 3. This might either indicate an underestimation of the lateral export from the wetland or significant contributions from other sources, such as the forest areas or anthropogenic inputs to the system from fish farms or waste water. In addition, emergent macrophytes that border both lakes and channels in the delta could play an important role since they fix carbon directly from the atmosphere but are decomposed in the water column.

5 Conclusions

The waterscapes in the Danube Delta differ significantly with respect to their carbon cycling. While the river is mainly influenced by the carbon signal provided by the upstream catchment, carbon loads and especially greenhouse gas concentrations in the channels are strongly affected by lateral inflow from adjacent wetlands. Local primary production and respiration, on the other hand, dominate the carbon dynamics in the delta lakes. Considering the spatial extent of the three different waterscapes and the seasonality of their effluxes, we estimate that 65 GgC yr⁻¹ (range: 30–120 GgC yr⁻¹) were emitted from the delta to the atmosphere in 2016. Considering the small surface area they cover (7 %), channels, in general, contributed disproportionately to the total flux (30 %). Small lakes without a direct connection to the main river could represent similar hot spots for greenhouse gas evasion to the channels. Overall, nearly 8 % of the total flux to the atmosphere was released as CH₄ and was mostly supplied by the lakes. Covering a full annual cycle and discriminating between the three dominant waterscapes of the delta, we reduce the uncertainty linked to seasonal and spatial variability. However, spatial estimates could be further improved by investigating the extent of hot spots, gradients between discrete sampling stations, the effect of more isolated lakes and channels of the delta and the inter-annual variability, which especially CO₂ seems to show.

We estimate that the Danube Delta receives about 850 GgC yr⁻¹ from the upstream catchment. The export surpasses these inputs with the net carbon source from the delta to the Black Sea, amounting to about 160 ± 280 GgC yr⁻¹. However, compared to the overall carbon transfer from the Danube catchment (8490 ± 240 GgC yr⁻¹) to the Black Sea, the contribution from the delta is about 2 % and will not significantly alter the bulk carbon composition of the river water. In terms of carbon yield, the contribution from the delta is about 4 times higher (45.6 gC m⁻² yr⁻¹) than the riverine carbon yield (10.6 gC m⁻² yr⁻¹).

In order to sustain the observed carbon fluxes from Danube Delta freshwaters to the atmosphere and the Black Sea while assuming a net balanced system, a minimum of 290 GgC yr⁻¹ would be required to be provided by the wetland realm or other sources within the Danube Delta.