Methane, carbon dioxide, and nitrous oxide emissions from two clear-water and two turbid-water urban ponds in Brussels (Belgium)

Shallow ponds can occur either in a clear-water state dominated by macrophytes or a turbid-water state dominated by phytoplankton, but it is unclear if and how these two alternative states affect the emission of greenhouse gases (GHGs) such as carbon dioxide (CO₂), methane (CH₄), and nitrous oxide (N₂O) to the atmosphere. We measured the dissolved concentration of CO₂, CH₄, and N₂O from which the diffusive air–water fluxes were computed, in four urban ponds in the city of Brussels (Belgium): two clear-water macrophyte-dominated ponds (Silex and Tenreuken), and two turbid-water phytoplankton-dominated ponds (Leybeek and Pêcheries) on 46 occasions over 2.5 years (between June 2021 and December 2023). Ebullitive CH₄ fluxes were measured with bubble traps in the four ponds during deployments in spring, summer, and autumn, totalling 48 d of measurements. Measured ancillary variables included water temperature, oxygen saturation level ( %O₂), concentrations of chlorophyll-a (Chl-a), total suspended matter (TSM), soluble reactive phosphorus (SRP), nitrite (NO${}_{\mathrm{2}}^{-}$), nitrate (NO${}_{\mathrm{3}}^{-}$), and ammonium (NH${}_{\mathrm{4}}^{+}$). The turbid-water and clear-water ponds did not differ significantly in terms of diffusive emissions of CO₂ and N₂O. Clear-water ponds exhibited higher values of ebullitive CH₄ emissions compared to turbid-water ponds, most probably in relation to the delivery of organic matter from macrophytes to sediments, but the diffusive CH₄ emissions were not significantly different between clear- and turbid-water ponds. Across seasons, CH₄ emissions increased with water temperature in all four ponds, with ebullitive CH₄ fluxes having a stronger dependence on water temperature (Q₁₀) than diffusive CH₄ fluxes. The temperature sensitivity of ebullitive CH₄ fluxes decreased with increasing water depth, implying that shallow sediments would respond more strongly to warming (e.g. heat waves). Total annual CH₄ emissions (diffusive + ebullitive) in CO₂ equivalents equalled those of CO₂ in turbid-water ponds and exceeded those of CO₂ in clear-water ponds, while N₂O emissions were negligible compared to the other two GHGs. Total annual GHG emissions in CO₂ equivalents from all four ponds increased from 2022 to 2023 due to higher CO₂ diffusive fluxes, likely driven by higher annual precipitation in 2023 compared to 2022 (leading putatively to higher inputs for organic or inorganic carbon from run-off), possibly in response to the intense El Niño event of 2023. The findings of this work suggest that it might be necessary to account for the presence of submerged macrophytes when extrapolating ebullitive CH₄ fluxes in ponds at a larger scale (regional or global) (particularly if Chl-a is used as a descriptor), although it might be less critical for the extrapolation of diffusive CH₄, CO₂, and N₂O fluxes.

1 Introduction

Greenhouse gas (GHG) emissions from inland water (rivers, lakes, and reservoirs) to the atmosphere such as carbon dioxide (CO₂), methane (CH₄), and nitrous oxide (N₂O) are quantitatively important for global budgets (Lauerwald et al., 2023). Global GHG emissions from lakes are lower than from rivers for CO₂ (Raymond et al., 2013) and for N₂O (Lauerwald et al., 2019; Maavara et al., 2019). However, reported global emissions of CH₄ from lakes (Rosentreter et al., 2021; Johnson et al., 2022) are equivalent or even higher compared to rivers (Stanley et al., 2016; Rocher-Ros et al., 2023). Global emissions of CO₂ and CH₄ from lakes to the atmosphere represent 1.25 to 2.30 Pg CO₂ equivalents (CO₂-eq) annually, with a significant proportion from CH₄ emissions, and represent nearly 20 % of global CO₂ emissions from fossil fuels (DelSontro et al., 2018). The contribution of CO₂ and CH₄ emissions from small lentic water bodies (small lakes and ponds) can be disproportionately high compared to large systems (Holgerson and Raymond, 2016) as small lakes and ponds are the most abundant of all water body types in number (Verpoorter et al., 2014; Cael et al., 2017), and fluxes per square metre are usually higher from smaller water bodies. The emissions of GHGs from artificial water bodies such as agricultural reservoirs, urban ponds, and storm-water retention basins could be higher than those from natural systems (Martinez-Cruz et al., 2017; Grinham et al., 2018; Herrero Ortega et al., 2019; Gorsky et al., 2019; Ollivier et al., 2019; Peacock et al., 2019, 2021; Webb et al., 2019; Bauduin et al., 2024a). These higher GHG emissions seem to result from higher external inputs of anthropogenic organic carbon and dissolved inorganic nitrogen (DIN) into artificial systems but might also reflect other differences compared to natural systems such as in hydrology (Clifford and Heffernan, 2018). Among artificial systems, urban ponds are the subject of a growing body of literature on GHG emissions (Singh et al., 2000; Natchimuthu et al., 2014; van Bergen et al., 2019; Audet et al., 2020; Peacock et al., 2021; Goeckner et al., 2022; Ray and Holgerson, 2023; Ray et al. 2023; Bauduin et al., 2024a). Urban areas can have numerous small artificial water bodies mostly associated with green spaces such as public parks, and their number is increasing due to rapid urbanization worldwide (Brans et al., 2018; Audet et al., 2020; Gorsky et al., 2024; Rabey and Cotner, 2024). Urban ponds are generally small and shallow, and usually their catchment consists of a majority of impervious surfaces with a smaller contribution from soils (Davidson et al., 2015; Peacock et al., 2021). In general, the main function of urban ponds is for storm-water management, but they provide additional benefits including aesthetic/recreational amenities and habitats for wildlife (e.g. Tixier et al., 2011; Hassall, 2014).

Shallow ponds and lakes occur in two alternative states corresponding to systems with either clear waters (macrophyte-dominated) or turbid waters (phytoplankton-dominated) during the productive period of the year (spring and summer in mid-latitudes) (Scheffer et al., 1993). Submerged macrophytes and phytoplankton regulate CO₂ dynamics directly through photosynthesis that can be more or less balanced by community respiration in the water column (e.g. Sand-Jensen and Staehr, 2007). However, it is not clear whether the presence of macrophytes increases or decreases the net CO₂ emissions from ponds and lakes. Some studies have shown a decrease in CO₂ emissions with increasing macrophyte density (Kosten et al., 2010; Ojala et al., 2011; Davidson et al., 2015), but other studies showed the opposite pattern (Theus et al., 2023). In phytoplankton-dominated lakes, CO₂ concentrations depend in part on the developmental stage of phytoplankton, with the growth and peak phases generally coinciding with lower CO₂ concentrations due to photosynthesis (Grasset et al., 2020; Vachon et al., 2020).

CH₄ emissions have been reported to increase with the concentration of chlorophyll-a (Chl-a) in phytoplankton-dominated lakes (DelSontro et al., 2018; Borges et al., 2022). The presence of macrophytes strongly affects the production of CH₄ in freshwaters (Bastviken et al., 2023), and vegetated littoral zones of lakes exhibit higher CH₄ emissions than non-vegetated zones (Hyvönen et al., 1998; Huttunen et al., 2003; Juutinen et al., 2003; Desrosiers et al., 2022; Theus et al., 2023). Macrophytes influence organic matter decomposition processes in sediments depending on the quality and quantity of plant matter they release into their environment (Reitsema et al., 2018; Grasset et al., 2019; Harpenslager et al., 2022; Theus et al., 2023). Yet, few studies have consistently compared CH₄ emissions in clear-water and turbid-water ponds (Hilt et al., 2017). A study in Argentina reported higher dissolved CH₄ concentrations in clear-water ponds with submerged macrophytes compared to turbid-water phytoplankton-dominated ponds but no differences in measured CH₄ emissions (Baliña et al., 2023).

The production of N₂O predominantly occurs through microbial nitrification and denitrification that depend on DIN, O₂ levels, and temperature (Codispoti and Christensen, 1985; Mengis et al., 1997; Velthuis and Veraart, 2022). Competition for DIN between primary producers and N₂O-producing microorganisms can impact N₂O production. Additionally, the transfer of labile phytoplankton organic matter to sediments fuels benthic denitrification and impacts N₂O fluxes. Eutrophication is assumed to drive high N₂O emissions from lakes and ponds (Audet et al., 2020; Webb et al., 2021; Wang et al., 2021; Xie et al., 2024), but some lakes with elevated Chl-a concentrations can act as sinks of N₂O due to removal of N₂O by denitrification (Webb et al., 2019; Borges et al., 2022, 2023). The presence of macrophytes also strongly influences nitrogen cycling in sediments of lakes and ponds (Barko et al., 1991; Choudhury et al., 2018; Deng et al., 2020; Dan et al., 2021) and should in theory also affect N₂O emissions, although this is seldom investigated, and available studies provide contradictory conclusions. N₂O emissions have been shown to follow diurnal cycles of O₂ concentrations in areas dominated by submerged macrophytes in Lake Wuliangsuhai (China) (Ni et al., 2022) and the seasonal cycle of aboveground biomass of emerged macrophytes (Phragmites) in Baiyangdian Lake (China) (Yang et al., 2012). On the contrary, a study showed there was no significant difference of N₂O production in sediments of macrophyte-rich (n=10) and macrophyte-free (n=12) lakes in subtropical China (Liu et al., 2018). There has been a very limited number of studies systematically investigating how emissions differ between ponds dominated by phytoplankton and those dominated by macrophytes (Baliña et al., 2023) and none simultaneously investigating CO₂, CH₄, and N₂O emissions including both diffusive and ebullitive components.

