Passaúna Reservoir is located at the Passaúna River in the southern part of Brazil, near to the city of Curitiba (25.53°S and 49.39°W, Figure 1). The reservoir was constructed in 1989 for drinking water supply. Passaúna is a polymictic and mesotrophic reservoir (Xavier et al., 2017), with an average depth of 8.3 m. It is elongated in the North-South direction with approximately 10 km length and 0.6 km width. Its main inflow is the Passaúna river with an average annual discharge of 2 m³ s⁻¹ (Carneiro et al., 2016).

According to the updated Köppen climate classification, the region is characterized by temperate oceanic climate (Cfb), with monthly mean temperatures below 22°C for all months and no significant precipitation difference between seasons (Beck et al., 2018). The mean annual precipitation in the region ranges between 1,400 and 1,600 mm. During measurements covering a full annual cycle conducted in 2018–2019, the reservoir water temperature ranged from 16 to 28°C, with August being the coldest month with an averaged bottom temperature of 16.5°C (Ishikawa et al., 2021). Meteorological data for this study was provided by the company SANEPAR which manages the reservoir. Weather data is logged at a station installed near to the dam.

Ebullition flux and pore water methane concentration were measured at three sampling sites distributed along the reservoir (Figure 1B). The site P1 is located most closely to the main river inflow (Passaúna River) at a water depth of 8 m. The site P2 was placed in the central region of the reservoir, near to the water intake facility, where the water depth is 12.5 m. Site P3, was placed in the deepest region of the reservoir, approximately 810 m upstream of the dam, at a water depth of 14 m.

The ebullition measurements analyzed in the present study are a continuation of the measurements that started in 2017 and were described by Marcon et al. (2019). The ebullition flux was monitored continuously by automated bubble traps (ABT, Senect GmbH, Germany) from January to March 2019. The ABT were attached to a surface buoy and deployed at a depth of approximately 4–5 m above the sediment surface. The device consists of a 1 m diameter funnel, which collects rising gas bubbles, and a differential pressure sensor to monitor the gas volume that accumulated during fixed time intervals (30 s). The device also measured water temperature and pressure, which were used to convert the collected gas volume to standard pressure and temperature (1,013.25 mbar and 20°C). In addition, bubbles were collected from the sediment near to the ABT’s location to estimate the methane fraction within the bubbles. The captured gas was transferred to vials with saturated salt solution until its analysis in the laboratory with an Ultraportable Greenhouse Gas Analyzer (UGGA, Los Gatos Research Inc.).

Dissolved methane concentration in the sediment pore water was measured with a dialysis pore water sampler (DPS) as described by Hilgert (2014) and Hölzlwimmer (2013). The DPS is a perforated aluminum frame of 70 cm height divided into 15 chambers of 4 cm width each. A cellulose membrane bag of 50 ml filled with ultra-pure water was added to each chamber. On 02 February 2019 three DPS’s were deployed by divers in the bottom sediment at each ABT location (P1, P2, and P3) and positioned vertically such that seven chambers were inside the sediment and eight chambers in the overlaying water column. The devices were deployed for 5 days, to allow the ultra-pure water of the bags to equilibrate with the ambient water or pore water. Directly after recovery, a sample of 5 ml was extracted from each bag with a plastic syringe, a headspace of 5 ml was created in the syringe and after rigorous shaking, the headspace gas was transferred to vials containing saturated salt solution, previously sealed with a rubber stopper and crimp-capped. The vials were stored upside-down until analysis in the laboratory, where the methane concentration was measured with an Ultraportable Greenhouse Gas Analyzer (UGGA, Los Gatos Research Inc.) in a closed loop arrangement, as described by Wilkinson et al. (2018); Wilkinson et al. (2019b). The corresponding methane concentration in the water sample ( C C H 4 in mol L⁻¹) was calculated as proposed by Bossard et al. (1981) as: C C H 4 = ( V H s V W ( X e q − X b a c k ) R T + X e q K H ) × P a t m   , where V H s is the headspace volume in L, V W is the water sample volume in L, X e q and X b a c k are the methane mole fractions measured with the gas analyzer after equilibration and in the initial headspace in ppm respectively. K H is the temperature dependent Henry gas solubility coefficient, which was calculated according to Goldenfum (2010) in mol L⁻¹ atm⁻¹, T is temperature in K, R is the gas constant ( R = 0.08205 L atm mol⁻¹ K⁻¹), and P a t m is the atmospheric pressure in atm.

