Quantifying bubbling emission (ebullition) of methane from a rice paddy using high-time-resolution concentration data obtained during a closed-chamber measurement

Methane (CH 4 ) produced in rice-paddy soil is transported to the atmosphere either via the rice plants or by bubbling events (ebullition); however, little is known about the frequency and intensity of bubbling CH 4 emissions and the factors that affect them. We developed a method to quantify ebullition using high-time-resolution (~1 Hz) CH 4 concentration data obtained by closed-chamber measurements. Field measurements were conducted in a Japanese rice paddy at different rice growth stages: panicle formation (PF), booting (BT), and heading (HD). A dataset of 132 chamber measurements was used to develop and evaluate the method. A scripting le written in R programing language was used to automatically determine CH 4 emissions via the two pathways. Plant-mediated CH 4 emission intensity was constant during chamber deployment and was reected as a steady linear increase in chamber [CH 4 ] with time or as a constant baseline in a ux time series. We found that the plant-mediated emission could be determined as the peak with the lowest ux intensity in the ux frequency distribution even if bubbling events occurred during the chamber deployment. The eld measurement results in combination with established data processing protocols showed that at PF, ebullition contributed only 4% of the total emission, whereas it accounted for 32% and 60% of the total emission at BT and HD, respectively. In contrast, the plant-mediated ux variation among growth stages was smaller. Both ebullition and plant-mediated emissions correlated signicantly with air temperature at HD, but the magnitude of the dependency was much higher for ebullition than for rice-mediated emission. These results demonstrate that ebullition occurs more frequently than has previously been thought, and the different transport pathways show varying degrees of dependency on plant phenological and environmental factors, thus underscoring the need to separately determine CH 4 emissions via each transport pathway.


Introduction
Rice paddies are a globally important source of atmospheric CH 4 (Ciais et al., 2013), a potent greenhouse gas and historically the second largest contributor to global warming after CO 2 (Myhre et al., 2013).
Methane produced in paddy soils enters the atmosphere either through the aerenchyma tissue of the rice plants (Nouchi et al., 1990;Wang et al., 1997) or by ebullition of CH 4 -containing gas bubbles (Komiya et al., 2020;Schütz et al., 1989;Wassmann et al., 1996). Entrapped gas bubbles (rather than dissolved CH 4 in soil solution) in rice paddies have been shown to represent a major CH 4 inventory component at later rice growth stages, even in soils regarded as water saturated (Tokida et al., 2013). Many studies have shown a very high CH 4 mixing ratio in gas bubbles in rice-paddy soils (Byrnes et al., 1995;Uzaki et al., 1991;Watanabe et al., 1994). Molecular diffusion of dissolved CH 4 across the water-atmosphere interface can also occur, but that contribution is usually marginal (Butterbach-Bahl et al., 1997;Schütz et al., 1989, Komiya et al., 2015 because CH 4 is a sparsely soluble gas and its diffusion in water is four orders of magnitude smaller than that in the gas phase. In addition, methanotrophic bacteria may consume 80-90% of the CH 4 diffusing through the soil-water interface before it reaches the atmosphere (Banker et al., 1995;Frenzel et al., 1992).
Quantifying CH 4 ebullition separately from the total CH 4 emission remains a challenging issue.
Conventionally, the ebullition rate has been estimated by trapping evolving gas bubbles in an inverted water-lled funnel placed over the soil surface (Holzapfel-Pschorn et al., 1986) or by using a chamber method where, to exclude rice-mediated emission, no rice plants are included in the chamber (Cheng et al., 2006;Komiya et al., 2015). These techniques intrinsically assume spatial homogeneity of bubbling emissions; however, more ebullition might occur in proximity to rice plants because they provide carbon sources for CH 4 production via rhizodeposition and degradation of plant tissues (Tokida et al., 2011).
Thus, the currently accepted notion that rice-plant mediated transport is the dominant emission pathway is open to question, especially when the rice plants are still small (Schütz et al., 1989;Wassmann et al., 1996), or at later growth stages, when root system decay may reduce the ability of rice to transfer CH 4 (Tokida et al., 2013).
It is obvious that even when ebullition is to be measured, it is preferable to include rice plants in the measurement chamber, but with conventional gas sampling and analysis protocols, it is not possible to separate ebullition from plant-mediated transport of CH 4 . Recently, spectroscopic gas analyzers have been developed that can measure the CH 4 concentration ([CH4]) with high precision and high time resolution (~1 Hz) (e.g. Gogo et al., 2011). Tokida (2021) has developed a modi ed closed chamber method for measuring rice-paddy CH 4 emission that uses a commercially available portable spectroscopic gas analyzer. With this method, ebullition is identi ed by a "jump" in the high-timeresolution record of [CH 4 ] in the chamber air. Motivated by these advancements in measurement techniques, in this study we aimed to develop systematic data handling and analytical procedures to distinguish CH 4 ebullition from total emissions. We also investigated two factors potentially affecting ebullition: i) plant growth stage, and ii) temperature, i.e., an environmental factor well known to affect CH 4 emission.