The emissions of CO₂ and N₂O from aquatic systems are exclusively through diffusion across the air–water interface (diffusive flux), while CH₄ can be additionally emitted as bubbles are released from sediments to the atmosphere (ebullitive flux). At annual scale, ebullitive CH₄ flux usually represents more than half of total (diffusive + ebullitive) CH₄ emissions from shallow lakes (Wik et al., 2013; Deemer and Holgerson, 2021), although the relative contribution of ebullitive and diffusive CH₄ emissions is highly variable seasonally (e.g. Wik et al., 2013; Ray and Holgerson, 2023; Rabaey and Cotner 2024). Ebullitive CH₄ fluxes are particularly high in the littoral zone of lakes at depths <5 m (Wik et al., 2013; DelSontro et al., 2016; Borges et al., 2022) and strongly increase in response to temperature (DelSontro et al., 2016; Aben et al., 2017; Rabaey and Cotner, 2024), as well as organic matter availability (DelSontro et al., 2016, 2018). Ebullitive CH₄ fluxes tend to be higher in small and shallow water bodies (Deemer and Holgerson, 2021) but are notoriously variable in time and space and are difficult to estimate reliably (DelSontro et al., 2011).

Here, we report a data set of CO₂, CH₄, and N₂O dissolved concentrations in four shallow and small urban ponds (Leybeek, Pêcheries, Silex, and Tenreuken) in the city of Brussels (Belgium) (Fig. 1), with data collected 46 times at regular intervals (between June 2021 and December 2023) on each pond. The air–water diffusive fluxes of CO₂, CH₄, and N₂O were calculated from dissolved concentrations and the gas transfer velocity, while the ebullitive CH₄ fluxes were measured with inverted funnels during eight deployments (totalling 48 d) in the four ponds. The four ponds have similar depth, surface area, and catchment urban coverage and mainly differ by the phytoplankton-macrophyte dominance, a clear-water state dominated by macrophytes and a turbid-water state dominated by phytoplankton (alternative states) (Fig. 1). We assess whether the differences in terms of (i) CO₂, CH₄, and N₂O dissolved concentration and diffusive emissions; (ii) ebullitive CH₄ emissions; and (iii) the relative contribution of CO₂, CH₄, and N₂O to the total GHG emissions in CO₂-eq between the four ponds are explained by the two alternative states.

Figure 1 Location of the four sampled urban ponds (black diamonds) in the city of Brussels (Belgium) delineated by the black line. Right panels indicate, for each pond, the shape of the ponds, surface area (ha), perimeter (m), average depth (m), mean ± standard deviation of chlorophyll-a (Chl-a, in µg L⁻¹), and total suspended matter (TSM, in mg L⁻¹) in summer (21 June to 21 September 2021, 2022, 2023) and summer total macrophyte cover (MC; in %) (Table S1).

2 Material and methods

2.1 Field sampling and meteorological data

Sampling was carried out at a single fixed station (pontoon) in each of the four ponds, on the same day between 09:00 am and 11:00 am CET, 46 times on each pond between June 2021 and December 2023 at a frequency ranging from one (winter) to three (summer) times per month. Water was sampled 5 cm below the surface with 60 mL polypropylene syringes for analysis of dissolved concentrations of CO₂, CH₄, and N₂O. Samples for CH₄ and N₂O were transferred from the syringes with a silicone tube into 60 mL borosilicate serum bottles (Wheaton), preserved with 200 µL of a saturated solution of HgCl₂, sealed with a butyl stopper and crimped with an aluminium cap, without a headspace, and stored at ambient temperature in the dark prior to analysis in the laboratory. The partial pressure of CO₂ (pCO₂) was measured directly in the field, within 5 min of sample collection, with a Li-Cor Li-840 infrared gas analyser (IRGA) based on the headspace technique with four polypropylene syringes (Borges et al., 2019). A volume of 30 mL of sample water was equilibrated with 30 mL of atmospheric air within the syringe by shaking vigorously for 5 min. The headspace of each syringe was then sequentially injected into the IRGA, and a fifth syringe was used to measure atmospheric CO₂. The final pCO₂ value was computed taking into account the partitioning of CO₂ between water and the headspace, as well as equilibrium with HCO${}_{\mathrm{3}}^{-}$ (Dickson et al., 2007) using water temperature measured in situ and after equilibration and total alkalinity (data not shown). Samples for total alkalinity were conditioned, stored, and analysed as described by Borges et al. (2019). The IRGA was calibrated in the laboratory with ultrapure N₂ and a suite of gas standards (Air Liquide Belgium) with CO₂ mixing ratios of 388, 813, 3788, and 8300 ppm. The precision of pCO₂ measurements was ±2.0 %. Water temperature, specific conductivity, and oxygen saturation level (%O₂) were measured in situ with VWR MU 6100H probe 5 cm below the surface. A 2 L polyethylene water container was filled with surface water for conditioning the samples for other variables at the laboratory in Université Libre de Bruxelles.

Three bubble traps were deployed 50 cm apart for measuring ebullitive CH₄ flux. The bubble traps consisted of inverted polypropylene funnels (diameter 23.5 cm) mounted with 60 mL polypropylene syringes, with three-way stop valves allowing the gas to be collected without contamination from ambient air. The polypropylene funnel was attached with steel rods to a polystyrene float. The volume of gas collected in the funnels was sampled with graduated polypropylene 60 mL syringes every 24 h. The value of the collected volume of gas was logged, and the gas was transferred immediately after collection to pre-evacuated 12 mL vials (Exetainers, Labco, UK) that were stored at ambient temperature in the dark prior to the analysis of CH₄ concentration in the laboratory. The time series of measurements was longer at the Silex pond than the other three ponds.

Surveys to identify and quantify visually the relative coverage of emerged and submerged macrophytes were conducted in summer 2023 (Table S1). The resulting list of macrophyte species agreed with past studies in Brussels' ponds (Peretyatko et al., 2007). Air temperature, precipitation, wind speed, and atmospheric pressure were retrieved from https://wow.meteo.be/en (last access: 18 January 2024) for the meteorological station of the Royal Meteorological Institute of St-Lambert (50.8408° N, 4.4234° E) in Brussels, located between 2.5 and 5.0 km from the surveyed ponds. Air temperature, wind speed, and atmospheric pressure were averaged over 24 h to obtain a daily mean value. Precipitation was integrated each day to obtain cumulated daily value.

2.2 Laboratory analysis

2.2.1 Chlorophyll-a, total suspended matter, and dissolved inorganic nutrients

Water was filtered through Whatman GF/F glass microfiber filters (porosity 0.7 µm) with a diameter of 47 mm for total suspended matter (TSM) and Chl-a. Filters for TSM were dried in an oven at 50 °C, and filters for Chl-a were kept frozen (−20 °C). The weight of each filter was determined before and after filtration of a known volume of water using an Explorer™ Pro EP214C analytical microbalance (accuracy ±0.1 mg) for determination of TSM concentration. Chl-a concentration was measured on extracts with 90 % acetone by fluorimetry (Kontron model SFM 25) (Yentsch and Menzel, 1963) with a limit of detection of 0.01 µg L⁻¹. Filtered water was stored frozen (−20 °C) in 50 mL polypropylene bottles for analysis of dissolved nutrients. Soluble reactive phosphorus (SRP) was determined by the ammonium molybdate, ascorbic acid, and potassium antimony tartrate staining method (Koroleff, 1983), with a limit of detection of 0.1 µmol L⁻¹. Ammonium (NH${}_{\mathrm{4}}^{+}$) was determined by the nitroprusside–hypochlorite–phenol staining method (Grasshoff and Johannsen, 1972), with a limit of detection of 0.05 µmol L⁻¹. Nitrite (NO${}_{\mathrm{2}}^{-}$) and nitrate (NO${}_{\mathrm{3}}^{-}$) were determined before and after reduction of NO${}_{\mathrm{3}}^{-}$ to NO${}_{\mathrm{2}}^{-}$ by a cadmium–copper column, using the Griess acid reagent staining method (Grasshoff et al., 2009), with a detection limit of 0.01 and 0.1 µmol L⁻¹, respectively. Concentration of dissolved inorganic nitrogen (DIN) was calculated as the sum of NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{2}}^{-}$, and NO${}_{\mathrm{3}}^{-}$ concentrations (in µmol L⁻¹).

2.2.2 CH₄ and N₂O measurements by gas chromatography