The following sections describe the re-analysis and additional processing of data from acoustic surveys conducted in Passaúna Reservoir during former studies.

We developed statistical and data driven models to predict the spatial distribution of the estimated sediment gas content from maps of water depth, sediment thickness, and organic matter content (LOI 550°C) as predictor variables. The input variables and the prediction of gas content were analyzed within a spatial grid created for the reservoir based on the acoustic survey conducted in 2019. The analysis grid was created using the Geographical Information System (GIS) software ArcGIS (v. 10.2.2) with the geoprocessing tool Fishnet, in which we divided the surface area of Passaúna reservoir in 8,454 rectangular grid cells with dimensions of 33 m by 33 m (Figure 1E). Additional information about the grid selection is provided in the supplementary material (Supplementary Figure S2). The size of the grid cells was chosen to be sufficiently small to resolve the spatial heterogeneities of the measured parameters occurring in the sediment (longitudinal and transversal variations), and large enough to analyze all parameters in terms of spatial averages, that integrate small-scale (unresolved) structures. Mean values of all parameters (water depth, bottom slope, sediment thickness, organic content through LOI 550°C, and gas content) were calculated for all grid cells for which survey data are available.

To reduce the uncertainty of parameters derived from acoustic measurements, we excluded grid cells with an average bottom slope larger than 10° from all subsequent analysis. This is justified by the fact that the acoustic backscatter strength of sediment surfaces and sediment layers depends on the grazing angle of the soundwave (angle between incident wave and the tangent to the surface). For incidence angles of the soundwave near to the normal direction in relation to the sediment surface, scattering (attack values) dominate the backscattering, in contrast for inclined conditions (i.e., grazing angles in the range 30–60°) the volume scattering dominates (Fonseca et al., 2002). Sternlicht and de Moustier (2003) showed that the effect of slopes is lengthening the echo, resulting in prolongated rising and decaying parts of the echograms, leading to a reduction of the maximum backscatter independent from the sediment composition.

Steep slopes in Passaúna Reservoir occurred mainly near the banks or along the old Thalweg of the Passaúna River as described by Sotiri et al. (2021). The slope threshold resulted in removal of 1,034 grid cells, thus 4,651 grid cells were classified as valid cells (with data of all parameters measured and in agreement with the slope criterion).

The potential CH₄ production (PMP), dissolved CH₄ concentration in pore water, and ebullition were not extrapolated to the entire reservoir, as they were measured at only three locations. Nevertheless, the measurements are used to support our discussion on the spatial distribution and dynamics of the estimated sediment gas content.

Three multiple linear regression models were tested for gas content prediction in which the water depth, sediment thickness, and LOI 550°C were the predictors. The first model (MLR I) was a simple multiple linear regression, for the second model (MLR II) the predicted value (gas content) was log transformed, and for the third statistical model a stepwise multiple regression (SMR) with untransformed values was performed. The main difference of the stepwise multiple regression to the two other models is that predictors are included sequentially in the model and accepted if a p-value criterion is met for a significance level of 5%. Furthermore, interactions of the input variables are tested in this model as predictors.

The data driven model is a supervised artificial neural network (ANN). The ANN architecture had three layers, one input layer with three neurons, one hidden layer with 10 neurons, and one output layer with one neuron (Figure 2). The input variables were water depth, sediment thickness, and LOI 550°C. For improving the performance of the ANN, the input variables were normalized to range from 0 to 1. The Hidden layer is a processing layer where the transfer function was a hyperbolic tangent sigmoid function, in which the values transferred to the output layer will vary between − 1 to 1. The output layer is the gas content per unit area in L m⁻².

We used the Levenberg—Marquardt backpropagation for training of the ANN. During the training step, the algorithm randomly divided the data into three parts. 70% of the data points are used for the actual training of the neural network, 15% are used for the validation of the ANN during the training calculations to avoid overfitting, and the remaining 15% are not included for the training and are used for testing the trained model.

The model’s result for gas content prediction were evaluated considering the coefficient of determination (R ²) between observed and predicted gas content and through the relative error. The relative error, which was expressed in percentage, was calculated for each grid cell as E r r o r r e l a t i v e =   ( | y p r e d i c t e d − y o b s e r v e d | y o b s e r v e d ) × 100 where y p r e d i c t e d and y o b s e r v e d are the predicted and estimated gas content, respectively.

