A boost on the final stretch: intense river metabolism and wetland discharge increase aquatic CO₂ dynamics in the Danube Delta

Many river deltas are aquatic hot spots for carbon dioxide (CO₂) emissions to the atmosphere. Their patchwork of wetlands, lakes, channels, and river reaches often complicates the analysis of CO₂ sources such as ecosystem respiration or lateral water transfer. Sensing techniques offer the opportunity of measuring the CO₂, O₂ and DIC concentrations at high temporal resolution for periods from days to months. Such time-series allow quantification of diurnal and seasonal cycles of river metabolism and lateral exchange. This study addresses the following general hypotheses: (1) Ecosystem metabolism intensifies when river water enters the slower flow paths through channels and lakes of a delta. (2) Wetland discharge from a delta significantly alters the oxygen and carbon dynamics in a river. (3) In such a case, average aquatic CO₂ emissions increase on the final stretch before a river reaches the sea. We tested these hypotheses based on measurements of the oxygen and carbon dynamics at 15 min time resolution obtained from sensor packages. They were deployed for two years in the three main river reaches of the Danube Delta in Romania and for an additional year in three channels within the delta. By combining covariance analysis and monthly averaging of 24 h cycles we found a factor 100 difference in the amplitude of daily O₂ and CO₂ fluctuations across different stations and seasons. Channels with slow flow paths within the delta exhibited 4–8 times larger median amplitudes in daily metabolic cycles compared to the upstream river station. Correspondingly, metabolic intensity was on average 3–4 times more sensitive to changes in water temperature and cloud cover within the delta compared to the main river. Discharge of O₂-depleted and CO₂-rich wetland water into the downstream river sections was most pronounced during spring floods with apparent mixing rations of up to 13 %–25 % depending on the station. In a delta channel draining wetland waters, average CO₂ supersaturation was almost an order of magnitude higher than in the Danube inflow. The combined effects of intense metabolism within the delta and wetland discharge doubled the CO₂ emissions near the river mouth compared to an upstream Danube station. At the landscape level, however, carbon drawdown is likely five times larger than aquatic CO₂ emissions. Based on a high-resolution timeseries spanning three years, this study demonstrates how connected wetlands enhance aquatic metabolism and associated CO₂ dynamics in a large, lowland river.

1 Introduction

Over the last decades, the paradigm of river systems as active reactors (Cole et al., 2007) replaced a pipe model for carbon and nutrient transfer from land to ocean (Degens et al., 1984). Global estimates of CO₂ emissions from rivers and streams are 1.8±0.25 petagrams of carbon per year (Raymond et al., 2013). Including wetlands and inundated areas to river corridors (Abril and Borges, 2019) has further improved the analysis of carbon (Zuijdgeest et al., 2016) and oxygen dynamics in rivers (Zurbrügg et al., 2012). Reviews of river metabolism (Hotchkiss et al., 2015) covering the land-ocean aquatic continuum from the headwaters to the coastal zone have compiled an extensive database of carbon fluxes and transformations for river systems but available data remain more limited for large rivers and for the heterogeneous environments of the coastal and deltaic zones. Estuaries and their coastal vegetation were recently recognized as an effective global carbon sink (Rosentreter et al., 2023) but lateral carbon transfer from wetlands plays a key role in CO₂ emissions of large river systems (Hastie et al., 2019). Our study aims to clarify the role of a large river delta in modifying the intensity of the regional river metabolism and the associated CO₂ emissions.

Because river deltas often enclose large wetlands, the final stretch of a large river before entering the sea may exhibit a markedly different carbon dynamics compared to the upstream part. Three hypotheses guided this study on the delta of the Danube River in Romania: (1) Ecosystem metabolism intensifies when river water enters the slower flow paths through channels and lakes of a delta. (2) Wetland discharge from a delta significantly alters the oxygen and carbon dynamics in a river. (3) In such a case, average aquatic CO₂ emissions increase on the final stretch before a river reaches the sea.

Testing these hypotheses requires monitoring and analysis of ecosystem metabolism in large river sections and small wetland channels, each with specific characteristics and challenges (Bernhardt et al., 2017): wide river sections with extensive open water areas are typically not shaded by riparian vegetation which reduces the seasonal effects of leaf cover. Furthermore, hydrological changes occur more gradually and with lower amplitude compared to headwater streams where disturbance regimes may induce abrupt ecosystem change such as the abrasion of productive biofilms (Sabater et al., 2016). On the other hand, lateral groundwater exchange in large river corridors and dynamic flooding of adjacent wetlands may play an important role for CO₂ emissions (Borges et al., 2019; Marzolf et al., 2022). Wetlands are biogeochemical hotspots with high rates of daily gross primary production (Rabaey et al., 2024). Therefore, estimates of CO₂ concentrations and their potential drivers in lowland rivers (Reiman and Xu, 2019) and river deltas (Huertas et al., 2017) can be improved by high-frequency measurements that allow resolving diel and seasonal cycles of photosynthesis, respiration, and transfer across system boundaries (Battin et al., 2023).

Progress in sensor technology (Rode et al., 2016) now facilitates the autonomous registration of water quality parameters at high frequency over time periods of months. Seasonal and multi-annual observations are feasible with appropriate maintenance and calibration. Specifically, wetland connectivity (Maier et al., 2022), the role of seasonal flood pulses (Dalmagro et al., 2018) and differences between daytime and nocturnal emissions (Gomez-Gener et al., 2021) can be monitored by in situ sensors. Covariance analysis of paired O₂ and CO₂ measurements from sensor deployments provides insights into the biological, chemical and physical forcing of aquatic ecosystem processes over time (Vachon et al., 2020). Analysing the temporal dynamics of paired O₂-DIC measurements offers the advantage of removing buffer effects and simplifying carbon budget calculations (Shangguan et al., 2025). This combination of sensor measurements with statistical analysis has proven highly valuable in recent studies (Rocher-Ros et al., 2025) for identifying drivers of CO₂ dynamics in US rivers (DelVecchia et al., 2023) and in the lower Ganges (Haque et al., 2022), for assessing the balance between river metabolism and CO₂ emissions (Solano et al., 2023), for identifying the contribution of lateral inputs to the carbon budget of rivers (Marzolf et al., 2022), and for calculating the fraction of bicarbonate that supports river photosynthesis (Aho et al., 2021).