Measurements of N₂O and CH₄ concentrations dissolved in water were made with the headspace technique (Weiss, 1981) with a headspace volume of 20 mL of ultra-pure N₂ (Air Liquid Belgium) and a gas chromatograph (GC) (SRI 8610C) with a flame ionization detector for CH₄ and an electron capture detector for N₂O calibrated with ${\mathrm{CH}}_{\mathrm{4}}:{\mathrm{N}}_{\mathrm{2}}\mathrm{O}:{\mathrm{N}}_{\mathrm{2}}$ gas mixtures (Air Liquide Belgium) with mixing ratios of 1, 10, and 30 ppm for CH₄ and 0.2, 2.0, and 6.0 ppm for N₂O. The precision of measurement based on duplicate samples was ±3.9 % for CH₄ and ±3.2 % for N₂O. Measurements of CH₄ concentration in the gas samples from bubble traps were also made by GC with the same set-up and calibration as for the determination of the dissolved concentrations in water samples.

The CO₂ concentration is expressed as partial pressure (pCO₂) in parts per million (ppm) and CH₄ as dissolved concentration (nmol L⁻¹), as frequently used in the topical literature. CH₄ concentrations were systematically and distinctly above saturation level (2–3 nmol L⁻¹), and pCO₂ values were below saturation only five times out of the 187 measurements. The N₂O concentrations fluctuated around atmospheric equilibrium, so data are presented as a percentage of saturation level (%N₂O, where atmospheric equilibrium corresponds to 100 %). The equilibrium with atmosphere for N₂O was calculated from the average air mixing ratios of N₂O provided by the Global Monitoring Division (GMD) of the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) (Dutton et al., 2023) and using Henry's constant given by Weiss and Price (1980).

2.3 Calculations

2.3.1 Diffusive GHG emissions

The diffusive air–water CO₂, CH₄, or N₂O fluxes (F_G) were computed according to

$$\begin{array}{}\text{(1)}& F_{\mathrm{G}}=k\times \mathrm{\Delta }\left[\mathrm{G}\right],\end{array}$$

where k is the gas transfer velocity, and Δ[G] is the air–water gas concentration gradient.

The atmospheric pCO₂ was measured in the field with the Li-Cor Li-840. For CH₄, the global average present-day atmospheric mixing ratio of 1.9 ppm was used (Lan et al., 2024). k was computed from a value normalized to a Schmidt number of 600 (k₆₀₀) and from the Schmidt number of CO₂, CH₄, and N₂O in freshwater according to the algorithms as a function of water temperature given by Wanninkhof (1992). k₆₀₀ was calculated from the parameterization as a function of wind speed of Cole and Caraco (1998). CH₄ and N₂O emissions were converted into CO₂ equivalents (CO₂-eq) considering a 100-year time frame, using global warming potentials of 32 and 298 for CH₄ and N₂O, respectively (Myrhe et al., 2013).

2.3.2 Ebullitive flux

Bubble flux (mL m⁻² d⁻¹) was calculated according to

$$\begin{array}{}\text{(2)}& F_{\mathrm{bubble}}={\displaystyle \frac{V_{\mathrm{g}}}{A\times \mathrm{\Delta }t}},\end{array}$$

where V_g is the volume of gas collected in the inverted funnel (mL), A is the cross-sectional surface area of the funnel (m²), and Δt is the collection time (d).

A multiple linear regression model of F_bubble dependent on water temperature (T_w in °C) and drops of atmospheric pressure (Δp in atm) was fitted to the data according to

$$\begin{array}{}\text{(3)}& {\log }_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathit{\alpha }\times T_{\mathrm{w}}+\mathit{\beta }\times \mathrm{\Delta }p+\mathit{\gamma },\end{array}$$

where γ is the y intercept, and α and β are the slope coefficients of the multiple linear regression model.

Δp was calculated according to Zhao et al. (2021):

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }p=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}{\int }_{\mathrm{0}}^{t}p-p_{\mathrm{0}};\forall p<p_{\mathrm{0}},\end{array}$$

where p is the atmospheric pressure (atm), p₀ a threshold value fixed at 1 atm, and Δt the time interval between two measurements (d) (Fig. S1).

A linear regression model of F_bubble dependent on T_w alone was fitted to the data according to

$$\begin{array}{}\text{(5)}& {\log }_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)={\mathit{\alpha }}^{\prime }\times T_{\mathrm{w}}+{\mathit{\gamma }}^{\prime },\end{array}$$

where α^′ is the slope coefficient, and γ^′ is the y intercept.

To evaluate the relative importance of T_w and Δp combined and T_w alone in driving F_bubble, the values modelled based on Eqs. (3) and (5) were compared to the observations in the Silex pond alone or in all four ponds together, for three T_w ranges (T_w<15 °C, T_w>15 °C, and the full T_w range).

Ebullitive CH₄ fluxes ($E_{{\mathrm{CH}}_{\mathrm{4}}}$ in mmol m⁻² d⁻¹) were calculated according to

$$\begin{array}{}\text{(6)}& E_{{\mathrm{CH}}_{\mathrm{4}}}=\left[{\mathrm{CH}}_{\mathrm{4}}\right]\times F_{\mathrm{bubble}},\end{array}$$

where [CH₄] is the measured CH₄ concentration in bubbles (mmol mL⁻¹).

The CH₄ content in bubbles expressed as a percentage of total gas (%CH₄) was fitted with a linear regression model dependent on T_w in the Silex pond. The correlation of %CH₄ with F_bubble was tested on the merged data of all the four ponds.

A linear regression model of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ dependent on T_w was fitted to the data according to

$$\begin{array}{}\text{(7)}& {\log }_{\mathrm{10}}\left(E_{{\mathrm{CH}}_{\mathrm{4}}}\right)=b\times T_{\mathrm{w}}+c,\end{array}$$

where b is the slope, and c is the y intercept of the linear regression.

Equation (7) is used to predict $E_{{\mathrm{CH}}_{\mathrm{4}}}$ in each pond from measured time series of T_w, allowing matching to each diffusive CH₄ flux estimate derived from Eq. (1).

The ratio of ebullitive CH₄ flux to total (diffusive + ebullitive) CH₄ flux $\left(\frac{\mathrm{Ebul}}{\mathrm{Tot}}\right)$ was fitted as a function of T_w according to DelSontro et al. (2016):

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{Ebul}}{\mathrm{Tot}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+f\times e^{\mathrm{g}\cdot T_{\mathrm{w}}}}}.\end{array}$$

The methane ebullition Q₁₀ represents the change in $E_{{\mathrm{CH}}_{\mathrm{4}}}$ per 10 °C change in T_w and was computed according to DelSontro et al. (2016):

$$\begin{array}{}\text{(9)}& Q_{\mathrm{10}}={\mathrm{10}}^{\mathrm{10}b},\end{array}$$

where b is given by Eq. (7).

The Q₁₀ of diffusive CH₄ fluxes was also computed from Eq. (9) but using the b value derived by replacing $E_{{\mathrm{CH}}_{\mathrm{4}}}$ by the diffusive CH₄ flux in Eq. (7).

2.3.3 Statistical analysis

Generalized linear mixed models (GLMMs) were used (1) to find differences in variables among ponds and (2) to relate GHG variables to their putative controls across all ponds and within individual ponds. GLMMs were computed with the lme4 package (Bates et al., 2015) in R version 4.4.1 (R Core Team, 2021).

GLMMs allowed us to compare pCO₂, dissolved CH₄ concentration, %N₂O, F_bubble, %CH₄ in bubbles, $E_{{\mathrm{CH}}_{\mathrm{4}}}$, and diffusive CH₄ fluxes among the four ponds, using the sampling date as a random effect, and post hoc tests were performed using estimated marginal means (emmeans package in R) to assess pairwise differences between ponds. This analysis aimed at investigating if patterns in GHG concentrations and emissions differed among ponds, in particular with regards to clear-water and turbid-water states. Data were also compared among the four ponds separated by seasons, but GLMMs did not converge due to an insufficient number of data points. Comparisons on log₁₀-transformed data were then made using repeated measures analysis of variance (ANOVA) (ezANOVA package in R) with Tukey's honestly significant difference (HSD) post hoc tests. This analysis aimed at investigating if patterns in data shown by the data analysis with the full data set (four seasons merged) were also observed when analysing the data separated by seasons.

For the data sets covering the whole sampling period (four seasons merged), GLMMs were constructed for pCO₂, dissolved CH₄ concentration, %N₂O, F_bubble, %CH₄ in bubbles, $E_{{\mathrm{CH}}_{\mathrm{4}}}$, and diffusive CH₄ fluxes that included T_w, precipitation, %O₂, Chl-a, TSM, DIN, and SRP as fixed effects and “pond” and “sampling date” as a random effect to account for repeated measurements with the lme4 package in R. This analysis aimed at investigating (1) the impact of photosynthesis–respiration on CO₂ concentrations and emissions based on the relationships with Chl-a, DIN, and SRP; (2) the impact of the response of methanogenesis to warming on CH₄ concentrations and diffusive and ebullitive emissions based on the relationships to T_w; and (3) the impact on DIN availability and T_w on N₂O concentrations and diffusive emissions. A linear regression model was used to assess the relationship between Chl-a, TSM, %O₂, SRP, and DIN versus T_w. This analysis aimed at investigating if some of the patterns between GHG variables versus Chl-a, TSM, %O₂, SRP, and DIN might in fact have reflected a relation between these variables and T_w indirectly. The relationships between the annual means of CH₄, CO₂ and N₂O fluxes and the annual means of a subset of variables (Chl-a, macrophyte cover, pond surface area, and depth) were assessed with linear or quadratic regressions. The modelled F_bubble values in the Silex pond and in the four ponds were compared to measured F_bubble values with a linear regression in order to evaluate the model performance. A correlation analysis was used to assess the relationship between %CH₄ and F_bubble, as well as between anomalies in annual air temperature and annual precipitation.