Lastly, the gas content in the sediment derived from the hydro acoustic surveys performed in 2016 and in 2019 was compared to analyze the temporal changes in different regions of the reservoir and to check the application of the prediction model from 2019 against the measurements of 2016. The temporal change in gas content between 2016 and 2019 was tested for each valid grid cell using a non-parametric hypothesis test (Wilcoxon rank sum test), considering only cells that contained 30 or more pings (sound pulses), which resulted in 1,321 cells for the comparison.

The highest values of potential methane production for all three sampling locations were found in the top 10 cm sediment layer. The maximum PMP values ranged from 3.4 mgCH₄ L⁻¹ d⁻¹ at P3 to 5.9 mgCH₄ L⁻¹ d⁻¹ at P2. From the near-surface layer, PMP decayed approximately exponentially with increasing depth in the sediment (Figure 3A). The integrated temperature corrected PMP for the top 10 cm resulted in a potential methane flux at the water-sediment interface of 306.4, 450.2, and 91.1 mgCH₄ m⁻² d⁻¹ at locations P1, P2, and P3 respectively.

Dissolved methane concentration was lowest in the overlying water with a strong increase in concentration towards the water-sediment interface (Figure 3B). In the sediment, the dissolved methane concentration continued to increase with increasing depth with comparable vertical gradients at all three locations. The maximum porewater methane concentrations were measured in the deepest DPS chamber at 25–30 cm sediment depth (1.4 mmol L⁻¹ at location P1, 1.3 mmol L⁻¹ at P2, and 1.2 mmol L⁻¹ at P3).

Gas ebullition flux was continuously measured during 45 days, starting prior to the acoustic survey conducted in 2019. The largest amount of gas was collected by the bubble trap at location P1 (9.3 L), followed by locations P2 (6.0 L), and P3 (3.7 L) (Figure 3C). For the 45 measurement days, the recorded volume represents a methane ebullition flux of 118.1, 77.5, and 48.1 mgCH₄ m⁻² d⁻¹ at the locations P1, P2, and P3 respectively (with CH₄ fraction in bubbles collected from the three locations of 68.9 ± 6.8%).

At locations P2 and P3 (near to the water intake and in the dam region) the dynamics of gas accumulation is characterized by stepwise increases, in which periods of several days without ebullition are interrupted by ebullition events. At P1 in contrast, more continuous ebullition events occurred throughout the measurement period, as it can be observed by the prolonged rise of the accumulated gas volume curve (Figure 3C). At all locations, pronounced ebullition was observed on February 5th during the period when the echo-sounding surveys were conducted. This ebullition event was associated with a weather change; from February 04th to February 07th, there was a reduction of the air temperature by approximately 10°C in comparison to the previous days (Supplementary Figure S3). Atmospheric pressure strongly decreased starting from February 05th. During the surveys, the water level in the reservoir decreased by 2 cm (from February 05th to February 07th) and the mean wind velocity was 1.6 ± 0.4 m s⁻¹ (average ± standard deviation).

According to Sotiri et al. (2021), the reservoir bottom is overlaid by an unconsolidated fine-grained low-density material layer. The analysis of the sediment cores showed that the bottom sediment is dominated by silt-clay grain sizes and an average loss on ignition (LOI) of 17 ± 8.5%. The averaged LOI in the sediment estimated from the acoustic parameter attack (Figure 4C) was 14.7 ± 4.9% for the whole reservoir, in which 3.9% of the valid grid cells had LOI of less than 10% and a minimum value of 1.8%.

The sediment magnitude estimated by Sotiri et al. (2021) varied from 0 to a maximum of 1.8 m, with highest sediment accumulation in the upstream region near to the river inflow and in the region near to the dam (see Figure 4B), where the water depth varies from 10 to 15 m. Average sediment thickness in the analysis grid ranged from 0 to 1.5 m, with a mean value of 0.5 ± 0.2 m.

The overall mean value of the maximum acoustic backscatter in the analyzed grid cells was −6.6 ± 2.0 dB. According to the empirical relationship proposed by Anderson and Martinez (2015), this corresponds to a mean sediment gas content for the whole reservoir of 4.6 ± 3.2 L m⁻². The largest values of sediment gas content were estimated for the upstream region of the reservoir, whilst the smallest values were found near the banks and in the deepest region of the reservoir in front of the dam (Figure 4D). Elevated (above average) gas content was also estimated for the central part of the deeper region near the dam. At the locations where the automated bubble traps were deployed, the averaged gas content in the sediment was 6.2 ± 2.1 L m⁻², 4.1 ± 1.6 L m⁻², and 5.6 ± 2.0 L m⁻², for the P1 to P3 respectively (Table 2).