This study in the Danube Delta builds on previous analyses of discrete monthly sampling campaigns during two years (Maier et al., 2021) which revealed median CO₂ fluxes in the main branches of the Danube of 25 $\mathrm{mmol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and four times higher values in the canals connecting wetlands and lakes. A subsequent on-site mapping study of dissolved gases at high spatial resolution by membrane-inlet mass spectrometry (Maier et al., 2022) revealed hot-spots in emission rates caused by wetland discharge to the canal systems, by plant mediated gas transfer via O₂ ebullition in lakes, and by excess air formation in reed stands. Here we report results from a two-year deployment of commercially available multiprobe sensors (YSI EXO2) recording every 15 min in the three main stems of the Danube in the Romanian part of the delta near the Black Sea coast which was followed by a one-year deployment in selected delta channels. Time series of dissolved O₂, CO₂ and DIC obtained from these campaigns allowed us to address the following specific research questions: (1) How does the intensity of river metabolism differ between the open waters of the delta compared to the upstream Danube River? (2) How do lateral inputs from lake systems and reedbeds modify the O₂ and CO₂ dynamics of the main river reaches crossing the delta? (3) What are the overall effects of river metabolism and wetland discharge on aquatic CO₂ emission rates in the Danube Delta? The first question addresses the intensity and characteristics of coupled O₂–CO₂ dynamics as indicators for differences in metabolic activity between river reaches and deltaic waters. The second question analyses O₂, CO₂ and DIC indicators for lateral inputs from lakes, channels and reed stands. Building on the evidence for metabolic activity and wetland discharge, the third question evaluates the consequences for aquatic supersaturation of CO₂ and resulting emission rates. By addressing these questions our study proceeds from the analysis of local O₂–CO₂ and DIC patterns and processes to ecosystem-scale emission rates. We show how high-frequency observations can identify the controls of carbon processing and aquatic CO₂ emissions in large, heterogeneous river deltas. The results fill a gap in our understanding how wetland–river interactions may define carbon dynamics on the final stretch of the land–ocean aquatic continuum.

2 Materials and Methods

2.1 Study site and monitoring stations

The Danube River is the second largest river in Europe and its international catchment of 817 000 km² covers parts of 19 European countries. Originating in the Black Forest Mountains in Germany, the Danube flows southeast for over 2857 km before discharging into the Black Sea. The long-term mean discharge is 6360 m³ s⁻¹ and ranged from 2930–11 300 m³ s⁻¹ in the period 1921–2015 at the Ukrainian station of Reni (Romanova et al., 2019). Receiving meltwater from the Alps and the Carpathian Mountains, the hydrology of the Danube River shows a pronounced seasonality. In general, peak discharge lasts from April–June and low flow conditions prevail between September and November. Driven by rainfall in the lower catchment, a secondary discharge maximum usually occurs from December through January. During our study, seasonal air temperature changes varied from around 0 °C in winter to 20–25 °C in summer. The probability of cloud cover was higher in winter than in summer (Fig. 1a).

Figure 1 (a) Three years of environmental records for the Danube Delta 2016–2018. Top: mean daily discharge of the Danube at Reni, Ukraine, upstream of the Delta (available at https://www.danubehis.org/stations, last access: 27 October 2025). Gray area indicates the range of daily discharge for the period 1998–2017. Middle: average daily air temperature; bottom: average daily cloud cover at Tulcea Airport (available at https://rp5.lv/Weather_in_Tulcea_(airport), last access: 27 October 2025). (b) Monitoring locations for EXO2 probes in main branches and within the Danube Delta. Map: © OpenStreetMap contributors 2024. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Close to the mouth, the river splits into three arms that form the Danube Delta: Chilia in the north constitutes the border with Ukraine, Sulina in the middle was modified for maritime commercial navigation while Sfântu Gheorghe (St. George) limits the delta area in the south (Fig. 1b). The surface area of the Danube Delta between these river reaches is approximately 4150 km², of which more than 80 % lies within Romania. As the youngest and the longest among the three arms, the 120 km long Chilia branch carries more than 50 % of the Danube water, the 64 km long Sulina channel about 27 % and the 70 km long river reach of St. George contributes about 20 % to the total water discharge. Travel-time analysis based on our sensor data resulted in average flow velocities in the range of 0.7 and 1.0 m s⁻¹.

An intricate network of approximately 470 lakes, swamps, shallow water pools and about 3500 km of channels connects the delta's wetlands with the Danube River (Gomez-Baggethun et al., 2019). Although considered the most pristine delta in Europe, a series of hydrotechnical works in the Danube Delta, mostly for agriculture and fishery, led to doubling the length of the internal canals after 1980 (Gatescu, 1993). This, in turn, gradually increased the water flow through the delta from 260 m³ s⁻¹ between 1951 and 1960 to 620 m³ s⁻¹ in the period 1981–1990 (Bondar, 1994). Today, approximately 10 % of the total Danube discharge is estimated to enter the delta, of which about 20 % ($\sim \mathrm{120}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) may be lost via evapotranspiration in the wetlands. Reed (Phargmites australis) is the dominant species in the Delta covering ∼1600 km² or almost 40 % of the study area (Oosterberg et al., 2000). Large stretches of the channels are bordered by riparian forests composed mainly of black alder (Alnus glutinosa) and silkvine (Periploca graeca) (Oprea et al., 2024). Most of the lakes are eutrophic and classified into three types: (1) lakes at short distance to the river with high flushing rates and dense submerged vegetation, (2) large, relatively deep lakes like Puiu and Rosu with longer residence time and high Potamogeton-cover that collapses during winter, and (3) small, isolated lakes receiving “black-water” inflows from floating reedbeds and dominated by Ceratophyllum demersum and Nitellopsis obtusa (Coops et al., 2008). A high diversity of submerged and floating macrophytes is present in the canals and old meanders (Oosterberg et al., 2000).

In the first phase of this study, we deployed four sensor packages (EXO2 Multiparameter Water Quality Sondes, YSI) along the main branches of the Danube at the stations labelled Tulcea, Chilia, Sulina, and St. George (Fig. 1b). Tulcea served as the upstream reference station for water entering the delta. The Chilia branch receives additional input from a Ukrainian catchment with shallow lakes to the north, the Sulina channel serves as the main shipping route with lateral inputs from adjacent wetlands, and the St. George branch collects inflow mainly from wetlands on its left side. The monitoring stations recorded data every 15 min between November 2015 and February 2018. For consistency, all timestamps were recorded in Easter European Time (EET) without a switch to summertime (EEST). For further details on time intervals and locations see Fig. A1 and Table A1. One multiprobe was damaged by floating ice in December 2016 which limited the time series at Chilia station to one year. In the second phase from February 2018–December 2018, the remaining three EXO2 probes were installed in channels within the delta at Puiu-Rosu, Balanova and Busurca. Puiu-Rosu is influenced by outflow of a large and relatively deep lake complex, Balanova receives mixed input from lakes wetlands while Busurca is mainly connected to reedbeds with negligible lake influence. Note that Balanova refers to a local name along the Central Channel that connects the Sulina and St. George branches (Fig. 1b).

2.2 Measured and derived parameters

The EXO2 multiprobe sensors recorded water temperature, pH, specific conductivity, and dissolved oxygen every 15 min. Data of sensors for turbidity and fluorescent dissolved organic matter was not analyzed for this study. Before each measurement, the sensor caps were automatically cleaned by a wiper. The EXO probes were moored at approximately 1 m below the water surface. Data gaps occurred sporadically due to unexplained energy drain from batteries and when probes were removed to prevent damage by floating ice sheets during cold spells in winter. Once a month, the probes were removed from the water for cleaning and pH and O₂ sensors were calibrated using standard buffer solution and air-saturated water, respectively. Post-deployment, the data was drift-corrected using the calibration data. Crosschecks with discrete in situ measurements using the YSI Optical Dissolved Oxygen and Professional Plus Multiparameter instruments (Maier et al., 2021) showed good agreement. We used a low conductivity threshold ($<\mathrm{200}\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}$) to detect and delete invalid data points recorded at times when the sensor was out of the water for servicing.

We derived the alkalinity time series from specific conductivity data using the linear correlation between laboratory analysis of alkalinity and specific conductivity (SpCond) measurement as: alkalinity $\left[\mathrm{mmol}{\mathrm{L}}^{-\mathrm{1}}\right]=\mathrm{0.0057}\times$ SpCond [µS cm⁻¹] + 0.60 (R²=0.70, Fig. B1). The same approach has been applied to derive multidecadal trends for CO₂ evasion in the Loire River (Nguyen et al., 2025) and to estimate global riverine CO₂ emission rates (Raymond et al., 2013). To calculate CO₂ from pH and alkalinity, we used the CO2SYS MATLAB code version 1.1 with Matlab R2017a and 2022b (van Heuven et al., 2011). These calculations were based on dissociation constants for a temperature range of 2–35 °C (Cai and Wang, 1998). To obtain an upper and lower estimate for the derived CO₂ time-series, we propagated the uncertainty of the alkalinity-conductivity regression estimates and obtained an error of ±13 %. The calculated CO₂ values correlated well with discrete headspace measurements taken at the same time and location (R²=0.82, Fig. B1). Finally, we used solubilities (Weiss, 1974) and the temperature recordings from the EXO2 probes to calculate the deviation from atmospheric equilibrium in µmol L⁻¹ of the O₂–CO₂ pairs. We defined the frequent excess of CO₂ and the typical deficit of O₂ with regards to atmospheric equilibrium as exCO₂ and exO₂ [µM] with a positive or negative sign, respectively. For the correlation of daily O₂–CO₂ covariance data with potential physical drivers, we used averaged daily water temperature from the EXO2 sensors and cloud cover from a meteorological station in Tulcea (Fig. 1a). The time-shift of one hour during daylight saving time was removed to synchronize the meteorological data with the continuous recording of EXO probes.

Some data-filtering was necessary for the time series of Chilia (2016) which was sporadically influenced by water from an adjacent wetland and for Sulina (2017) which recorded disturbances from the Black Sea, respectively. We flagged these spikes in conductivity by comparing the time series with upstream Tulcea station and removed data if they deviated more than 40 µS cm⁻¹ from the baseline. This approach missed the beginning and the end of the spike events. To refine the cleanup, we took 5 h before and after an identified spike event and eliminated periods with changes in specific conductivity larger than 2 µS cm⁻¹ and completed the task with few manual adjustments.

2.3 Statistical analysis

Vachon et al. (2020) proposed a set of parameters obtained from covariance analysis of paired exO₂–exCO₂ measurements to gain insights into aquatic ecosystem metabolism. Because buffer effects may distort the exO₂–exCO₂ data clouds at low CO₂ concentrations (Diamond et al., 2025; Nguyen et al., 2025), we included the DIC timeseries and paired exO₂–DIC measurements in our analysis. We developed a workflow in Mathematica 14.3 to analyze the Danube Delta time series at monthly and daily timescales. Days with fewer than 85 data pairs were excluded from analyses, as typically 96 pairs were recorded daily at 15 min intervals. Monthly analyses usually involved more than 2800 pairs and months with less than 500 pairs were omitted. Results of the covariance calculations included the centroids which represent the mean of the exO₂ and exCO₂ data over the period of days or months and the offset of the mean from an ideal line with slope −1 for photosynthesis and respiration (Fig. 2a). Stretch and width were obtained from the Eigenvalues of the covariance matrix as the major and minor axes length of the 95 % confidence ellipse and we obtained the slope s_cov of the stretch-line from the Eigenvectors of the covariance matrix. Parameters for the 95 %-ellipses obtained from covariance analysis provide valuable insights into the intensity of photosynthesis and respiration (stretch), the ecosystem quotient ($\text{EQ}=\mathrm{1}/\text{slope}$), and the impact of lateral inflow (centroid position, offset).

To analyse monthly day-night dynamics, we calculated the average 24 h day-night cycles for each month of sensor deployment. For instance, monthly averages for 08:00 am in June were obtained by taking the mean of all June measurements at 08:00, 08:15, 08:30 and 08:45 EET. This approach resulted in average timing and amplitudes of the exO₂, exCO₂ and DIC cycles. The time-of-day averages reduced the influence of outliers and resulted in slopes s₂₄ that may reflect EQs more closely than the tilt of the monthly covariance ellipse (s_cov) (Fig. 2b). The following analysis, however, focuses on amplitudes and stretch as indicators of metabolic intensity and on the offset and centroid positions as indicators of wetland discharge.

For a more detailed analysis, boxplots of the amplitudes of exO₂, exCO₂, DIC and stretch parameters were based on all available monthly data at the seven stations (Table A1). The boxes show the median and quartiles (25 % and 75 %), whiskers indicate the range of values (excluding outliers). The statistical difference of these parameters for metabolic intensity between river and delta stations was evaluated by Kruskal–Wallis tests. To assess the influence of temperature and daylight on the aquatic ecosystem metabolism, stretch parameters from daily covariance analyses were correlated by linear regression and covariance analysis with average daily cloud cover data from Tulcea Airport and water temperature measured in-situ by the EXO2 probes.

Figure 2 (a) Example of covariance analysis of exO₂–exCO₂ data from July 2018 at Balanova station. Small dots are individual data points acquired every 15 min. The ellipse encloses 95 % of the datapoints and is defined by the centroid and the two main axes (125 and 29 µM) derived from the Eigenvalues of the covariance matrix. The offset corresponds to the nearest distance of the centroid from the grey dotted line with a slope of −1 for a theoretical photosynthesis-respiration process. (b) Using the same data as in (a), the monthly averages over 24 h of the day define a mean daily cycle for July at Balanova station. The average daily cycle shows an exCO₂ amplitude of 36 µM but a wider exO₂ amplitude of 88 µM.

3 Results

3.1 Intensity of river metabolism

The intensity of photosynthesis and respiration, as recorded by CO₂ oversaturation with respect to atmospheric equilibrium, O₂ deficits and DIC concentrations, varied strongly in space and time during the observation period. The results allow testing of the hypothesis that aquatic metabolism is more intense within the delta compared to the main river reaches. The inflowing Danube at Tulcea showed DIC minima of 2.6 mM in the warm season and maxima of 3.5 mM during the cold months (Fig. C1). Averaging the 24 h cycles for O₂, CO₂ and DIC at the monthly scale provided insights into the intensity and seasonality of river metabolism (Figs. 3 and C2). In 2017 the sensors at Tulcea recorded small and rather constant supersaturation (exCO₂) and undersaturation (negative exO₂) in the range of ±40 µM with diel amplitudes from ∼1 µM in the cold months and below 10 µM in the warm season. Within the delta, the 24 h cycles at Balanova channel showed maximum daily fluctuations of more than 100 µM for exO₂ in August. Overall, the 24-h amplitudes spanned two orders of magnitude between the extremes of 1 µM for the cold season in the upstream Tulcea station and 100 µM for warm season amplitudes within the delta.

The Balanova observations were characterized by maxima in the early morning (06:00–08:00 am EET) and minima in the afternoon (04:00–06:00 pm EET) between April and November (Fig. 3). Compared to exCO₂, the oxygen cycles represented a mirror image where the strong undersaturation in the morning hours relaxed to near-equilibrium conditions in the late afternoon. Compared to CO₂ and O₂, the 24 h DIC curves showed more stochastic deviations from the sinusoidal pattern, likely reflecting the rapidly changing DIC concentration in the inflow (Figs. C1 and C2).

Figure 3 Examples of 24 h monthly averages of exCO₂ (solid lines) and exO₂ (dashed lines) at Tulcea from March–December 2017 and at Balanova station from February–December 2018. The averaged daily cycles define characteristic amplitudes and timing of river metabolism for each month. See Fig. C2 for corresponding 24-cycles of DIC.

A comparison of median 24 h amplitudes based on four indicators of metabolic intensity (amplitudes of exO₂, exCO₂, DIC and the stretch parameter from covariance analysis) revealed significant differences between sites (Fig. 4). The two downstream sites of St. George and Chilia showed consistently higher median amplitudes compared to Tulcea. Amplitudes of exCO₂ were generally lower than those of exO₂, whereas the 24-h DIC variability was typically higher. Overall, the daily metabolic cycles intensified on the final $\sim \mathrm{60}-\mathrm{120}\mathrm{km}$ of river flow between Tulcea and the three near shore stations.

Figure 4 Box plots of monthly mean 24 h amplitudes of exO₂, exCO₂, DIC and the stretch parameter. Stations from left to right: Tulcea, Sulina and St. George (2016–2017), Chilia in 2016 and Puiu-Rosu, Busurca and Balanova in 2018. The number of monthly means per station ranged from 11–21 (Table A1). White horizontal lines mark median values, boxes show the 25 %–75 % quartiles and whiskers indicate the range except for outliers marked with dots.

Compared to the river stations, amplitude and stretch values recorded in the delta were significantly larger (p<0.01 for exO₂, exCO₂ and stretch, p∼0.01 for DIC). Overall, the average daily fluctuations revealed increasing median values in the amplitude of exO₂ from the apex to the inner delta, with low values from 6.5 and 4.5 µM at Tulcea and Sulina, respectively, to higher amplitudes of 29 and 9.3 µM at St. George and Chilia and consistently high values within the delta (28, 50, and 29 µM in Puiu-Rosu, Busurca and Balanova). In terms of median amplitudes, metabolic activity per unit volume was thus about 4–8 times higher at stations within the delta than at the upstream station of Tulcea. A similar pattern was observed for the median DIC amplitudes with 13 and 12 µM in Tulcea and Sulina, 50 and 23 in St. George and Chilia and elevated values within the delta: 23, 32 and 40 µM at Puiu-Rosu, Balanova and Busurca, respectively. Using DIC as an indicator, metabolic activity within the delta was 2–3 times higher than at the upstream Danube station. The results demonstrate how river metabolism intensified in the slow-flowing waters of the delta.

Plotting the 24 h cycles in the exO₂–exCO₂ plane allowed for a closer look at diel dynamics and the intensities of autotrophic versus heterotrophic ecosystem metabolism over the course of a year (Fig. 5). In the following, we focus on differences between the cold season (October–March) and the warm season (April–September). At Tulcea and Sulina, these plots illustrate minimal metabolic activity during cold months with exO₂ amplitudes in the range of 1–2 µM and comparatively small signals of 7–9 µM during the warm months. In contrast, Chilia and St. George, were characterized by significantly stronger daily cycles with mean warm season exO₂ values of 17 and 32 µM, respectively.

The delta stations reveal different types of dominant metabolic activities and a broad dynamic range of O₂ and CO₂. At the lake-outlet of Puiu-Rosu, autotrophic activity significantly exceeded heterotrophic metabolism with a pronounced buildup of dissolved oxygen and positive daytime exO₂ values during the warm season. In contrast, the wetland discharge defining the water characteristics at Busurca was dominated by heterotrophic signatures with extremely high exCO₂ concentrations of 300–500 µM during the warm season, while exO₂ oscillated between −150 and −250 µM and approached anoxic conditions. Balanova, showed an intermediate signature between the lake and wetland waters with warm-season exO₂-values ranging between 0 and −100 µM.

At high CO₂ supersaturation observed for instance at Busurca, Balanova and St. George, the slopes of linear regression lines often converged towards idealized values of −1 for the O₂:CO₂ stoichiometry of photosynthesis and respiration (Fig. 5). At Puiu-Rosu, however, steeper slopes were observed due to the buffer effect of the carbonate system when CO₂ values approach atmospheric equilibrium (Shangguan et al., 2025). The Busurca data followed a heuristic trajectory along a slope of −0.5, extending towards ${\mathrm{exO}}_{\mathrm{2}}\approx -\mathrm{200}\mathrm{\mu }\mathrm{M}$, and ${\mathrm{exCO}}_{\mathrm{2}}\approx +\mathrm{400}\mathrm{\mu }\mathrm{M}$.

Figure 5 Paired exO₂–exCO₂ plots of mean 24 h cycles. Dots represent monthly averages at one hour resolution and are connected by linear regression lines. Diagonals mark theoretical O₂–CO₂ stoichiometry of photosynthesis and respiration. Note extended x–y scale at Busurca.

To assess the effect of potential drivers for the intensity of river metabolism such as temperature and cloud cover, we related the daily stretch parameters with mean water temperature measured by the in-situ probes and the mean daily cloud cover observed at Tulcea Airport (Table 1). Cloud cover blocking sunlight was negatively correlated to the diel stretch parameters for the seven analysed time series. The Pearson's coefficients ρ ranged from −0.26 at Busurca to −0.46 at Chilia. On average, the cloud-cover sensitivity of the stretch factor was ∼3 times larger at the delta stations compared to the river sections. The R² values were typically <0.2, but the correlations were significant at the p<0.05 level.

Water temperature has been identified as a key variable for respiration rates in streams (Perkins et al., 2012). In this study, the correlation between the stretch parameters and mean daily water temperatures was rather strong (range $\mathrm{0.37}<\mathit{\rho }<\mathrm{0.77}$) and significant (p<0.05) for the seven stations (Table 1). On average, increasing water temperature enhanced the indicator of metabolic activity by 0.7 µM °C⁻¹ at the river stations and by 2.8 µM °C⁻¹ within the delta. For a seasonal shift in water temperature of 20 °C, the corresponding change in the stretch parameter amounted to 14 and 56 µM for typical Danube River reaches and delta channels, respectively. This factor 4 difference in the temperature-sensitivity of the daily stretch parameter between river and delta sites mirrors the metabolic response to cloud cover and is in line with the higher metabolic amplitudes in the delta.

Table 1 Linear regression slopes of daily stretch values as a function of mean cloud cover at Tulcea Airport and local water temperature. R² is the coefficient of determination.

3.2 Wetland discharge

The following analyses address our second hypothesis that lateral inputs modify the downstream O₂–CO₂ dynamics. The exO₂–exCO₂ patterns revealed significant differences in the stretch (elongation) of the 95 % covariance ellipses between the main reaches of the Danube in 2016 and the stations within the delta in 2018 (Fig. 6). In November, data from the river reaches were enclosed by rather circular shapes with minimal offset from the 1:1 diagonal. The June records, however, resulted in elongated ellipses characteristic for higher amplitudes. The data patterns within the covariance ellipses formed distinct monthly fingerprints which were only slightly distorted after moving from Tulcea to the downstream stations of Sulina and Chilia but appeared randomized at St. George.

Figure 6 Examples of covariance ellipses for exO₂–exCO₂ data pairs at 95 % confidence level for the Danube River reaches (upper panels) and stations within the delta (lower panels). Diagonal lines represent the theoretical slope of −1 for photosynthesis and respiration. Dots mark the centroids. Expanded axes for Busurca.

The delta stations showed smaller seasonal differences in stretch between June and November but marked shifts in the position within the exO₂–exCO₂ diagram. While the Puiu-Rosu data indicated intense photosynthesis with O₂ supersaturation and equilibrium-level CO₂, Busurca represented an extreme case of wetland discharge from the surrounding reed stands with CO₂-rich and O₂-depleted water masses. Balanova located downstream of Puiu-Rosu reflected a mixture of lake- and wetland water. Its monthly fingerprints resembled those of Puiu-Rosu but with a downward shift along the theoretical line of slope −1.

In the downstream Danube River reaches, wetland discharge led to increasing offsets and a shift along the heuristic slope −0.5 (Busurca, Fig. 5). Centroids were displaced towards O₂ undersaturation by −40 µM and CO₂ supersaturation of +35 to +85 µM. While offsets in June remained similar between Tulcea and Chilia, the two other downstream stations showed offset increases of 21 and 16 µM at Sulina and St. George, respectively. In November, covariance clouds remained compact and close to the $-\mathrm{1}:\mathrm{1}$ line except for St. George where the offset exceeded that of Tulcea by more than 20 µM. Overall, these examples from two contrasting months provide clear evidence of wetland discharge affecting dissolved gas dynamics between Tulcea and the stations near the Black Sea coast.

A more systematic analysis of the centroids allows tracking how the monthly means of exO₂–exCO₂ evolved over time. The results reveal distinct annual journeys at the four Danube stations in 2016 and 2017 (Fig. 7). The trajectories at Tulcea were confined to a narrow quadrant between 0 and ±40 µM for both exO₂ and exCO₂. At the downstream stations, however, the warm-season centroids were typically more depleted in O₂ and enriched in CO₂ compared to Tulcea. The trajectories often moved along a heuristic slope of −0.5 for wetland discharge (Busurca, Fig. 7). The increase in average exCO₂ between Tulcea and downstream stations reached factors as high as 4 (June 2016, Sulina) and 18 (June 2017, St. George). In contrast, the downstream depletion of exO₂ at Sulina and St. George was less pronounced than the exCO₂ loading and maximum exO₂-differences compared to Tulcea were limited to factors of ∼3 (Sulina, July 2017) and 2 (St. George, October 2016). The Chilia station also showed consistently lower exO₂ and higher exCO₂ values compared to Tulcea, but extreme additions of CO₂ rich waters were not observed in this more natural river reach which remains less affected by dredged channels.

Figure 7 Top: Annual trajectories of monthly centroids at the four Danube River stations in 2016 and 2017. Extended x-axis for St. George. Bottom: Monthly centroids at three stations within the Danube Delta in 2018. Extended x–y axes for Busurca.

Flood dynamics appear to play an important role at the stations within the delta. The largest flood peak during the three years of our study occurred in early April 2018 (Fig. 1) and was followed by a strong flood recession until mid-June. All three delta stations showed the highest CO₂ loads with the lowest O₂ levels in May 2018, when the flooded reed stands were likely exporting carbon-rich water into the adjacent channels and lakes. Higher discharge in July–August reversed the situation and moved the Busurca centroids towards photosynthesis-respiration line with slope −1. However, when exceptionally dry conditions re-established the wetland discharge in October–November, the Busurca centroids moved back towards the slope −0.5 line (Fig. 7).

3.3 CO₂ emissions

Finally, we address the third hypothesis that a boost in aquatic metabolism in the delta and additional wetland discharge will increase the average aquatic CO₂ emissions on the final stretch before a river reaches the sea. To compare CO₂ emission rates across the delta, we calculated average exCO₂ values at the seven stations. To avoid sampling bias caused by interruptions during winter, we interpolated missing months by average cold season values and obtained mean emission rates over two years for Tulcea, Sulina and St. George and over one year for the other stations (Table 2). The average exCO₂ at the three downstream stations was 54 µM corresponding to a 1.9 times higher supersaturation level than at Tulcea.

Table 2 Average exCO₂ supersaturation recorded in 2016–2017 at Tulcea, Sulina and St. George, in 2016 at Chilia and in 2018 at Puiu-Rosu, Balanova and Busurca. Warm = warm season, April–September, cold = cold season,October–March.

In 2018 the lake outflow at Puiu-Rosu showed only 31 % of the exCO₂ supersaturation observed at Tulcea in the two preceding years, whereas the wetland discharge at Busurca reached almost an order of magnitude higher exCO₂ levels. In a previous study based on discrete samples and floating chamber measurements, Maier et al. (2021) determined the median values of the gas transfer coefficients k₆₀₀ [m s⁻¹] for the different waterscapes of river reaches, lakes and channels. These results allowed the estimation of updated CO₂ emission rates based on the continuous sensor observations. For the river reaches, we assumed an emission scenario of a linear increasing between Tulcea and the three downstream stations. For comparison with an earlier study (Maier et al., 2021), the Puiu-Rosu station was used as representative estimate for lakes and the average between Balanova and Busurca was assumed to represent the mean channel emission rate (Table 3).

Despite substantially higher CO₂ concentration in channels compared to the river reaches and lakes, the waterways of the delta contribute a similar amount as lakes to total emissions because of their small surface area. According to data from global meta-analyses, the vast area of reed stands holds the potential for significant CO₂ drawdown (Li et al., 2025; Lu et al., 2017).

Table 3 CO₂ emission rates across aquatic ecosystems in the Danube Delta based on average exCO₂ values and median gas transfer rates. The 25 % and 75 % quartiles are shown in brackets. Gg=10⁹ g.

^* Data from Maier et al. (2021), ${}^{**}$ CO₂ uptake data from Li et al. (2025) and Lu et al. (2017)

4 Discussion

4.1 Variability of ecosystem metabolism

Why is it relevant to study the differences in river metabolism between lowland river reaches and wetland channels? Our study provides evidence for significantly faster metabolic processes in the slowly moving waters of the delta wetlands compared to the main river reaches of the Danube. Monthly averages of 24 h amplitudes and the stretch of exO₂–exCO₂ covariance ellipses consistently indicate a 4–8 times more intense aquatic ecosystem metabolism within the delta. This distinction is particularly relevant in the context of ongoing wetland restoration efforts across the Lower Danube (Csagoly et al., 2018), where changes in hydrological connectivity may substantially alter carbon processing, nutrient retention, and greenhouse-gas dynamics (Hemes et al., 2018; Tschikof et al., 2022).

Additional analyses of daily DIC variability support the evidence for more intense carbon processing within the delta, although the DIC patterns exhibit higher stochasticity due to hydrological fluctuations. In contrast the direct response of O₂ to photosynthesis and respiration, DIC concentrations are influenced by a broader range of interacting processes, like groundwater inputs, sediment-water exchange, and episodic flooding events. These hydrological perturbations can alter both the concentration and residence time of dissolved carbon, partially masking the metabolic signal at daily time scales.

The mean daily exO₂ amplitude as an indicator for metabolic activity covered two orders of magnitude from cold season values in the upstream station of the Danube with cold-season values around 1 µM at Tulcea to a warm season maximum of 100 µM observed at Balanova station in August 2018. The dense aquatic vegetation in the delta lakes and channels including floating reed, riparian vegetation, macrophytes and high chlorophyll concentration translate into high substrate availability for ecosystem respiration (Coops et al., 2008; Oprea et al., 2024). A hydrological model for the Delta estimated travel times across the network of channels and lakes on the order of 2–3 weeks for the Puiu-Rosu and Balanova stations (Oosterberg et al., 2000), facilitating a high biomass concentration in prolonged contact with moving water. Average exCO₂ concentrations close to 0.5 mM at Busurca in May and July 2018 may therefore result from the accumulation of mineralization products over timespans of weeks in water parcels traversing the reed stands (Fig. 5). By contrast, water residence times in the main river reaches are much shorter: sensor data revealed travel times between Tulcea and in the three downstream stations on the order of one day for Sulina and St. George and two days for Chilia essentially limiting the time for accumulating changes in O₂ and CO₂.

The significant differences between the three river reaches, can likely be attributed to their morphology and hydrological connectivity. Patterns exO₂–exCO₂ pairs within covariance ellipses retained close similarity between Tulcea and the downstream stations of Sulina and Chilia, whereas they appeared “randomized” at St. George (Fig. 6). The Sulina channel has been strongly modified for navigation with straight geometry and dykes that reduce lateral exchange, while the Chilia branch retains a more natural meandering geometry. The St. George reach consist of natural meanders and straight cuts for navigation, resulting in a broader distribution of water travel times and more intense exchange with adjacent wetlands and aquifers. This close interaction with wetlands at St. George likely superimposes multiple water masses with different metabolic histories, thereby disrupting the coherent diel covariance patterns that were preserved at Sulina and Chilia. Consequently, the 24 h amplitudes of exO₂, exCO₂ and DIC all follow the sequence St. George > Chilia > Sulina consistent with decreasing lateral connectivity (Fig. 4).

We used the daily stretch factors to test hypotheses that water temperature and cloud cover would influence the intensity of aquatic metabolism (Table 1). Overall, the correlation between cloud cover and covariance stretch was weaker than for water temperature as a forcing factor. The relations between cloud cover, photosynthetically active radiation, PAR, (Kathilankal et al., 2014) and metabolic response are inherently complex in multiyear (Dodds et al., 2013) or multisite analyses (Mulholland et al., 2001). As cloud cover does not directly translate to in-situ light conditions due factors like the seasonality of daylight and water transparency, further assessments require local PAR analyses at the monitoring stations. Such local data were available for determining the temperature sensitivity of metabolic amplitudes. The regression of daily mean in-situ temperatures with diel stretch factors resulted in a temperature sensitivity of the aquatic metabolism of 0.7 µM °C⁻¹ in river reaches and 2.8 µM °C⁻¹ in delta channels. These results suggest that a warming climate would affect the wetland's aquatic metabolism about four times more strongly than main river branches. Previous studies have shown that biofilm respiration exhibits complex temperature sensitivity (Dybdahl et al., 2024) while ecosystems with higher quantities of organic substrate are subject to stronger temperature sensitivity of ecosystem respiration (Jankowski et al., 2014).

The timing of daily cycles indicates a preferred sampling window around noon for obtaining representative values of O₂, CO₂ and DIC, whereas the early hours of 05:00–08:00 am EET and the late afternoon of 03:00–06:00 pm would introduce the largest sampling biases in the case of the lower Danube aquatic ecosystems (Figs. 3 and D1). For logistical reasons, field work is often scheduled between the late morning and early afternoon hours, which likely explains, why the CO₂ emission rates based the continuous monitoring of this study (Table 3) agrees very with results from grab samples collected during field campaigns during the same years (Maier et al., 2021).

4.2 Wetland discharge

What are the effects of wetland discharge on river biogeochemistry and how can we identify its origin and magnitude? Large-scale approaches may analyse the CO₂ dynamics within an entire river system: sampling transects across the Congo River a study demonstrated that lateral connectivity with riparian wetlands was driving spatial and temporal variability of CO₂ concentrations and the river's greenhouse gas emissions (Borges et al., 2019). More targeted study designs focus directly on wetland drainage pathways: recent studies revealed disproportionally high contributions of ditches to CO₂ emissions from Dutch peatlands (van der Knaap et al., 2025) and from a Phragmites dominated Chinese wetland (Xue et al., 2025). In the present study, we combined both approaches, although logistical constraints prevented synchronous observations of river reaches and delta channels.

The order of magnitude of exCO₂ found in the delta channels compares well with a global meta-analysis of wetland ecosystem metabolism (Richardson et al., 2024). These authors reported mean values of dissolved pCO₂ for fens and bogs of 3400 and 4100 ppm, respectively, which corresponds to approximately 100 µM exCO₂ and falls within the range of downstream Danube stations with exCO₂ of ∼50 µM and the extreme case of Busurca with 500 µM (Fig. 7). The mixing ration of wetland discharge in the Danube branches can be estimated based on the mean warm-season offset at Busurca (∼100 µM). If we interpret offsets of 0–100 µM as a proxy for lateral inflow of 0 %–100 %, then the warm-season offsets observed in 2016 at Chilia, Sulina, and St. George would correspond to contributions of wetland discharge of 13, 22 % and 25 % to total flow in these River reaches. These estimates should be considered as upper limits considering that the diversion of Danube water through the delta was estimated at 10 % (Oosterberg et al., 2000). The discrepancy can be attributed to incomplete horizontal mixing, and the positioning of the sensors close to the riverbanks. At the St. George station the inflow from the Tartaru Channel probably contributed the large offsets observed between May and July 2017 (Fig. 7).

A quantitative assessment of the different lateral interactions affecting the three Danube River reaches was based on monitoring the average monthly difference in DIC concentrations between the downstream stations and the inflow at Tulcea (Fig. C3). Integrating these differences from April–December 2016 reveals an annual DIC addition of 0.5 and 1.2 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ at Sulina and St. George, respectively. The observations in 2017 provided a similar result for St. George, but a higher accumulation rate for the Sulina stations. The seasonal balance for the Chilia branch, however, resulted in a negative value ($-\mathrm{0.1}\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$). These contrasting observations can be explained by differences in hydrological connectivity and the origin of lateral inflows. Sulina and St. George are receiving carbon-rich wetland water. The important northern inflow to Chilia, however, originates from Lake Yuluph, the largest freshwater lake in Ukraine. Water from this system is transported via Lake Kuhurlui and connecting channels to the Chilia branch and appears to deliver DIC-depleted water.

A final question remains: how can we explain the absolute stoichiometry of up to 1:2 for exO₂ : exCO₂ during the warm season at the station downstream of the delta and within the Balanova and Busurca channels receiving wetland discharge? A mapping campaign in October 2017 with in-situ CH₄ sensors reported median CH₄ concentrations of 0.54, 0.70 and 2.2 µM for the Danube reaches, lakes and channels and an extreme value of 16 µM at a hotspot site near Busurca station (Canning et al., 2021). These data provide indirect evidence for an important source of CO₂ co-generated by methanogens in anaerobic sediments of the reed stands. Although lateral water fluxes were not measured directly, the combination of CO₂ enrichment, O₂ depletion and deviations from 1:1 stoichiometry provides a consistent biogeochemical signature of wetland-derived inputs (Borges et al., 2019; Maier et al., 2021; Zuijdgeest et al., 2016). This anaerobically generated CO₂ contributes to the offset towards a diagonal of slope −0.5 (Fig. 7). A portable membrane inlet mass spectrometer for O₂, N₂, He and Ar used during the same campaign in May and October 2017 revealed additional processes affecting O₂ : CO₂ ratios including excess air formation in riverbanks and root ventilation in Phragmites stands. These processes may fuel methane oxidation but will remain “invisible” to sensors in the water column (Maier et al., 2022). Therefore, anaerobically generated CO₂ accumulating in the water column of reed beds, provided a distinct biogeochemical signature for identifying wetland discharge to channels and river reaches in the delta.

4.3 CO₂ emissions by the wetland pump

In their perspective paper, Abril and Borges (2019) argued that riparian wetlands should be included in an expanded version of the reactive pipe concept for rivers (Cole et al., 2007). With their high productivity and slow gas-exchange at the air–water interface, riparian wetlands release large amounts of dissolved CO₂. Additional fluxes of plant litter, and DOC add to the function of wetland discharge as a carbon pump. There is now increasing evidence that the wetland CO₂ pump is increasing CO₂ emissions of inland waters across different climatic zones from the tropical Congo River (Borges et al., 2019) to the arctic River Ob (Vorobyev et al., 2024).

The exCO₂ records (Fig. 5 and Table 2) allow constraining magnitude and timing of the wetland pump in the Danube Delta. Overall, the downstream stations exhibited exCO₂ values and emission rates which were 1.7–2.1 times higher than measured at Tulcea (Fig. 8). Within the delta the high productivity of the Lake Puiu-Rosu complex reduced the average CO₂ supersaturation to 30 % with reference to the inflowing Danube water indicating substantial CO₂ uptake. In contrast, wetland discharge increased the average exCO₂ by factors of 1.6 and 8.6 at Balanova and Busurca. The excess CO₂ levels remained relatively constant at Tulcea and in the near natural Chilia branch receiving mainly lake outflow. Within the delta channels and at the Sulina and St. George stations, however, the wetland pump was most active during the warm season with exCO₂ ratios between the warm and cold season of 1.3–2.1 (Table 2). This pronounced seasonality of the carbon pump is likely driven by the significant temperature dependence of the aquatic metabolism (Table 1) and the hydrological dynamics with flood peaks in spring an early summer (Fig. 1a).

Figure 8 Schematic overview of average CO₂ emission rates in main river reaches upstream and downstream of the Danube Delta for 2016-2017 compared to stations within the delta in 2018 and global median net ecosystem production (NEP) uptake by wetlands from Li et al. (2025) and Lu et al. (2017).

The day-night cycles of exCO₂ reached their maximum typically in the morning hours and showed a minim in the late afternoon and early evening (Fig. 3). Individual water samples were taken across the delta in 2016 and 2017 during daytime (Maier et al., 2021). The analysis of the three river reaches revealed no significant difference between individual water samples and the average values from the high-frequency records in this study (Table 3). This result seems to contradict a global study (Gomez-Gener et al., 2021) which reported evidence for higher night-time CO₂ emissions from rivers based on an extensive compilation of day and night CO₂ data. Analysing the 24 h cycles at the four Danube River stations in more detail (Fig. D1) reveals that sampling around midday captured values relatively close to mean daily conditions because exCO₂ maxima and minima were cantered around the early morning and late afternoon hours, respectively. For systems affected by the wetland pump, these results indicate that the day-night sampling bias will depend on factors like lateral connectivity and water exchange rates.

Finally, we must address the “elephant in the delta” – namely the role of CO₂ uptake by wetland vegetation and specifically by Phragmites stands. Direct carbon uptake by emergent wetland plants remains largely invisible to aquatic sensing, yet meta-analyses of net ecosystem production (NEP) report average global values of 168 and 208 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ for marshes (Li et al., 2025) and coastal wetlands (Lu et al., 2017), respectively. Applied to a 1600 km² area of Phragmites dominated reed stands in the Danube Delta, these rates correspond to a carbon uptake on the order of 270–330 Gg C yr⁻¹. This estimated uptake is roughly five time higher than the diffuse CO₂ emissions from water surfaces in the delta region estimated as 49 Gg C yr⁻¹ (this study, Table 3) and 60 Gg C yr⁻¹ (Maier et al., 2021). This scenario estimates indicate that the wetland CO₂ pump observed in the Danube delta is driven by the highly productive reed stands which act as strong net sink for CO₂ at the landscape scale. Much of the CO₂ released into channels and river branches originates from the respiration of recently fixed organic matter within reed stands and wetland sediments. At the same time, the sustained primary production of Phragmites seasonally removes atmospheric CO₂ and transfers carbon into submerged biomass and accumulating sediments. The high CO₂ emissions from water surfaces in the delta are therefore part of an efficient carbon processing ecosystem. A complete assessment of the carbon budget would require ecosystem-scale CO₂ and CH₄ exchange measurements, for example by eddy-covariance techniques, which have provided the basis for many wetland carbon budgets (Zou et al., 2022).

5 Conclusion

This study combined covariance analysis with averaging 24 h cycles of O₂–CO₂ pairs from sensor measurements at seven locations in the Danube Delta region. The dataset covered 105 monthly and about 3000 daily cycles. The results provide answers to the guiding questions outlined in the introduction:

• How does the intensity of river metabolism differ between the open waters of the delta compared to the upstream Danube River? Analysing the intensity of exO₂–exCO₂ dynamics revealed two orders of magnitude difference in the amplitude of average diel cycles among stations and seasons with the lowest values during the cold season in the upstream stretch of Danube and the highest intensity in August at the the Balanova channel. The monthly average of 24 h amplitudes and the stretch of exO₂–exCO₂ covariance ellipses consistently indicated a 4–8 times more intense aquatic ecosystem metabolism within the delta compared to the main river branches.

• How significant is the contribution of lateral inflows from lake systems and reedbeds to the O₂ and CO₂ dynamics of the main river reaches crossing the delta? Wetland discharge could be identified by CO₂ rich water causing an offset with respect to the $-\mathrm{1}:\mathrm{1}$ exO₂ : exCO₂ stoichiometry. Based on this indicator local wetland discharge was estimated as 13 %–25 % of total discharge. These numbers exceed the total estimates of river diversion through the delta and should therefore be considered an upper limit. Balancing DIC increases along the final stretches of the main Danube branches revealed clear differences between low-carbon inflow from lakes in Ukraine and high-carbon wetland discharge in the southern part of the Delta.

• What are the overall effects of river metabolism and wetland discharge on aquatic CO₂ emission rates in the Danube Delta? The exCO₂ records helped quantifying the magnitude of the wetland CO₂ pump in the Danube Delta. Overall, the downstream stations near the Black Sea exhibited exCO₂ values which were 1.7–2.1 times higher than at Tulcea. The estimated annual diffuse CO₂ emissions from water surface in the delta of 49 Gg C yr⁻¹ were compatible with an earlier estimate based on more stations but less frequent sampling of 60 Gg C yr⁻¹. Estimates of the landscape-level CO₂ sink by the highly productive wetlands, however, are ∼5 times higher.

In summary, the well-connected delta wetlands provide a boost in ecosystem metabolism which leaves clear signals of the wetland CO₂ pump in the lateral discharge reaching the main Danube branches. At the landscape scale, however, there is evidence from the literature that the highly productive wetland areas act as a net carbon sink which retains organic carbon in the sediments and exports POC and DOC to the Black Sea (Durisch-Kaiser et al., 2010; Maier et al., 2021).

Appendix A: Monitoring stations

The monitoring stations were located upstream of the Delta at Tulcea and downstream along the three branches of the Danube at St. George, Sulina and Chilia (Fig. 1). During the final year, the sensors were relocated to the Puiu-Rosu, Balanova, and Busurca channels within the Delta. The time coverage during the field campaign ranged between 11 and 21 months per station (Table A1). Spike removal at Sulina and Busurca reduced the number of observation days (Sect. 2.2). Due to logistical constraints the sensors were typically moored from anchored boats (Fig. A1).

Table A1 Time series of O₂, CO₂ and DIC observations.

^* Months with >500 data, ${}^{**}$ days with >85 data.

Figure A1 Characteristics of monitoring sites. (a) The Chilia site was located on the Romanian side of a small branch that defines the border with Ukraine. Despite its location in a sidearm, the water showed the same chemical signature as the main branch at the time of installation. (b) EXO2 probe at the Chilia station during winter conditions. (c) Clear water from the Busurca channel entering the Sulina reach of the Danube from the right. The Busurca monitoring station was located 25 m upstream of the confluence with the Sulina branch. (d) Sulina station located about 1 km downstream of Busurca junction and moored at an old, anchored ship. (e) View from Tartaru channel close to Balanova station. The Danube branch of St. George in the background with its turbid water is flowing from left to right towards the Black Sea. (f) For logistical reasons, the monitoring station at St. George branch was located about 500 m downstream of the confluence with the Tartaru channel. The Black Sea is visible on the horizon. (g) The Balanova multiprobe was deployed in the Central Channel connecting Sulina and St. George branches. The channel runs parallel to the gravel road linking two communities. The location is influenced by water originating from Lake Rosu. Balanova is a local name. (h) The Puiu-Rosu sensor recorded water quality at the outlet of Lake Puiu. The station was located on a small channel connecting Lake Puiu with Lake Rosu. (i) The Busurca multiprobe was moored on a small, anchored houseboat. Not shown: The Tulcea station monitoring Danube water at the apex of the delta. The station was located 300 m downstream of the bifurcation between the Tulcea and Chilia branches.

Appendix B: Calculation of CO₂

The EXO2 sensors measured pH, temperature and specific conductivity at a frequency of 15 min. Using grab samples we correlated the continuous conductivity data with discrete lab-based alkalinity measurements (Raymond et al., 2013) and used the COSYS code (van Heuven et al., 2011) to calculate dissolved CO₂ concentrations (Sect. 2.2).

Figure B1 (a) Linear regression of alkalinity analyses from water samples versus specific conductivity. The specific conductivity recorded by the EXO2 sensors was correlated with alkalinity measurements obtained from titration analysis of water samples taken at the same time and location. Labels mark the four different EXO-2 probes. Alkalinity $\left[\mathrm{mM}\right]=\mathrm{0.0057}\times$ specific conductivity [µS cm⁻¹] +0.60, with a correlation coefficient of R²=0.70. (b) Measured CO₂ data based on the headspace technique compared to dissolved CO₂ calculated from EXO2 data Alkalinity was estimated by correlation (panel a) and pH plus temperature records were combined to calculate dissolved CO₂. With few exceptions, the CO₂ estimates compared well with direct CO₂ measurements (R²=0.82). Diagonal marks the 1:1 relationship with lines at ±1000 ppm. Symbols represent the same sensors as in panel a).

Appendix C: DIC dynamics at Danube Delta stations

Figure C1 Timeseries of DIC [mM] at the seven stations from January 2016–December 2018. DIC based on measurements of conductivity, temperature, pH, and the calibration with titration alkalinity (Fig. B1a).

Figure C2 Example of 24 h monthly averages for DIC at Tulcea from March–December 2017 and at Balanova station from February–December 2018. The DIC [micromolar] cycles show stronger random variability compared to the exO₂ and exCO₂ cycles in Fig. 3. Rapid changes in DIC at monthly timescales are likely and important source of this stochasticity (Fig. C1).

Figure C3 Cumulative increase in DIC concentration during seasonal cycles calculated from the difference between the mean monthly DIC at a station and the upstream data at Tulcea.

Appendix D: Overview of 24 h cycles

Figure D1 Monthly averaged time-of-day plots showing exCO₂ (solid lines) and exO₂ (dashed) during 2016 and 2017 at the Danube stations during 2018 at the Delta stations. Extended y-axis for Busurca.