The linear regression model and the Pearson correlation coefficient (r) were computed in R using the functions lm and cor(method = “pearson”), respectively. Statistical significance was set at p<0.05 for all analyses. Significant differences between groups presented in box plots are indicated by lower-case letters in the figures.

3 Results

3.1 Seasonal variations of meteorological conditions and GHG concentrations

The city of Brussels experiences a temperate climate with mild weather year-round and evenly distributed abundant precipitation, totalling 837 mm on average annually for the reference period 1991–2020. The average annual air temperature was 11 °C, with a summer average of 17.9 °C and a winter average of 4.1 °C for the reference period 1991–2020. During the sampling period, from June 2021 to December 2023, T_w in the surface of the four sampled ponds (Leybeek, Pêcheries, Silex, and Tenreuken; Fig. 1) tracked the air temperature closely, which ranged between −1.5 and 30.0 °C following the typical seasonal cycle at mid-latitudes in the Northern Hemisphere (Fig. S2). The years 2022 and 2023 were about 1 °C warmer than the average for the period 1991–2020 (11 °C), while the year 2021 was closer to the long-term average (Fig. 2). The year 2022 was warmer and drier than 2021 and 2023 (Fig. 2), with positive air temperature anomalies observed evenly throughout the year (9 months out of 12) and negative precipitation anomalies in summer, autumn, and early winter (Fig. S2). The year 2021 had warmer and drier months in June and September and colder and wetter months in July and August and was wetter and colder than 2022 overall (Fig. 2). The year 2023 was marked by both positive air temperature and precipitation anomalies (Fig. S2), resulting in a wetter and warmer year than normal and compared to 2021 and 2022 (Fig. 2). Daily wind speed was generally low (<1 m s⁻¹) except for a windier period in spring 2022 (up to 5.8 m s⁻¹, corresponding to the Eunice storm) and in autumn 2023 (up to 9.7 m s⁻¹, corresponding to the Ciarán storm) (Fig. S2).

Figure 2 Anomaly of annual air temperature (°C) as a function of anomaly of annual precipitation (%) from 2003 to 2023 with respect to the average of the 1991–2020 period (11 °C and 837 mm, respectively). Small grey hexagons represent anomaly values for the years 2003–2020 and larger blue hexagons represent anomaly values for the years 2021–2023. The air temperature anomaly was correlated to the precipitation anomaly for the years 2003–2020 (Pearson $r=-\mathrm{0.57}$, p=0.0147, n=18). Note the anomalous rainy year in 2023 relative to the pattern of precipitation as a function of temperature anomalies for the other years, possibly in response to the strong El Niño event of 2023 (Chen et al., 2024).

The four sampled ponds are situated in the periphery of the city of Brussels, with the Silex pond being bordered by the Sonian Forest (Fig. 1). The four ponds are relatively small (0.7–3.2 ha) and shallow (0.6–1.4 m) and have not been drained or dredged since at least 2018 (Table S2). The four studied ponds had significantly different Chl-a concentration values during summer, with the Leybeek pond having higher Chl-a (78.8 ± 49.5 µ g L⁻¹), followed by the Pêcheries pond (19.1 ± 13.7 µg L⁻¹), the Tenreuken pond (3.3 ± 2.4 µg L⁻¹), and the Silex pond (1.0 ± 1.2 µg L⁻¹) (Figs. 1, 3, Table S3). The Leybeek and Pêcheries ponds with higher summer Chl-a concentration had turbid water (summer TSM = 48.7 ± 36.2 and 13.7 ± 10.7 mg L⁻¹, respectively) and undetectable submerged macrophyte cover in summer (Fig. 1, Table S1). The Tenreuken and Silex ponds with lower summer Chl-a concentrations had clear water (summer TSM = 4.9 ± 3.2 and 4.0 ± 3.2 mg L⁻¹, respectively) and a high total macrophyte cover during summer (68 % and 100 %, respectively, Fig. 1, Table S1). Seasonally, the highest values of Chl-a were observed in summer in the turbid-water Leybeek and Pêcheries ponds, related to algal blooms. Conversely, lowest values of Chl-a were observed in summer in the clear-water Tenreuken and Silex ponds (Figs. 1, 3), probably related to competition for dissolved inorganic nutrients with macrophytes.

The %O₂ values ranged from 11 % to 191 % (Fig. 3). The highest %O₂ values in the four ponds were observed in spring and summer compared to autumn and winter owing to aquatic primary production. In summer, %O₂ was significantly higher in the Leybeek pond (109 ± 46 %), characterized by higher Chl-a concentration compared to the Pêcheries pond (75 ± 23 %, p=0.0212, Table S3). The lowest average %O₂ was observed in autumn in the Pêcheries pond (46 ± 22 %) and was significantly lower than in the Leybeek (85 ± 34 %, p=0.0146, Table S3) and Silex ponds (81 ± 19 %, p=0.0130, Table S3).

Figure 3 Seasonal variations of Chlorophyll-a (Chl-a, in µg L⁻¹), total suspended matter (TSM, in mg L⁻¹), oxygen saturation (%O₂, in %), partial pressure of CO₂ (pCO₂ in ppm), dissolved CH₄ concentration (CH₄, in nmol L⁻¹), and N₂O saturation level (%N₂O, in %) in four urban ponds (Leybeek (Leyb), Pêcheries (Pech), Tenreuken (Trk), and Silex (Slx)) in the city of Brussels (Belgium) from June 2021 to December 2023. Box plots show median (horizontal line), mean (cross), and 25 %–75 % percentiles (box limits). Whiskers extend from minimum to maximum values. Grey and white bands in the plots on the right correspond to the autumn/winter and spring/summer periods, respectively, and dotted vertical bars indicate the first day of each season. Lower-case letters indicate significant differences between ponds (Tables S3 and S4).

The pCO₂ values ranged from 40 to 13 804 ppm (Fig. 3). Undersaturation of CO₂ with respect to atmospheric equilibrium (∼ 410 ppm) was only observed on five occasions out of the 187 measurements, three times in the turbid-water Leybeek pond in summer (40 ppm on 13 August 2021, 220 ppm on 27 June 2022 and 149 ppm on 13 June 2023) and twice in the clear-water Tenreuken pond in spring and summer (383 ppm on 13 August 2021 and 55 ppm on 2 May 2022). Low values of pCO₂ were generally observed in spring and summer, and high values of pCO₂ were observed in autumn in the four ponds (Fig. 3). In summer, pCO₂ was lower in the Leybeek pond (2187 ± 2012 ppm) than in the Pêcheries (3427 ± 1672 ppm, p=0.0015, Table S3) and the Silex (3222 ± 1175 ppm, p=0.0002, Table S3) ponds. When data were pooled, pCO₂ was negatively influenced by %O₂ and positively by DIN, SRP, and precipitation (Table S4). In individual ponds, pCO₂ was negatively influenced by %O₂ and positively by precipitation in the four ponds, positively by DIN in the Leybeek pond, positively by DIN and SRP in the Tenreuken pond, and negatively by Chl-a in the Silex pond (Table S5).

The CH₄ dissolved concentrations ranged from 194 to 48 380 nmol L⁻¹ (Fig. 3) and were always above saturation (∼ 2 nmol L⁻¹). High values of CH₄ dissolved concentrations were generally observed in spring and summer, and low values of CH₄ dissolved concentrations were generally observed in winter in the four ponds (Fig. 3). In summer, CH₄ dissolved concentration was higher in the Silex pond (4898 ± 3384 nmol L⁻¹) than in the Pêcheries (2518 ± 2105 nmol L⁻¹, p=0.0385, Table S3) and the Tenreuken (2189 ± 1365 nmol L⁻¹, p=0.0055, Table S3) ponds. When data were pooled, dissolved CH₄ concentration was influenced positively by T_w (Table S4). In individual ponds, CH₄ dissolved concentration was also influenced positively by T_w in each of the four ponds (Table S5). Additionally, CH₄ dissolved concentration was positively influenced by precipitation in the Leybeek pond and by SRP in the Silex pond and negatively by DIN in the Pêcheries pond and by Chl-a in the Tenreuken and the Silex ponds. (Table S5). These relationships between CH₄ and other variables (SRP, DIN, Chl-a) probably indirectly reflected the seasonal variations of these other variables that were also influenced by T_w. Indeed, DIN was negatively influenced by T_w in the Pêcheries pond, Chl-a was negatively influenced by T_w in the Tenreuken and the Silex ponds, and SRP was positively influenced by T_w in the Silex pond (Table S6).