The variables water depth, which ranged from 1.4–15.35 m (8.9 ± 3.3 m), sediment thickness with values in the range of 0.03–1.5 m (0.5 ± 0.2 m), and LOI 550°C varying from 1.8–53.5% (14.7 ± 4.9%) were tested as predictor variables for the gas content in the sediment, which varied from 0.1 to 40.4 L m⁻² (4.6 ± 3.2 L m⁻²). In an exploratory analysis of the variables, a Spearman rank correlation test was applied to check statistical correlation among the parameters (see Supplementary Figure S4). Considering the individual correlations between gas content and the predictor variables, the strongest correlation was found between gas content and LOI (Spearman correlation r s = − 0.58 and p = 0), whereas weaker, yet significant, correlations were found between gas content and water depth ( r s = − 0.12 , p < 0.05 ) and gas content and sediment thickness ( r s = − 0.18 p < 0.05 ).

Multiple linear regression (model MLR I) resulted in a coefficient of determination (R ²) of estimated gas content of 0.24 (p < 0.05 for the F-test on the model). As the distribution of gas content in the sediment deviates from a normal distribution (Supplementary Figure S4), a multiple linear regression with log transformed gas content (MLR II) was tested, resulting in a R ² of 0.16 (p < 0.05 for the F-test on the model). The stepwise multiple regression (SMR) resulted in a R ² of 0.54 (p = 0 for the F-test on the model). The full model equations are presented in Table 3. The trained artificial neural network (ANN) reproduced 89% of the gas content variance (R ² = 0.89 for all data points). Comparable coefficients of determination were obtained for the test data set which was not included during the training of the ANN (R ² = 0.91, Supplementary Figure S5).

The largest relative errors (>80%) for all models occurred near the banks at deepest region near to the dam (Figure 5 panels (B), (C), and (D)). The averaged relative error for the 4,651 grid cells was 42.5, 61.2, and 72.96% for the ANN, SMR, and MLR I model, respectively. In addition, the linear regressions applied to the predicted and estimated gas content indicate a systematic underestimation of all three models compared to observations (Figure 5A). The frequency distributions of the relative error showed that for all models the most frequent relative errors obtained were smaller than 25% (see Supplementary Figure S6 panels (A) and (B)). On the other hand, the largest relative errors values (above 500%) occurred for the grid cells with gas content per unit area of less than 4 L m⁻².

The gas content in the sediment was also calculated with the hydro acoustic data recorded during the 2016 survey. In both surveys (March 2016 and February 2019), the water column was thermally stratified, with a warmer top layer of ∼26°C and colder bottom water with the lowest temperature of ∼21°C at location P3 (Supplementary Figure S7 panels (A) to (C)). Small variations of the reservoir water level were recorded during the hydro acoustic surveys of both years (3 cm in 2016 and 2 cm in 2019). However, for the 2016 survey the water level was ∼1 m higher than the water level recorded during the survey of 2019.

In both years, the estimated gas content in the sediment of the reservoir was similar with mean values of 4.4 ± 3.1 L m⁻² in 2016 and 4.6 ± 3.2 L m⁻² in 2019. Nevertheless, the range of the spatial variation of the gas content was slightly different with larger values in 2016 (2016: 0.1–57.0 L m⁻²; 2019: 0.1–40.4 L m⁻²). In both years, the largest gas content was estimated for sediments in the upstream part of the reservoir. Nevertheless, the coarser spatial resolution of the echo-sounding transects conducted in 2016 (partial loss of data) did not allow for capturing the spatial structure of the gas content hotspots in this region (Figure 6A).

As the spatial coverage of the acoustic transects during the 2016 survey was coarser than in the 2019 measurements, only 1,321 grid cells could be considered for comparing the gas content between both years. A non-parametric hypothesis test was performed for each grid cell to verify the occurrence of a significant difference in the median gas content between both years (Figure 6). 48% of the cells had no significant difference in the median gas content and 52% were found to have different median gas content. Significant differences between both years occurred mainly at the shallower upstream part of the reservoir (near monitoring location P1), where the gas content was lower in 2019 compared to 2016 (Supplementary Figure S8). Higher gas content was detected in the reservoir stretch between the monitoring locations P1 and P3.