Experimental eld and rice cultivation
We used an experimental paddy eld located at the National Agriculture and Food Research Organization in Japan (36°01'28" N, 140°06'30" E). We sowed pregerminated seeds of rice (Oryza sativa L., cv. Koshihikari) on April 15, 2020, in seedling trays, and raised the seedlings in an open eld. Koshihikari seedlings (5-leaf age) were transplanted at a spacing of 15 cm × 30 cm with three seedlings per hill on May 14. We prepared a total of 44 measurement plots (11 plots × four replicates) within a single eld (10 m × 50 m), where each plot was an independent canopy consisting of 24 hills of rice plants (1.2 m × 0.9 m). For each replicate, 11 plots were placed every 5 m (center to center) along the longer side of the eld, and the four replicates were placed every 2.3 m across the shorter dimension of the eld. The eld was continuously ooded from transplanting until the middle of the grain-lling stage. Heading initiated on July 26. In June, July, and August 2020, the monthly average air temperature was 22.4, 23.5, and 27.8 °C, respectively, and monthly precipitation was 199, 229, and 59 mm, respectively. Air temperature (2 m above the ground level) was observed at a meteorological tower 200 m away from the experimental eld.

Measurement of CH 4 emission
Methane emission was measured successively in each plot by a closed-chamber method (Minamikawa et al., 2015) but with a substantial modi cation of the [CH 4 ] measurement procedure. [CH 4 ] was monitored on site at high time resolution (at intervals of 0.9 s) by using a mobile gas analyzer (G4301, Picarro Inc., Santa Clara, CA, USA) instead of a gas chromatograph in the laboratory (Tokida, 2021). The chamber and analyzer were connected by two PTFE tubes; a chamber air sample was withdrawn from the chamber via one tube, and then, after it passed through the gas analyzer, it was re-injected into the chamber via the other tube. We used two chambers and adopted an alternating chamber strategy to make the most e cient use of the gas analyzer (Tokida 2021). A closed-top acrylic chamber (basal dimensions, 30 cm × 60 cm) was used to enclose four hills of rice plants at the center of each plot. Measurements were made over four days (1 replicate per day) at each of three growth stages: panicle formation (PF, leaf number index 84-87; June 29, 30, July 2, 3 for Reps 1, 2, 3, 4, respectively), booting (BT, ag leaf just fully expanded; July 21, 22, 24, 25), and full heading (HD; August 4, 5, 7, 8). Thus, a total of 132 measurements (11 plots × 4 replicates × 3 stages) were obtained during this study. At the PF stage, only a 60-cm-high chamber equipped with a fan was used, whereas at the BT stage, a double-deck chamber, consisting of a 20-cm-high lower chamber without a fan and a 60-cm-high upper chamber equipped with a fan was used.
At the HD stage, the lower chamber was replaced with a 60-cm-high chamber equipped with a fan. Actual duration of the chamber closure differed from one measurement to another (1.6-34 minutes). We kept the chamber closed until we obtained a steady [CH 4 ] increase, not interrupted by the bubbling event, for at least 1 minute. Therefore, the chamber deployment time was basically longer as the bubbling event occurred more frequently (section 2.3).

Separate quanti cation of CH 4 emissions via rice plants and ebullition
We estimated plant-mediated emission and ebullition separately by analyzing the high-time-resolution data of chamber [CH 4 ] (Fig. 1). Two patterns were observed in each [CH 4 ] time series: a steady linear increase, and sudden rapid increases (Fig. 2a, d). We assumed that the steady linear increase corresponded to the emission via rice plants, which occurs at a constant rate, and that the stepwise increases occurring over a short time were due to the release of bubbles, that is, ebullition (indicated by arrows in Fig. 2a, d). Ebullition was observed frequently, especially at the HD stage, and the intensity of ebullition could vary substantially from one bubbling event to another, even during a single chamber deployment. The total CH 4 ux was calculated by the traditional regression analysis approach (Minamikawa et al., 2015), whereas the rate of plant-mediated emission was determined by a series of data processing steps, starting with noise reduction (see section 2.3.2) and followed by analysis of the ux frequency distribution (section 2.3.3). Finally, the bubbling emission (ebullition) was calculated as the difference between the two (section 2.3.4). Each data processing step is explained in detail below.

Cleaning raw [CH 4 ] data
Because the gas analyzer was operated continuously, the raw time-series data contained [CH 4 ] measurements in ambient air introduced in between the measurements from different plots. Therefore, we removed the data from unstable periods (gray bands in Fig. 2a, d) by analyzing the raw [CH 4 ] and smoothed ux data (obtained by differentiation of [CH 4 ] with time and applying an adequate smoothing procedure (i.e., by applying the 10-point moving average ve times) to obtain "processed" [CH4] data ( Fig.  1). The cleaning procedure is described in detail in the Supplementary Information.

Moving-average protocols for noise reduction
The plant-mediated emission, which is represented by a linear increase in [CH 4 ] (Fig. 2a, d), can be identi ed as the " at background" when expressed as a time series of d [CH 4 ]/dt or in units of ux (Gogo et al., 2011). Therefore, we differentiated [CH 4 ] measured at 0.9 s intervals with respect to time and carried out the calculations necessary to transform the results into CH 4 ux (mg-C m -2 h -1 ) data using the air volume and temperature in the chamber (Minamikawa et al., 2015) (Fig. 1: "Convert to ux" on the right under "Processed [CH4] data"). However, the time-series ux data (at 0.9 s intervals) showed frequent uctuations (Fig. 2b, e, light gray lines), indicating that the high-frequency measurements included substantial errors due to incomplete mixing of chamber air and to the analytical error of the gas analyzer. Therefore, to correctly determine the plant-mediated ux, we needed to reduce the noise (by smoothing) (Fig. 1: "Moving average" under "Convert to ux").
To smooth the data, we calculated moving averages. To identify the most appropriate protocol, we compared moving-average results obtained by applying different numbers of data points (5, 10, 20, or 30 points) and repetitions (1-9). As a result, we obtained 36 different smoothed ux datasets for each chamber measurement. Then, we computed the frequency distribution of each ux dataset by estimating the kernel density distribution and using it to determine the plant-mediated ux for that dataset (Fig. 2c, f) (see section 2.3.3). We then investigated the appropriateness of each moving-average protocol by examining the robustness of the plant-mediated ux (Fig. 2b, e, blue baseline, and Fig. 2c, f, solid red vertical line) estimated by using each combination (number of data points in the averaging window and number of applications of the moving average; see section 3.1).

Plant-mediated emission: estimation from the kernel density distribution
To determine the intensity of the plant-mediated ux (blue baseline in Fig. 2b, e), we converted each ux time series into the ux frequency distribution by estimating the kernel density distribution (Fig. 2c, f). We hypothesized that the baseline emission could be determined from the peak with lowest ux intensity in the frequency distribution (solid red lines in Fig. 2c, f). We then evaluated whether the estimated ux agreed with the baseline in the ux time series by visual inspection (blue lines in Fig. 2b, e).
In the kernel density estimation, the bandwidth (smoothing parameter) was set to 0.30 for all plots; this value was the median bandwidth selected for all 132 ux measurements (11 plots × 4 replicates × 3 stages) by the bw.nrd function in R ver 4.0.3 (R Core Team, 2020). We then identi ed local maxima (i.e., peaks) in the kernel density distribution by using the peaks function (with default parameters) in the splus2R package in R (Constantine and Hesterberg, 2021).
When little ebullition occurred, the plant-mediated ux could be identi ed as the sharp peak with the highest density in the kernel density distribution (solid red vertical line in Fig. 2c). When considerable bubbling emission occurred, the kernel density distribution showed multiple peaks, but the plant-mediated ux could be still determined as the local maximum (peak) with the lowest ux intensity (solid red vertical line in Fig. 2f).

Total emission and ebullition
The total emission was determined by the traditional regression approach; we estimated the rate of increase of [CH 4 Fig. 2a, d) and

] (processed [CH 4 ]) by a linear regression analysis (dashed lines in
converted it to the units of ux (mg-C h -1 m -2 ) using the air volume and temperature in the chamber (Minamikawa et al., 2015). The rate of ebullition averaged over each chamber deployment was calculated as the difference between the total and plant-mediated emission (total emission minus plant-mediated emission). Note that we only estimated the intensity of bubbling averaged over each chamber deployment period, not that of each bubbling event.

Statistical analysis
We investigated the relationship between CH 4 uxes via rice plants and via ebullition at each growth stage and air temperature by a simple correlation analysis using the lm function in R ver 4.0.3.

Results And Discussion
3.1 Selection of the moving-average protocol and estimation of plant-mediated CH 4 emission Comparison of the results obtained by using different smoothing protocols led us to conclude that application of a 10-point moving average ve times was the most appropriate protocol for our dataset. The most stable estimates were obtained with a 10-point window; standard deviations of estimates obtained by applying moving averages with different window sizes one to nine times increased in the order 10, 20, 5, and 30 data points (the median standard deviation for the 132 ux measurements was 0.05, 0.09, 0.29, and 0.35 mg CH 4 m -2 h -1 , respectively).
The estimated plant-mediated ux generally increased with the number of times the moving average was applied, but with a 10-point window, the estimated values reached a plateau with two applications for most plots and with ve applications for all plots except one (Fig. 3, 10-point results). In contrast, adopting other averaging windows caused the values to increase continuously with the number of applications, especially at the HD stage (Fig. 3). Thus, with a 5-point moving average, the noise might not be su ciently removed, whereas using a 20-or 30-point moving average might obscure uxes resulting from bubbling events. Visual inspection con rmed that the uxes estimated by using a 10-point moving average ve times agreed well with the baselines of the ux time series (blue lines in Fig. 2b, e; see Fig. S2 for the results for all plots), indicating the validity of this protocol.
Note that the protocol for the moving average reported here may be speci c to the current dataset. The best smoothing protocol may depend on many factors, including the analytical error of the gas analyzer, the chamber design, the e ciency of chamber air mixing by the fan, and the chamber deployment time. Therefore, investigators need to identify the best settings for themselves. Here, we argue only that with an appropriate window size and number of repeated applications of the moving average, the plant-mediated ux can be accurately determined.

Temporal changes in the CH 4 concentration during measurement
The pattern of the [CH 4 ] increase in the measurement chamber changed as the rice plants grew (Fig. S1).
The [CH 4 ] values increased linearly at the PF stage, whereas sudden rapid increases, superimposed on the steady linear increase, were observed in several plots at the BT stage and in most plots at the HD stage, indicating more frequent ebullition events.

Contribution of ebullition to total emissions
The eld measurement results in combination with the established protocols for noise reduction (applying a moving average) and ux determination showed that total CH 4 emissions were almost the same between the PF and BT stages but increased drastically from the BT to the HD stage (Fig. 4). Little ebullition occurred at the PF stage, accounting for only 4% of the total emission, but ebullition increased as the rice plants developed and accounted for 32% and 60% of the total emissions at the BT and HD stages, respectively. In contrast, the plant-mediated ux marginally differed among the growth stages investigated. A rough estimation of the seasonal CH 4 emissions indicated that bubbling emission contributed 26% of the total emission during the growing season (see Supplementary information 3, Table S1). This value is a conservative estimate because we assumed that no bubbling emission occurred at night, but in fact, bubbling emission can occur at night, at least to some extent (Komiya et al., 2020). Therefore, ebullition contributed substantially to the total seasonal emissions and occurred more frequently than has previously been inferred (Cicerone and Shetter 1981;Nouchi et al., 1990;Butterbach-Bahl et al., 1997).

Effect of temperature on plant-mediated emission and ebullition
Both plant-mediated and bubbling emissions correlated signi cantly with air temperature at the HD stage (P < 0.001 for both) but no correlation was observed at the PF and BT stages (Fig. 5). These observations are consistent with previous studies of CH 4 emissions from rice paddies in the same region of Japan, which reported the occurrence of diel variation in CH 4 emissions only during the late reproductive stage (Minamikawa et al., 2012;Tokida et al., 2014). The positive correlation observed at the HD stage may be due to both the accumulation of trapped gas bubbles in the soil and an increase in the gas volume associated with increased air temperatures. The higher the temperature, the lower the water-solubility of CH 4 would be, leading to an increase in the gas phase in the soil. Simultaneously, the gas phase would expand in accordance with Charles's law, leading to the release of bubbles that, because of buoyancy, cannot be retained in the soil matrix (Tokida et al., 2009).
An increase in gas-phase volume may also enhance plant-mediated emission. The greater the bubble volume, the larger the root surface area in direct contact with bubbles would be, leading to an enhanced diffusive uptake of CH 4 via the root and aerenchyma tissues of rice (Tokida et al., 2013). This possible mechanism might be re ected in the weak but signi cant correlation between plant-mediated emission and temperature at the HD stage (r 2 = 0.26, P < 0.001), but the temperature dependency of plant-mediated uxes was much lower than that of bubbling emission, as shown by the smaller slope of the former (0.33 vs 1.63) in the correlation plot (Fig. 5).

Conclusion
We developed a method to separately quantify CH 4 emissions via rice plants and ebullition in a high-timeresolution dataset of [CH 4 ] obtained by a closed-chamber method. Plant-mediated emission showed a small growth-stage dependency, whereas bubbling emissions drastically increased from the BT to the HD stage. The temperature dependency of the CH 4 ux was clear for both pathways at the HD stage, but the temperature sensitivity of bubbling emissions was much higher. These different responses between the two pathways clearly indicate the need to determine pathway-dependent emissions for accurate observation and modeling of rice-paddy CH 4 emissions.   replicate 2, HD stage). In (a) and (d), data in the shaded bands, which correspond to non-deployment and unstable measurement periods, were not used in the subsequent data processing. The dashed lines represent the tted linear regression lines used to calculate the total CH4 ux. In (b) and (e), light gray, red, and blue lines indicate the original ux, the smoothed ux (by applying a moving average ve times), and the plant-mediated ux, respectively. Dark gray lines, which indicate a less-smoothed ux curve obtained by applying the moving average only once, can be useful for identifying bubbling emissions (black arrows). In (c) and (f), vertical solid and dotted red lines indicate local peaks in the kernel density distribution of the ux; the plant-mediated ux was identi ed as the peak with the lowest ux intensity (red solid line).

Figure 3
Plant-mediated uxes estimated by using different moving-average protocols: averaging windows of 5, 10, 20, or 30 data points and one to nine repeated applications of the moving average. Rep1 to Rep4 refer to the experimental replicates, each comprising 11 plots, indicated by the different colors, in the paddy eld. The vertical red dashed lines indicate the protocol selected for use in this study (10-point moving average applied ve times).