The %N₂O values ranged from 32 % to 826 % (Fig. 3). Undersaturation of N₂O with respect to atmospheric equilibrium was observed 66 times out of the 187 measurements. Low values of %N₂O were generally observed in spring and summer, and high values of %N₂O were generally observed in autumn and winter in the four ponds (Fig. 3). During spring, %N₂O was lower in the Pêcheries pond (90 ± 11 %) than the Leybeek (138 ± 30 %, p=0.0043, Table S3) and the Tenreuken (138 ± 41, p=0.0057, Table S3) ponds. During summer, %N₂O was lower in the Pêcheries pond (78 ± 17 %) than the Leybeek (191 ± 104 %, p<0.0001, Table S3) and the Silex (126 ± 49 %, p=0.001, Table S3) pond and lower in the Tenreuken pond (133 ± 106 %) than the Leybeek pond (p=0.0219, Table S3). During autumn, %N₂O was lower in the Pêcheries pond (103 ± 33 %) than the Leybeek pond (190 ± 70 %, p=0.0174, Table S3). For the whole sampling period, %N₂O was lower in the Pêcheries pond (94 ± 28 %) than the Leybeek (178 ± 82 %, p<0.0001, Table S7), the Tenreuken (140 ± 77 %, p<0.0001, Table S7), and the Silex (144 ± 113 %, p<0.0001, Table S7) ponds and was lower in the Tenreuken pond than the Leybeek pond (p=0.0038, Table S7). When data were pooled, %N₂O was influenced negatively by T_w and positively by DIN and NH${}_{\mathrm{4}}^{+}$ (Table S4). In individual ponds, %N₂O was influenced negatively by T_w in the Leybeek, the Pêcheries, and the Tenreuken ponds (Table S5). %N₂O was influenced positively by NO${}_{\mathrm{3}}^{-}$ in the Leybeek pond and by NH${}_{\mathrm{4}}^{+}$ in the Pêcheries and Tenreuken ponds (Table S8). %N₂O was influenced positively by Chl-a and TSM in the Tenreuken pond and negatively by Chl-a in the Leybeek pond (Table S5), probably reflecting the negative influence on Chl-a and TSM by T_w in the Tenreuken pond and the positive influence on Chl-a by T_w in the Leybeek pond (Table S6).

3.2 Drivers of bubble flux

The F_bubble measured with inverted funnels in the four sampled ponds in the city of Brussels ranged between 0 and 2078 mL m⁻² d⁻¹ and was influenced positively by T_w in all four systems (Fig. 4). The mean %CH₄ of the bubbles in the four sampled ponds in the city of Brussels was 31 ± 21 %, and values were influenced positively by T_w in the Silex pond (Fig. 4). The %CH₄ of the bubbles was correlated with F_bubble (Fig. S3) as both variables were positively influenced by T_w (Fig. 4).

Figure 4 Bubble flux (F_bubble in mL m⁻² d⁻¹) and the relative CH₄ content in bubbles (%CH₄, in %) as a function of surface water temperature (T_w in °C) in four urban ponds (Leybeek, Pêcheries, Tenreuken, and Silex) in the city of Brussels (Belgium) from June 2021 to December 2023. Bubbles fluxes were measured with three bubble traps in spring, summer, and autumn of 2022 and 2023, totalling 8 d in the Leybeek, Pêcheries, and Tenreuken ponds and 24 d in the Silex pond. Given the shallowness of the sampled systems (<1.5 m, Fig. 1), we assumed that sediments experience the same temperature as surface waters. In panel (a), solid lines represent linear regression of log₁₀(F_bubble) as a function of T_w for the Leybeek (log${}_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathrm{0.0664}\times T_{\mathrm{w}}+\mathrm{1.3095}$, n=22, p<0.0001), the Pêcheries (log${}_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathrm{0.0486}\times T_{\mathrm{w}}+\mathrm{1.3257}$, n=22, p=0.0146), the Tenreuken (log${}_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathrm{0.0492}\times T_{\mathrm{w}}+\mathrm{1.7039}$, n=19, p<0.0001), and the Silex (log${}_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathrm{0.0945}\times T_{\mathrm{w}}+\mathrm{1.0373}$, n=72, p<0.0001) ponds. In panel (b), the solid line represents the linear regression of %CH₄ as a function of log₁₀(T_w) for the Silex pond (%CH${}_{\mathrm{4}}=\mathrm{101.11}\times {\log }_{\mathrm{10}}\left(T_{\mathrm{w}}\right)-\mathrm{87.8}$, n=72, p<0.0001).

The time series at the Silex pond allowed the effects of T_w and atmospheric pressure variations on F_bubble to be investigated in more detail (Fig. 5). In spring 2022, the F_bubble at the Silex pond increased during drops in atmospheric pressure (depressions) (Fig. 5). There was no relation between wind speed and peaks of F_bubble (r²=0.01, p=0.463), suggesting a more important role of changes in atmospheric pressure than of wind speed in triggering bubble fluxes in the Silex pond in spring 2022. The F_bubble values at the Silex pond were higher in summer (1152 ± 433 mL m⁻² d⁻¹) than during spring (198 ± 170 mL m⁻² d⁻¹), and the temporal changes of F_bubble tracked those of T_w (Fig. 5). In order to evaluate the relative importance of changes of atmospheric pressure and water temperature in triggering bubble fluxes, the F_bubble was modelled as a function of T_w alone or as a function of both T_w and Δp (Figs. 5, S4). For periods of low T_w (<15 °C), the inclusion of the term of Δp in the model improved the performance of the model by comparison to the measurements (Figs. 5, S4). But for warmer periods (>15 °C), when bubbling fluxes were quantitatively more important, the inclusion in the model of the term of Δp did not improve the performance of the model (Figs. 5, S4). For the full T_w range (<15 and >15 °C), the inclusion of the term of Δp only improved the performance of the model very marginally (Fig. S4).

Figure 5 Time series of surface water temperature (T_w, °C), wind speed (m s⁻¹), atmospheric pressure (atm), and measured and modelled bubble flux (F_bubble in mL m⁻² d⁻¹) in the Silex pond from 29 March 2022 to 15 April 2022 and from 18 July 2022 to 23 July 2022. The bubble flux was modelled from a fit to data based on T_w alone (${\log }_{\mathrm{10}}\left(F_{\mathrm{bubble}}\right)=\mathrm{3973}\times {\log }_{\mathrm{10}}\left(T_{\mathrm{w}}\right)-\mathrm{215}$, p<0.0001, n=72) and based on both T_w and drops in atmospheric pressure (Δp) $\left({\log }_{\mathrm{10}}\right(F_{\mathrm{bubble}})=\mathrm{4.551}\times {\log }_{\mathrm{10}}(T_{\mathrm{w}})+\mathrm{1.962}\times \mathrm{\Delta }p-\mathrm{3.006}$, p<0.0001, n=72).

3.3 Drivers of methane ebullitive fluxes

The $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the four ponds ranged between 0 and 59 mmol m⁻² d⁻¹ and were positively related to T_w (Fig. 6). The fitted relations between $E_{{\mathrm{CH}}_{\mathrm{4}}}$ and T_w were specific to each pond and encompassed the fitted relations established in similar systems: four small ponds in Quebec (DelSontro et al., 2016) and a small urban pond in the Netherlands (Aben et al., 2017). The Q₁₀ of CH₄ ebullition values ranged between 4.4 in the deeper Pêcheries pond and 26.9 in the shallower Leybeek pond (Table S9). The Q₁₀ values of CH₄ ebullition in the four studied ponds of the city of Brussels, in Quebec (DelSontro et al., 2016), and in the Netherlands (Aben et al., 2017) were negatively related to water depth (Fig. 6).

Figure 6 (a) Measured ebullitive CH₄ fluxes (mmol m⁻² d⁻¹) as a function of surface water temperature (°C) in four urban ponds (Leybeek, Pêcheries, Tenreuken, and Silex) in the city of Brussels (Belgium), in spring, summer, and autumn of 2022 and 2023, totalling 8 d in the Leybeek, Pêcheries, and Tenreuken ponds and 24 d in the Silex pond, with three bubble traps. Solid lines represent exponential fit for the Leybeek ($Y=\mathrm{0.01}\cdot e^{\mathrm{0.32}\cdot X}$, n=22), Pêcheries ($Y=\mathrm{0.16}\cdot e^{\mathrm{0.15}\cdot X}$, n=22), Tenreuken ($Y=\mathrm{0.10}\cdot e^{\mathrm{0.23}\cdot X}$, n=19), and Silex ($Y=\mathrm{0.54}\cdot e^{\mathrm{0.18}\cdot X}$, n=72) ponds (Table S7). Dashed lines represent the published exponential fit established in similar systems: four small ponds in Quebec ($Y=\mathrm{0.06}\cdot e^{\mathrm{0.25}\cdot X}$) (DelSontro et al., 2016) and a small urban pond in the Netherlands ($Y=\mathrm{0.51}\cdot e^{\mathrm{0.17}\cdot X}$) (Aben et al., 2017). (b) Each exponential curve allows us to determine a Q₁₀ of CH₄ ebullition, plotted against water depth; the solid line represents linear regression (Y=30.64–19.67⋅X, n=6).