Lastly, we compared the predicted gas content from the artificial neural network model, which was based on 2019 data, with estimated gas content in 2016 and 2019. The gas content in the sediment in 2016 agreed with the predicted and estimated gas content of 2019, as shown in Figure 6D (R ² = 0.45 for a linear fit with no intercept between predicted gas content from the ANN model and estimated gas content in 2016, p = 0 for the estimated coefficient). As previously mentioned, the upstream areas of potential high gas content in the sediment were not resolved in the 2016 survey. Nevertheless, grid cells with higher gas content in 2016 were observed in the vicinity of that region.

A high spatial resolution sediment gas content mapping was obtained as a function of the acoustic parameter maximum backscatter using the empirical relationship proposed by (Anderson and Martinez, 2015). The average maximum backscatter in Passaúna Reservoir ( − 6.6 ± 2.0 dB) is in the range of the values obtained by Anderson and Martinez (2015) for the eutrophic reservoir Hodges. The gas content ranged between 0.1 and 40.4 L m⁻², which for a 10 cm sediment layer represents a volumetric gas content of 0.1–40.4% (vol). The addition of gas voids to the sediment matrix results in a decrease in the sediment wet bulk density leading to instability of the sediment matrix and release of gas bubbles (Van Kessel and Van Kesteren, 2002; Liu et al., 2018). Van Kessel and Van Kesteren (2002) found that for muddy sediments, gas fraction in the range of 25–37% results in instabilities in the sediment matrix. For Passaúna reservoir less then 1% of the cells had gas content larger than 25% (Supplementary Figure S9). Moreover, considering the sediment wet bulk density mapped by Sotiri et al. (2019b), the regions with lower wet bulk density had relatively less gas content.

For lake Kinneret the upper limit of gas fraction in the sediment from acoustic measurements were 0.2% (Katsnelson et al., 2017) and 3.8% (Uzhansky et al., 2020). Nevertheless, for laboratory conditions Liu et al. (2016) found a maximum gas content in incubated clay sediment of 46.8% with a depth-average of 18.8% for this fine sediment. Thus, the few highest fractions of gas content found (<1% cells) at Passaúna are in the same range as the values reported by Liu et al. (2016), whereas the averaged values of gas content are one order the values found in Lake Kinneret (Uzhansky et al., 2020).

Transversal and longitudinal variation of gas content in the sediment were observed along the reservoir. Longitudinally, the largest amount of free gas in the sediment was detected in the upstream part of the reservoir, closest to the main inflow (Figure 4D), in a region with the largest sediment thickness and low water depth. We found that the gas content in the sediment tended to be higher in regions of preferred sedimentation (i.e. large sediment thickness), nevertheless with gas content magnitude being additionally affected by the water depth. For instance, at the deepest region of the reservoir where elevated sediment deposition was also mapped, the gas content was lower (up to 10 L m⁻²) compared to the shallow upstream region (up to 30 L m⁻², Supplementary Figure S10C).

The occurrence of temporal change in gas content in the sediment was verified by analyzing two available echo-sounding surveys conducted 3 years apart from each other. The two surveys were done during summer at comparable water temperature and thermal stratification, and with small water level variation during both surveys, see Supplementary Figure S7. The mean sediment gas content was similar during both surveys (4.4 ± 3.1 L m⁻² in 2016 and 4.6 ± 3.2 L m⁻² in 2019). The range of its spatial variability, however, differed and the maximum gas content in 2016 (56.9 L m⁻²) was almost 50% higher than the maximum value in 2019 (40.4 L m⁻²). We found that 52% of the analyzed sediment area had statistically different gas content between both surveys and lower gas content in 2019 was mainly estimated for the shallower upstream region (around location P1) (see Supplementary Figure S8). The lower gas content can be related to the 1 m higher water level in 2016. Higher water level would represent additional hydrostatic pressure at the water-sediment interface, and thus increases the burden for bubble formation. A strong linear decrease of sediment gas content with increasing water level has been observed in Lake Kinneret, where acoustic surveys were conducted across several years (Uzhansky et al., 2020). In addition to water depth, the proximity to the main inflow Passaúna river makes this region more affected by floods events which can alter the bottom sediment conditions.