3.4 Relative contribution of methane ebullitive and diffusive fluxes

Diffusive CH₄ fluxes computed from dissolved CH₄ concentration and k derived from wind speed ranged between 0.1 and 19.7 mmol m⁻² d⁻¹ (Fig. 7). The diffusive CH₄ fluxes tended to be higher in summer and spring than in autumn and winter, owing to the strong positive influence on CH₄ dissolved concentration by T_w (Fig. 3; Tables S4, S5). In addition, wind speed only showed small seasonal variations during sampling (0.6 ± 0.6 m s⁻¹ in spring, 0.3 ± 0.2 m s⁻¹ in summer, 0.7 ± 0.7 m s⁻¹ in autumn, and 0.6 ± 0.2 m s⁻¹ in winter) (Fig. S2). $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values were calculated from the relations with T_w for each pond given in Figure 6 from the T_w data coincident with the diffusive CH₄ fluxes (Fig. 7). The resulting calculated $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values allowed us both to compare and integrate components of CH₄ emissions to the atmosphere seasonally and to calculate the relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to total (diffusive + ebullitive) CH₄ emissions. The relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to total CH₄ emissions ranged between 1 % and 99 % in the four sampled ponds in the city of Brussels (Fig. 7) and was influenced positively by T_w (Fig. S5). The values of Q₁₀ of diffusive CH₄ fluxes were lower than those for $E_{{\mathrm{CH}}_{\mathrm{4}}}$ in each pond and less variable seasonally (1.2 in the Pêcheries pond to 2.9 in the Silex pond) (Table S9).

Figure 7 Seasonal variations of (a, b) diffusive and (c, d) ebullitive CH₄ fluxes (mmol m⁻² d⁻¹) and (e, f) the ratio of ebullitive CH₄ flux to total (ebullitive + diffusive) CH₄ flux (%) in four urban ponds (Leybeek (Leyb), Pêcheries (Pech), Tenreuken (Trk), and Silex (Slx)) in the city of Brussels (Belgium) from June 2021 to December 2023. Diffusive fluxes were calculated from CH₄ concentration and gas transfer velocity derived from wind speed. Ebullitive CH₄ fluxes were calculated from the relations with water temperature for each pond (Fig. 6; Table S7) from the water temperature data coincident with the diffusive CH₄ fluxes. Note that the relations of ebullitive CH₄ fluxes as a function of water temperature were established over a temperature range (7.0 to 26.3 °C) that is consistent with the range of water temperature values (2.0–25.9 °C) over which the ebullitive CH₄ fluxes were modelled. Box plots show median (horizontal line), mean (cross), and 25 %–75 % percentiles (box limits). Whiskers extend from minimum to maximum values. Grey and white bands in the plots on the right correspond to the autumn/winter and spring/summer periods, respectively, and dotted vertical bars indicate the first day of each season. Lower-case letters indicate significant differences between ponds (Tables S3 and S4).

The annually averaged diffusive and ebullitive fluxes of CH₄ in the four ponds in the city of Brussels were plotted against annually averaged Chl-a concentration, total macrophyte cover in summer, water depth, and lake surface area (Fig. 8) that are frequent predictors of variations of CH₄ fluxes among lakes and ponds (Holgerson and Raymond, 2016; DelSontro et al., 2018, Deemer and Holgerson, 2021; Casas-Ruiz et al., 2021; Borges et al., 2022). The annually averaged $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values were significantly higher in the two clear-water ponds (7.3 ± 2.9 and 13.4 ± 3.7 mmol m⁻² d⁻¹ in the Tenreuken and the Silex ponds, respectively) than the two turbid-water ponds (3.8 ± 3.2 and 2.5 ± 1.4 mmol m⁻² d⁻¹ in the Leybeek and the Pêcheries ponds, respectively) (Table S7). The annually averaged $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values were significantly higher in the Silex pond, which showed a higher macrophyte cover during summer (100 % in the Silex pond and 68 % in the Tenreuken pond) than the Tenreuken pond (p<0.0001, Table S7), and were not significantly different in the two turbid-water Leybeek and Pêcheries ponds (p=0.0617, Table S7) that showed similar low macrophyte cover during summer (6 % and 9 % in the Leybeek and Pêcheries ponds, respectively) (Fig. 8). The annually averaged $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values were overall influenced positively by macrophyte cover and negatively by Chl-a (Fig. 8).

Figure 8 Mean diffusive and ebullitive CH₄ fluxes (mmol m⁻² d⁻¹) and mean ratio of ebullitive CH₄ flux to total (diffusive + ebullitive) CH₄ flux (%) versus chlorophyll-a (Chl-a, in µg L⁻¹), total macrophyte cover in summer (%), water depth (m), and lake surface area (ha) in four ponds (Leybeek, Pêcheries, Tenreuken, and Silex) in the city of Brussels (Belgium) from June 2021 to December 2023. Error bars indicate the standard deviation. Solid lines indicate either linear or polynomial fits. Statistical comparisons between the four ponds are summarized in Table S3.

In the four sampled urban ponds, annually averaged CH₄ diffusive fluxes were higher in the pond with the highest total macrophyte cover in the clear-water ponds and higher in the pond with the highest Chl-a concentration in the turbid-water ponds (Fig. 8). The annually averaged relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to total CH₄ emissions was higher in the two clear-water ponds than the two turbid-water ponds (Table S7). The relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to the total CH₄ flux was influenced positively by macrophyte cover and negatively by Chl-a (Fig. 8).

The annually averaged diffusive fluxes of CO₂ ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) and N₂O ($F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$) in the four ponds in the city of Brussels were also plotted against annually averaged Chl-a concentration, total macrophyte cover in summer, water depth, and lake surface area (Fig. S6). Annually averaged $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values were lower in the Leybeek pond than the Pêcheries and the Silex ponds (Table S7). $F_{{\mathrm{CO}}_{\mathrm{2}}}$ was not significantly influenced by the other variables (Chl-a concentration, total macrophyte cover, water depth, and lake surface area) (Fig. S6). Annually averaged $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ was not significantly different between clear-water and turbid-water ponds. $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ was significantly lower in the deeper Pêcheries pond than the two shallower Leybeek and Silex ponds (Table S7), and $F_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ showed a significant negative relationship with water depth (Fig. S6).

3.5 Relative contribution of CO₂, CH₄, and N₂O emissions

The emissions in CO₂-eq for the three GHGs averaged per season for both 2022 and 2023 peaked seasonally in summer in the Silex (2.9 mg CO₂-eq m⁻² d⁻¹), the Tenreuken (1.7 mg CO₂-eq m⁻² d⁻¹), and the Leybeek (1.1 mg CO₂-eq m⁻² d⁻¹) pond (Fig. 9) but peaked in autumn in the Pêcheries pond (1.3 mg CO₂-eq m⁻² d⁻¹). The higher value of the total GHG emissions in autumn compared to other seasons in the Pêcheries pond was due to an increase of CO₂ emissions in autumn that surpassed the peak of CH₄ emissions in summer. The GHG fluxes were lowest in winter in the Silex (1.3 mg CO₂-eq m⁻² d⁻¹), the Tenreuken (0.9 mg CO₂-eq m⁻² d⁻¹), the Pêcheries (0.8 mg CO₂-eq m⁻² d⁻¹), and the Leybeek (0.6 mg CO₂-eq m⁻² d⁻¹) ponds. The relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ peaked in summer in the Silex (73.8 %), the Tenreuken (70.9 %), the Pêcheries (23.6 %), and the Leybeek (58.3 %) ponds. The relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ was lowest in winter in the Silex (22.1 %), the Tenreuken (10.0 %), the Pêcheries (6.7 %), and the Leybeek (1.0 %) ponds.

Figure 9 Seasonal and year-to-year variations of the emissions to the atmosphere of CO₂ (diffusive), CH₄ (diffusive and ebullitive), and N₂O (diffusive) expressed in CO₂ equivalents (in mg CO₂-eq m⁻² d⁻¹) and their relative contribution (in %), in four urban ponds (Leybeek, Pêcheries, Tenreuken, and Silex) in the city of Brussels (Belgium) from June 2021 to December 2023. Seasonal averages include data from 2021, 2022, and 2023. The annual precipitation was higher in 2023 (1011 mm) than in 2022 (701 mm).

The annual emissions in CO₂-eq of the three GHGs (CO₂, CH₄, and N₂O) in 2022 and 2023 were higher in the two clear-water ponds (1.3 ± 0.5 and 1.8 ± 0.9 mg CO₂-eq m⁻² d⁻¹ in the Tenreuken and Silex ponds, respectively) than in the two turbid-water ponds (1.0 ± 0.2 and 0.9 ± 0.5 mg CO₂-eq m⁻² d⁻¹ in the Leybeek and Pêcheries ponds, respectively) (Fig. 9) due to higher total CH₄ emissions (diffusive + ebullitive) in clear-water ponds (0.7 ± 0.4 and 1.2 ± 0.5 mg CO₂-eq m⁻² d⁻¹ in the Tenreuken and Silex ponds, respectively) than in turbid-water ponds (0.2 ± 0.2 and 0.4 ± 0.3 mg CO₂-eq m⁻² d⁻¹ in the Leybeek and Pêcheries ponds, respectively). The contribution of N₂O to the total GHG emissions was marginal and did not affect the differences in total GHG fluxes between ponds, with the highest contribution observed in the Leybeek pond, with a contribution of 1.7 %.

The majority of GHG emissions in CO₂-eq were related to CO₂ and CH₄ (diffusive + ebullitive) in the four ponds. In turbid-water ponds CO₂ represented the largest fraction of GHG emissions (68.5 % (2022) and 79.3 % (2023) in the Pêcheries pond, and 49.0 % (2022) and 58.3 % (2023) in the Leybeek pond). In clear-water ponds CH₄ represented the largest fraction of GHG emissions (66.5 % (2022) and 63.3 % (2023) in the Silex pond and 60.8 % (2022) and 50.0 % (2023) in the Tenreuken pond). The higher annual GHG emissions in CO₂-eq from the two clear-water ponds than from the turbid-water ponds were related to the higher contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$.

The annual GHG fluxes increased from 2022 to 2023 due to an increase in the relative contribution of CO₂ diffusive emissions in all four ponds. Diffusive CO₂ emissions averaged 0.5 mg CO₂ m⁻² d⁻¹ annually in 2022 and 0.7 mg CO₂ m⁻² d⁻¹ annually in 2023 in all four ponds. Diffusive CO₂ emissions were 2.1 times higher in summer 2023 than in summer 2022 and 2.5 times higher in autumn 2023 than in autumn 2022 and showed similar values between 2023 and 2022 in spring and winter (1.1 higher and 1.1 lower, respectively).

4 Discussion

The Leybeek and Pêcheries ponds are turbid-water systems (high Chl-a and TSM values, low submerged macrophyte cover), and the Tenreuken and Silex ponds are clear-water systems (low Chl-a and TSM values, high submerged macrophyte cover) (Figs. 1, 3). All four ponds have a relatively similar size (0.7 to 3.2 ha) and depth (0.5 to 1.4 m) and are uniformly located in an urban landscape in the city of Brussels. It can be assumed that, among the four systems, the major difference that is expected to affect GHG emissions is the dominance of aquatic primary producer, either phytoplankton or macrophytes, corresponding to two alternative states sensu Scheffer et al. (1993). Our data set provides the opportunity to investigate the effect of the two alternative states on GHG emissions from small lentic systems.

The reported pCO₂ values (40 to 13 804 ppm) (Fig. 3) in the four ponds in the city of Brussels were within the range of values typically observed in ponds (Holgerson and Raymond, 2016; Peacock et al., 2019; Audet et al., 2020) (Fig. 3). The pCO₂ values were influenced negatively by %O₂ and positively by DIN and SRP across seasons (Tables S4, S5), showing that their seasonal variability was driven by aquatic primary production and degradation of organic matter (e.g. Holgerson, 2015). Accordingly, low values of pCO₂ were generally observed in spring and summer, probably due to uptake of CO₂ by primary production from either phytoplankton or submerged macrophytes. High values of pCO₂ were observed in autumn in the four ponds and probably reflect the release of CO₂ from degradation of organic matter due to the senescence of phytoplankton or macrophytes (Fig. 3). In all four ponds, pCO₂ values were influenced positively by precipitation (Tables S4, S5), suggesting an additional control of external inputs of carbon either as organic carbon sustaining internal degradation of organic matter or as soil CO₂ (e.g. Marotta et al., 2010; Ojala et al., 2011; Rasilo et al., 2012; Vachon and del Giorgio, 2014; Holgerson, 2015). The %N₂O values (32 % to 826 %) (Fig. 3) in the four ponds were within the range of values typically observed in ponds (Audet et al., 2020; Rabaey and Cotner, 2022). When all the data were pooled, the %N₂O was influenced positively by DIN (Table S4), as also frequently reported by other studies in ponds and interpreted as a control of nitrification and/or denitrification (hence N₂O production) by DIN levels (Audet et al., 2020; Webb et al., 2021; Wang et al., 2021; Xie et al., 2024). The negative influence on %N₂O by T_w (Table S4) might reflect the effect of the inhibition at low temperatures of the final step of denitrification leading to an accumulation of N₂O (Velthuis and Veraart, 2022) but could also indirectly result from the higher DIN concentrations observed at low T_w values (Table S6). The CH₄ dissolved concentrations (194 to 48 380 nmol L⁻¹) (Fig. 3) in the four ponds were within the range of values typically observed in ponds (Natchimuthu et al., 2014; Holgerson and Raymond, 2016; Peacock et al., 2019; Audet et al., 2020; Rabaey and Cotner, 2022; Ray et al., 2023) and were influenced positively by T_w in all four ponds individually and, when pooled (Tables S4, S5), most probably reflecting the increase in sedimentary methanogenesis with temperature (Schulz and Conrad, 1996).

Temperature also exerted a strong control on bubble flux from sediments and ebullitive CH₄ emissions. The F_bubble values (0 and 2078 mL m⁻² d⁻¹) in the four sampled ponds (Fig. 4) were within the range of values reported in lentic systems of equivalent size by Wik et al. (2013) (0 to 2772 mL m⁻² d⁻¹), DelSontro et al. (2016) (11 to 748 mL m⁻² d⁻¹), and Ray and Holgerson (2023) (0 to 2079 mL m⁻² d⁻¹). The F_bubble was influenced positively by T_w (Fig. 4), in agreement with previous studies (e.g. Wik et al., 2013; DelSontro et al., 2016; Aben et al., 2017; Ray and Holgerson, 2023). Bubbling events from lake sediments are known to also be triggered by a decrease in hydrostatic pressure on the sediments due to water level fluctuations or drops in atmospheric pressure (Tokida et al., 2007; Scandella et al., 2011; Varadharajan and Hemond, 2012; Wik et al., 2013; Taoka et al., 2020; Zhao et al., 2021). In the Silex pond, in spring 2022, some peaks in F_bubble were related to drops in atmospheric pressure (Fig. 5) but unrelated to wind speed as shown in Gatun Lake (Keller and Stallard, 1994). A statistical model of F_bubble that included the contributions of T_w and Δp was used to quantify the relative importance of each of these two drivers (Fig. S4) and showed that an air pressure drop only seemed quantitatively important at low T_w and that the intensity of bubble flux was mainly driven by temperature change at yearly scales, in agreement with previous studies (e.g. Wik et al., 2013; DelSontro et al., 2016; Aben et al., 2017; Ray and Holgerson, 2023).

The mean %CH₄ of the bubbles (31 ± 21 %) in the four sampled ponds in the city of Brussels was comparable to the values obtained by Wik et al. (2013) (35 ± 25 %), DelSontro et al. (2016) (58 ± 25 %), and Ray and Holgerson (2023) (25 ± 13 %) in lentic systems of similar size. The increasing pattern of %CH₄ of the bubbles with T_w (Fig. 4) was most probably related to the strong dependence of methanogenesis on temperature (Schulz and Conrad, 1996). The increase in methanogenesis with temperature leads to the build-up of gas bubbles in sediments that are richer in CH₄, and consequently to higher bubble fluxes with a higher CH₄ content at higher temperatures (Figs. 4, S3). Since both F_bubble and the %CH₄ of the bubbles increased with T_w (Fig. 4), $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the four ponds were also positively related to T_w (Fig. 6), as shown previously in other small lentic systems (e.g. Wik et al., 2013; DelSontro et al., 2016; Natchimuthu et al., 2016; Aben et al., 2017; Ray and Holgerson, 2023; Rabaey and Cotner, 2024). Yet, the dependency of CH₄ ebullition on temperature (Q₁₀) was different among the four ponds and was negatively related to depth including data from systems in Quebec (DelSontro et al., 2016) and the Netherlands (Aben et al., 2017) (Fig. 6). This implies that an increase in T_w leads to a smaller increase in CH₄ ebullitive fluxes (lower Q₁₀) in deeper ponds as the impact of hydrostatic pressure on sediments is higher in deeper ponds compared to shallow ponds, restricting bubble formation and release (e.g. DelSontro et al., 2016). This dependence of Q₁₀ of CH₄ ebullition on depth suggests that the response of CH₄ ebullition to heat waves (or longer-term warming) might be more intense the shallower the pond, in addition to other effects from heat waves on GHG emissions (e.g. Audet et al., 2017).

The values of Q₁₀ for diffusive CH₄ fluxes in the four ponds were lower than those for ebullitive CH₄ fluxes (Table S9), as reported by other studies in lentic systems (DelSontro et al., 2016; Xun et al., 2024). The lower dependence on T_w of diffusive CH₄ fluxes compared to ebullitive CH₄ fluxes might be related to a lower relative change of CH₄ concentrations and k with the variation of T_w. Dissolved CH₄ concentrations in surface waters of lentic systems are strongly affected by microbial methane oxidation (e.g. Bastviken et al., 2002). A relative increase in methanogenesis in sediments might lead to a stronger increase in CH₄ emission by ebullition than by diffusion because of a reduction of CH₄ diffusive emissions resulting from microbial methane oxidation. Additionally, k depends on wind speed, but in the four ponds, the warmer periods of the year (summer) tended to be less windy (∼ 0.3 m s⁻¹) than the other seasons (>0.6 m s⁻¹), also contributing to a lower dependence on T_w of CH₄ diffusive fluxes compared to ebullitive fluxes (lower Q₁₀ values).