Although small-scale variability of gas content in the sediment is not resolved in this study, either due to the echo response which is resulting from a bottom area (0.8 m² at 8.3 m depth) and from the averaging with the grid cell; measurements with closer distributed hydro acoustic transects were important for detecting spatial patterns in the gas content. For instance, in the 2016 survey, which covered a distance 3 times shorter than the survey performed in 2019, the areas of high gas content located at the upstream region of the reservoir were not captured due to missing data points.

The results from the incubated sediment cores showed that the top sediment layer has the largest potential for methane production (PMP). The finding that the top 10 cm sediment layer is most productive is in agreement with Isidorova et al. (2019), who considered sediment within the range of sediment age (<6–12 years) as still active for methane production. In Passaúna, the first 10 cm of sediment have an age of 5 years, based on the findings of Sotiri et al. (2021), who estimated a sedimentation rate of 1.9 cm yr⁻¹.

The integrated PMP over the 10 cm depth resulted in a potential methane flux at the sediment water interface of 306, 450, and 91 mgCH₄ m⁻² d⁻¹ at the locations P1, P2, and P3 respectively at an average sediment temperature. The PMP of Passaúna sediment is within the range found in sediments from impoundments in Germany (Wilkinson et al., 2015) and the potential fluxes at the water-sediment interface are higher than the values reported for tropical lakes in Uganda (15.4–144 mgCH₄ m⁻² d⁻¹) (Morana et al., 2020). The largest integrated PMP value at location P2 coincided with the highest organic matter content (indicated by LOI 550°C) from sediment cores (LOI 550°C 15% at P1, 22% at P2, and 16% at P3), which is in agreement with other studies that reported enhanced methane production for sediments with higher organic matter content (West et al., 2012; Grasset et al., 2021).

The averaged methane ebullition fluxes were 2.6, 5.8, and 1.9 times smaller than the corresponding potential flux at the water-sediment interface estimated from the integrated PMP at the three locations (see Table 2). The lower emission rates compared to sediment CH₄ production differ from observations by Wilkinson et al. (2015), in which the measured ebullition flux was comparable to the potential methane flux from PMP of shallow (<4 m) impoundments. We attribute this difference to the higher water depths of Passaúna Reservoir (average of 8.3 m), which can favour the diffusive flux of methane from the sediment (Langenegger et al., 2019). The transported methane can accumulate in the overlaying water where it is susceptible to oxidation or release to the atmosphere during mixing events (Vachon et al., 2019). In fact, from the dissolved methane concentrations in the water overlaying the sediment, the shallowest upstream location, which is more likely to have mixing due to its lower water depth, had methane concentrations that were one order of magnitude smaller than at the other two sites (0.02, 0.36, and 0.33 mmolCH₄ L⁻¹ for locations P1, P2, and P3).

Furthermore, although location P1 had a lower integrated PMP than location P2, the methane ebullition flux was 1.5 times higher, which we attribute in addition to the enhanced diffusion of CH₄ to the water, to the larger water depth and thus, to the total pressure at the water-sediment interface. As the solubility limit of methane in the porewater increases with higher hydrostatic pressure (Figure 3B), larger dissolved gas concentrations are required for gas bubble formation (Bazhin, 2003; Langenegger et al., 2019). As shown in Figure 1, the methane concentration measured in sediment porewater at locations P2 and P3 were lower than the methane saturation limit considering the local water depth at these two sites. On the other hand, at the shallower location P1, the methane concentration in the sediment porewater was above the saturation limit for the deeper sediment layers. Lastly, an important aspect to consider, is that the PMP was obtained for laboratory conditions which may differ from the dynamic environmental conditions in the reservoir.

Ebullition is widely characterized in the literature as being highly variable in space and in time (Wik et al., 2013; Maeck et al., 2014). For Passaúna Reservoir, it was previously observed that ebullition events at the three monitoring locations where synchronized on a daily basis in which the ebullition events were triggered by large-scale forcing, such as drops in the atmospheric pressure (Marcon et al., 2019). For the 45-days time-series analyzed in the present study, a larger release of gas from the sediment was observed for February 4th to February 7th, which we could associate with the beginning of a period of decreasing atmospheric pressure (Supplementary Figure S3). We suggest that the drop in atmospheric pressure during the echo-sounding survey in 2019 potentially triggered gas release from the bottom sediment. Nevertheless, considering the elevated potential methane production of the bottom sediments, we assumed that the spatial patterns of gas in the sediment estimated from the echo-sounding was valid, as the ebullition trigger, in this case the atmospheric pressure drop, acted over the whole reservoir area.

Heterogeneities in gas distribution in the sediment are caused by different mutually influencing factors. Methane production in the sediment is driven by organic matter supply (Grasset et al., 2021) and its degradability (Sobek et al., 2012; West et al., 2012; Praetzel et al., 2019), as well as by temperature (Aben et al., 2017; Wilkinson et al., 2019a) and the presence of alternate electron acceptors (Bastviken, 2009). Nevertheless, not all the produced methane escapes the sediment as ebullition, as also observed at Passaúna. The methane accumulation in the sediment is affected by CH₄ oxidation (Bastviken, 2009; Martinez-Cruz et al., 2018) and its transport out of sediment by diffusive fluxes (Langenegger et al., 2019). Once the methane concentration in porewater reaches saturation gas voids are formed. The distribution and persistence of gas voids is further dependent on sediment properties such as grain size (Boudreau et al., 2005; Algar and Boudreau, 2009) and the sediment capacity to hold the free gas (Van Kessel and Van Kesteren, 2002; Liu et al., 2016). The ebullition flux (i.e., release of free gas from the sediment matrix as gas bubbles), is then a result of triggers facilitating gas release provided that there is free gas accumulated in the sediment. Drops in hydrostatic (Scandella et al., 2011; Maeck et al., 2014) and atmospheric pressure (Casper et al., 2000; Natchimuthu et al., 2016) and bottom currents (Joyce and Jewell, 2003) have been reported as ebullition triggers.

The hotspot of gas content acoustically detected at the upstream region of Passaúna Reservoir indicates a higher potential for methane ebullition, which is confirmed by the highest ebullition flux recorded at the location P1. This is in accordance with numerous other studies, reporting high CH₄ fluxes in regions near to the main inflow river with high sedimentation rates (DelSontro et al., 2011; Grinham et al., 2018; Hilgert et al., 2019a; Linkhorst et al., 2021). Whereas in the deeper locations of the reservoir (monitoring sites P2 and P3), higher methane partial pressure is required for bubble formation, Figure 3B, which combined with the deposition of finer sediment particles may increase sediment cohesivity and capacity to hold the produced gas in the sediment, and thus, would explain the observed dynamics of cumulative ebullition fluxes (Figure 3C) at locations P2 and P3 with longer periods (days) of no ebullition.

The capability of predicting gas content in the sediment is a useful tool for estimating methane ebullition from inland water. Large parts of the spatial variations in the estimated sediment gas content in Passaúna Reservoir could be explained by variations in more readily accessible characteristics of the reservoir and its sediment, including water depth, sediment thickness, and organic matter content of the sediment. The latter was estimated from the attack phase of the bottom echo.

Considering the relationships between gas content and individual parameters, the gas content was positively correlated with sediment thickness ( r s = 0.2 p = 3.4   10 − 36 ) which may serve as a proxy for the deposition rate of fresh sediment (Supplementary Figure S4). Water depth was negatively correlated with gas content ( r s = − 0.1 p = 4.7   10 − 17 ), and accounts for hydrostatic pressure, as well as an indirect measure of sediment temperature (colder at largest depths), dissolved oxygen concentration in the hypolimnion (formation of bottom anoxic layer is favored during periods of stratification). Surprisingly, LOI was negatively correlated to gas content ( r s = − 0.6 p = 0). However, the correlation varied with water depth (Supplementary Figure S10) and the negative correlation was mainly caused by data from the upstream region of the reservoir, where large sediment thickness with high gas content occurred together with lowest LOI values (<10%) were found. Gas content and LOI were acoustically derived. Although they are calculated from different parameters (maximum backscatter and attack), there is a strong dependency between maximum backscatter and attack (Spearman correlation r s = 0.9 p = 0). On the other hand the acoustic derived LOI values agreed with the measurements from sediment samples (Sotiri, 2020). In addition, in this study the available data set didn’t allow to verify the relationship of LOI and organic matter content for the sediment characteristics of Passaúna Reservoir.