The difference in the Q₁₀ of diffusive and ebullitive CH₄ fluxes was consistent with a seasonally variable contribution of the diffusive and ebullitive CH₄ fluxes as a function of T_w. The contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to total (diffusive + ebullitive) CH₄ emissions strongly increased with T_w in the four ponds (Fig. S5). At an annual scale, $E_{{\mathrm{CH}}_{\mathrm{4}}}$ represented between 55 % and 83 % of the total CH₄ emissions in the Leybeek and Silex ponds, respectively. This finding is consistent with other studies showing that ebullitive CH₄ fluxes can account for more than half of total CH₄ emissions in small and shallow lentic systems (e.g. Wik et al., 2013; Deemer and Holgerson, 2021; Ray and Holgerson, 2023; Rabaey and Cotner, 2024). The averaged $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values were higher in the two clear-water ponds (10.4 mmol m⁻² d⁻¹) than the two turbid-water ponds (3.2 mmol m⁻² d⁻¹) (Fig. 7). The averaged $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the four ponds were influenced positively by macrophyte cover and negatively by Chl-a (Fig. 8). The higher $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values from the two clear-water ponds would suggest that the delivery of organic matter to sediments from macrophytes sustained a quantitatively larger methane production than from phytoplankton. This finding is consistent with the notion that vegetated littoral zones of lakes are hotspots of CH₄ production and emission (e.g. Hyvönen et al., 1998; Huttunen et al., 2003; Juutinen et al., 2003; Desrosiers et al., 2022). CH₄ fluxes in lentic systems have been extrapolated at a global scale assuming a dependency on aquatic productivity using Chl-a as a predictor (e.g. DelSontro et al., 2018). The negative relation between $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values with Chl-a (Fig. 8) shows that Chl-a concentration alone fails to predict ebullitive fluxes in macrophyte-dominated clear-water ponds.

The annually averaged diffusive CH₄ emissions in the four ponds seemed to respond positively to both increasing phytoplankton and macrophyte biomass, resulting in a U-shaped relation between diffusive CH₄ emissions and Chl-a, as well as macrophyte cover (Fig. 8). Higher values of annually averaged CH₄ diffusive fluxes occurred at the extreme values of Chl-a or of macrophyte cover (minimum or maximum), and lower values occurred at the intermediate values of Chl-a or macrophyte cover. Such a U-shape relation resulted from the inverse relationship between macrophyte cover and Chl-a (alternative states) and is consistent with the reported positive relation between diffusive CH₄ fluxes with both macrophyte cover (e.g. Ray et al., 2023; Theus et al., 2023) and phytoplankton biomass (e.g. DelSontro et al., 2018; Yan et al., 2019; Bartosiewicz et al., 2021). The relative contribution of $E_{{\mathrm{CH}}_{\mathrm{4}}}$ to the total annual CH₄ flux increased with the macrophyte cover (Fig. 8), in agreement with the idea of an increase in CH₄ ebullition relative to diffusive CH₄ emissions in vegetated sediments compared to unvegetated sediments (e.g. Desrosiers et al., 2022; Ray et al., 2023; Theus et al., 2023).

Fluxes of CH₄ and CO₂ have been reported to be negatively related to surface area and depth by numerous studies in ponds (e.g. Holgerson, 2015; Holgerson and Raymond, 2016; Ray et al., 2023; Theus et al., 2023) and lakes (e.g. Kankaala et al., 2013; DelSontro et al., 2018, Deemer and Holgerson, 2021; Casas-Ruiz et al., 2021; Borges et al., 2022). Annual diffusive $F_{{\mathrm{CH}}_{\mathrm{4}}}$ and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ were both unrelated to surface area and depth in the four studied ponds (Figs. 8, S6) resulting from the narrow range of variation of water depth (0.6 to 1.4 m) and surface area (0.7 to 3.2 ha). The lack of relationship between annual $F_{{\mathrm{CO}}_{\mathrm{2}}}$ and both Chl-a and macrophyte cover in the four ponds (Fig. S6) might be surprising since other studies have reported lower CO₂ fluxes in more productive lentic systems (e.g. Sand-Jensen and Staehr, 2007; Borges et al., 2022). We hypothesize that given that the four systems were either phytoplankton-dominated or macrophyte-dominated (alternative states), the ponds had an important submerged productivity, in both cases, resulting in a relatively invariant $F_{{\mathrm{CO}}_{\mathrm{2}}}$ as a function of either Chl-a or macrophyte cover.

Global average emissions of GHGs in CO₂-eq from inland waters are dominated by CO₂ followed by CH₄ with a small contribution from N₂O according to Lauerwald et al. (2023). However, in small lentic systems such as ponds, the CO₂-eq emissions from CH₄ can match or dominate those of CO₂ (e.g. Webb et al., 2023; Ray and Holgerson, 2023; Rabaey and Cotner, 2024). The meta-analysis of Holgerson and Raymond (2016) suggested that the CO₂ and CH₄ emissions in CO₂-eq are numerically close in small lentic systems such as ponds but become increasingly dominated by CO₂ emissions in larger lentic systems. In the four studied ponds, the GHG emissions in CO₂-eq were dominated by CO₂ and CH₄ with a marginal contribution (<2 %) from N₂O (Fig. 9). Annually, CO₂ represented the largest fraction of GHG emissions in CO₂-eq (∼ 60 %) in turbid-water ponds (Leybeek and Pêcheries), while CH₄ represented the largest fraction of GHG emissions in CO₂-eq (∼ 60 %) in clear-water ponds (Silex and Tenreuken) as a result of higher $E_{{\mathrm{CH}}_{\mathrm{4}}}$ values in the clear-water ponds (Fig. 7).

The annual GHG emissions in CO₂-eq increased from 2022 to 2023 due to an increase in the relative contribution of CO₂ diffusive emissions in all four ponds (Fig. 9) as a result of higher precipitation in 2023 (Fig. 2). Air temperatures were similar in both years (annual average of 12.2 °C in 2022 and 12.1 °C in 2023), and precipitation was 1.5 times higher in 2023 than in 2022. Higher precipitation is likely to increase the inputs of organic and inorganic carbon from soils to ponds, as previously shown in other lentic systems (e.g. Marotta et al., 2010; Ojala et al., 2011; Rasilo et al., 2012; Vachon and del Giorgio, 2014; Holgerson, 2015). This hypothesis is only based on the comparison of 2 years, but the increase in the relative contribution of CO₂ diffusive emissions in 2023 was observed in all four ponds. The synchronicity of the increase in CO₂ diffusive emissions in 2023 compared to 2022 suggests a common uniform driver in all four ponds that would be consistent with a large variation in weather between the 2 years, such as annual precipitation. The strong El Niño event in 2023 induced low-level cyclonic wind anomalies and higher precipitation over western Europe, including Belgium (Chen et al., 2024).

5 Conclusions

Ebullitive CH₄ emissions in 2022–2023 were higher in the two clear-water, macrophyte-dominated ponds (Tenreuken and Silex) than in the two turbid-water, phytoplankton-dominated ponds (Pêcheries and Leybeek) of the city of Brussels, although the diffusive CH₄ fluxes were not significantly different between the clear-water ponds and the turbid-water ponds. The annually averaged diffusive N₂O and CO₂ fluxes were not significantly different in the two clear-water ponds from those in the two turbid-water ponds. Other studies have found no difference in N₂O sedimentary production in lakes with a high and low density of submerged macrophytes. We hypothesize that CO₂ fluxes were relatively invariant among the four sampled ponds because of their similar size, depth, and putatively productivity (either from phytoplankton or submerged macrophytes). The total (diffusive and ebullitive) CH₄ emissions represented 58 % of total annual GHG emissions in CO₂-eq in the two clear-water ponds compared to 41 % in the two turbid-water ponds. CO₂ represented nearly all the remainder of total annual GHG emissions in CO₂-eq, and N₂O represented a very marginal fraction (<2 %).

The seasonal variations in GHG emissions were mainly driven by CH₄ ebullitive emissions that peaked in summer (both quantitatively and relatively), as CH₄ ebullition was related positively to T_w, resulting from an increase in both the flux of bubbles and CH₄ content of bubbles with warming. The pCO₂ values in the four sampled ponds increased with precipitation at a seasonal scale, probably in relation to higher inputs of organic and inorganic carbon from soils. The years 2022 and 2023 were abnormally dry and wet, respectively. The GHG emissions were higher in 2023 than 2022, mainly due to an increase in the relative contribution of CO₂ emissions, probably due to higher precipitation in response to a strong El Niño event of 2023. This would suggest that variations in precipitation also affected year-to-year variations in CO₂ emissions in addition to partly regulating seasonal variations in CO₂ emissions from the four studied ponds.