As discussed above, the gas content in the sediment depends on a combination of different parameters, thus the combination of available information was tested as predictors of gas content. The multiple regression (MR) models resulted in a lower agreement between predicted and estimated gas content in comparison to the artificial neural network (ANN) model (R ² < 0.55 for MR and R ² = 0.89 for the ANN). The ANN model has the capability of accounting for nonlinearities among the variables and to handle high-dimensional multi-scale systems, thus identifying hidden patterns in the data set (Fausett, 1994). In the present application, this was observed in the contrasting magnitude of gas content for comparable sediment thickness regions with differing water depths.

The largest relative errors between the predicted and estimated sediment gas content were found for low gas content (<4 L m⁻²), at the deepest region of the reservoir near to the dam, and towards the reservoir banks. Steep slopes are known to affect acoustic backscatter measurements of bottom sediments (Sternlicht and de Moustier, 2003). This supports the application of a slope threshold in our spatial analyses. On the other hand, no evident dependence of the relative error on the average slope of the respective grid cell was observed for slopes smaller than the threshold (Supplementary Figure S6C). In addition, compact sediments are reported to have higher maximum backscatter comparable to the acoustic response of gassy sediments (Hilgert et al., 2016; Sotiri et al., 2019a).

We contrasted the predicted gas content from the trained artificial neural network for 2019 with the estimated gas content derived from the 2016 and 2019 hydroacoustic surveys and found good agreement between both years (Figure 6D). In this way, we denote that even though small spatial scale heterogeneities occur, which were averaged within grid cells, the main underlaying spatial variability of gas content was maintained between the years and that the trained artificial neural network (ANN) model for gas content prediction is valid for both years. The prediction of the ANN model can be complemented and tested further with a stricter slope criterion and with the inclusion of additional relevant parameters. For instance, the origin of the organic matter (West et al., 2012), sediment exposure to dissolved oxygen concentrations (Yvon-Durocher et al., 2014), and potential of methane production in the sediment (Wilkinson et al., 2015) were reported to be relevant to methane ebullition. However, mapping of such parameters for the whole reservoir would require the combination of the in-situ measurements complemented with modelling or upscaling techniques.

In this study, we derived gas content in the sediment from acoustic measurements and investigate its spatial-temporal variability in a subtropical reservoir. The analysis was supported by comparing gas content estimates with spatial maps of sediment thickness, loss on ignition, and bathymetry and considering the three locations with estimated potential methane production, dissolved methane concentration in the pore water, and continuously measured ebullition flux.

The potential of using the echo-sounding approach for detecting gas content in the sediment and the need for further investigations of its spatial distribution and relation to methane fluxes was highlighted in previous studies (Anderson and Martinez, 2015; Katsnelson et al., 2017; Uzhansky et al., 2020) and corroborated in this study. One remaining challenge is the lack of direct measurements of sediment gas content under in-situ conditions that can serve for testing and calibration of acoustic approaches. Anderson and Martinez (2015) collected gas that was released from the sediment upon mechanical disturbance, which can be difficult when applied from a boat at larger water depth where an accurate definition of the disturbed area that contributes to the collected volume of gas is not possible. More recent studies analyzed sediment cores frozen under in-situ conditions [freezed cores (Dück et al., 2019a)], which are analyzed in an X-ray CT scanner to quantify the amount of free gas (Dück et al., 2019b; Liu et al., 2019). However, freezing of the sediment cores is also reported to cause mechanical disturbances in the sediment core and result in additional bubble formation (Dück et al., 2019b). Sampling of pressurized sediment cores (e.g., Wilkens and Richardson (1998)) certainly allow for most accurate estimates of column-integrated gas content in the laboratory, yet they require the support of divers for sediment sampling and are also affected by bubbles that escape during corer penetration.

The understanding of spatial variability and temporal dynamics of methane fluxes from inland waters can be improved by knowing the process that affect the production, transport, oxidation, and emission of CH₄, which include storage of free gas in the sediment. High resolution acoustic surveys can provide estimates of sediment gas content and its spatial and temporal dynamics. Nevertheless, additional sampling locations for echo-sounding and ebullition monitoring would be required to explore relationships between gas storage and ebullition. A main advantage of acoustic gas content measurements as a proxy for ebullition flux is the high potential areal coverage of echo-sounding in comparison to the limited area sampled by bubble traps (funnel diameter of 1 m) and the possibility to measure transects covering the entire reservoir. As an additional aspect, it remains to be investigated if the uncertainties of flux measurements that are associated with temporal variability of ebullition can be reduced by accessing the gas stored in the sediment.