Estimating internal dissolved methane loading in rivers using a mass balance approach

Study site and sampling

The Kokai River originates from Kokaigaike Pond (elevation: 140 m) in Tochigi Prefecture, flows through Tochigi and Ibaraki Prefectures, and eventually joins the Tone River. The river has a catchment area of 1,043 km² and a main channel length of 112 km. Land use in the Kokai River basin in 2014 was predominantly mountainous areas (51%), followed by farmland (46%) and residential areas (3%) (National Land Numerical Data, https://nlftp.mlit.go.jp).

Sampling was conducted at four bridges: Nagamine (36.11132°N, 140.00009°E), Shinfukurai (36.07801°N, 140.00464°E), Joso (36.02970°N, 140.02465°E) and Yamato Bridges (N36.01978°N, 140.00563°E) (Fig. 2). We focused on two short reaches, Nagamine Bridge to Shinfukurai Bridge (upstream reach) and Joso Bridge to Yamato Bridge (downstream reach). We applied a mass balance method to determine the internal methane loading within each reach. The flow distances of these reaches were 4.37 and 2.15 km, respectively, and the downstream time (calculated as flow distance divided by mean flow velocities) ranged from 1.63 to 14.0 h for the upstream reach and 0.641 to 2.47 h for the downstream reach. The riverbed areas of these reaches were 0.315 and 0.139 km². The Fukuoka weir (36.04287°N, 140.02682°E) is upstream of Joso Bridge. The weir is operated for agricultural irrigation during the spring and summer months (April to August), and remains open during the rest of the year. The irrigated area is about 2,800 ha, the maximum water withdrawal is about 13.6 m³ s^–1, and the water storage capacity is 2.75 million tons.

River water was collected using a stainless-steel beaker with a handle, attached to a rope from the bridges. In the water samples, we measured nutrients (total nitrogen: TN, total phosphorus: TP, total dissolved nitrogen: TDN, total dissolved phosphorus: TDP, nitrate+nitrite: NOx, nitrite: NO₂⁻, ammonium: NH₄⁺, and phosphate: PO₄³⁻), dissolved organic carbon (DOC), dissolved methane concentration (CH₄), methane oxidation rate and bacterial production (BP). We measured water temperature, conductivity, dissolved oxygen (DO), chlorophyll a concentration, and phycocyanin concentration using a YSI ProDSS water quality probe.

Water samples for dissolved nutrients (TDN, TDP, NOx, NO₂⁻, NH₄⁺, and PO₄³⁻) and DOC were immediately filtered through 0.45-µm syringe filters (Durapore, Millipore, Burlington, MA, USA) on the sampling site, and the filtrate was frozen for later analysis.

For measurement of BP, 10 mL of collected river water was poured into 15-mL centrifuge tubes and incubated with 50 nM of [¹⁵N₅]-2’-deoxyadenosine (¹⁵N-dA, NLM-3895-PK, Cambridge Isotope Laboratories, Inc., Tewksbury, MA, USA) in the dark at in situ temperature for 1–5 h depending on water temperature, according to a previous study (Tsuchiya et al., 2015, 2021). The incubation was started at each sampling station immediately after the sampling. The incubation was quenched by adding 99.5% ethanol (Special Grade, Wako, Osaka, Japan) to the sample (final concentration, >20%). After quenching, each water sample was filtered onto a 0.2-µm PTFE membrane filter (Omnipore, Millipore, Burlington, MA, USA), and the filter was stored at –20 °C until further analysis.

For measurement of dissolved methane concentration and methane oxidation rate, we prepared four 20-mL glass vials (GL science) of each station; two out of the four were for T = 0 samples (ambient concentration), and the others were for T = 24 and T = 48 samples in methane oxidation rate measurements. After the sampling, the water samples were immediately transferred to the four glass vials from the stainless-steel beaker. Vials were slowly filled to overflowing with several volumes of water and sealed without any headspace. We added 0.1 mL of 8 mol L^–1 potassium hydroxide (KOH) solution (Wako) to two glass vials of T = 0 for sample preservation (Magen et al., 2014). Each glass vial was stoppered with a butyl rubber septum and sealed with an aluminum crimp seal. After sealing, the vials were incubated in the dark. During the field sampling, they were kept in a cooler box filled with ambient river water to maintain in situ temperature conditions. After returning from the field, the incubation was continued in a laboratory incubator set to the in situ river water temperature. After 24 and 48 h, we added 0.1 mL of 8 mol L^–1 potassium hydroxide (KOH) solution to the other two vials to quench the methane oxidation.

We used a compact fish sonar to measure river water velocity and depth (Deeper CHIRP +2, Baltic Vision Co., Ltd., Tokyo, Japan). We cast the sonar attached to a fishing rod from the bridge at least three lines in each station and estimated the average flow velocity from GPS information. Depth was measured by moving the sonar along the bridge to determine the cross-sectional area of the river. The discharge was calculated by multiplying the average water velocity and cross-sectional area of the river. Note that the sonar device cannot measure water depths shallower than approximately 15 cm; therefore, extremely shallow shoreline areas were not included in the discharge calculation. However, such areas were spatially limited in our study sites and are expected to have a minimal influence on overall discharge due to their typically low flow velocities (see supporting information, Fig. S1).

Sample analysis

Nutrients were analyzed in our laboratory using a continuous-flow analyzer (QuAAtro, BLTEC) in technical triplicates (Nojiri, 1987; Otsuki et al., 1993). DOC measurements were conducted as nonpurgeable DOC with a TOC analyzer (TOC-V, Shimadzu) equipped with a Pt catalyst on quartz wool. At least three measurements were made for each sample, and analytical precision was typically less than 2%. Potassium hydrogen phthalate (Kanto Chemical) was used as a standard.

To assess bacterial production (BP), bacterial DNA was isolated from the filter samples using the Extrap Soil DNA Kit Plus ver. 2 (J-Bio21, Nippon Steel & Sumikin EcoTech Corp., Chiba, Japan). Cell disruption was achieved through bead beating with a Fast Prep FP120 Cell Disrupter (MP Biomedicals, Santa Ana, CA, USA) at a speed setting of 6.0 for 40 s. DNA was then purified with magnetic bead separation, following the manufacturer’s protocol. Based on prior validation (Tsuchiya et al., 2019), DNA extraction efficiency was assumed to be 100%. Subsequent quantification of ¹⁵N-labeled deoxyadenosine (¹⁵N-dA) incorporation followed the procedures described in previous studies (Tsuchiya et al., 2015, 2020). Briefly, extracted DNA was enzymatically digested into nucleosides using three enzymes (nuclease P1, Wako; phosphodiesterase I, Worthington Biochemical Corp.; and alkaline phosphatase, Promega Corp., Madison, WI, USA). After the enzymatic hydrolysis, the amount of ¹⁵N-dA (¹⁵N₅-dA + ¹⁵N₄-dA) incorporated during the incubation was analyzed by LC–MS/MS using ¹³C₁₀¹⁵N₅-deoxyadenosine (CNLM-3896-CA, Cambridge Isotope Laboratories, Inc., Tewksbury, MA, USA) as a surrogate (internal standard), with analyses performed in technical duplicate.

The dissolved methane concentration was measured with a gas chromatograph with a flame ionization detector (GC-FID). Before GC injection, we made a headspace of 2 mL helium gas in the sample glass vials and shook them vigorously for 2 min for gas-liquid equilibrium. We injected the 0.5 mL headspace gas to GC-FID. The dissolved methane concentration was then calculated according to a previous study (Magen et al., 2014) using the Bunsen coefficient β for methane at known pressure, temperature, and salinity (Yamamoto, Alcauskas & Crozier, 1976). The relative percent difference (RPD, %) between duplicate dissolved methane measurements at time zero (T = 0) was calculated to assess analytical precision, using the following equation:

(1) $$RPD\left(\mathrm{\%}\right)=\frac{|C_1-C_2|}{\left(C_1+C_2\right)/2}\times 100$$where C₁ and C₂ represent the two replicate values. In the present study, the average RPD was 5.2 ± 4.9%, indicating acceptable analytical precision.

Calculation: methane oxidation rate, emission to atmosphere, and mass-balance model

The specific methane oxidation rate was calculated by linear regression of the natural log of methane concentration against time (0, 24, and 48 h). The specific methane oxidation rate is the first-order rate constant for methane oxidation (in units of h^–1). The volumetric methane oxidation rate (nmol L^–1 h^–1) was calculated by multiplying the specific oxidation rate at any station by the measured ambient methane concentration (T = 0). The areal methane oxidation rate (µmol m^–2 h^–1) was calculated by multiplying the average water depth to the volumetric methane oxidation rate. The precision of methane oxidation rate estimates was evaluated using the standard error (SE) of the regression slope obtained from concentration vs. time plots. The average SE was 30.1 ± 33.1%. A tendency was observed for the SE to increase as the slope decreased, indicating reduced measurement precision at lower oxidation rates.

The percentage of saturation of methane in the water samples was calculated as:

(2) $$Methanesaturation={\displaystyle \frac{C_w}{C_{eq}}\times 100}$$where C_w is the measured methane concentration in water and C_eq (nmol L^–1) is the corresponding equilibrium methane concentration in river water that is in equilibrium with the ambient atmosphere at the in situ pressure and temperature (Wiesenburg & Guinasso, 1979), assuming atmospheric CH₄ concentration was 1.9 ppm.

The diffusive flux of methane (F_outgassing, mmol m^–2 d^–1) from the river surface to the atmosphere was calculated as:

(3) $$F_{outgassing}=k\times \left(C_w-C_{eq}\right)$$where k is the integrated gas transfer coefficient (m s^–1) for methane that incorporates physical processes. The Schmidt number (Sc) for CH₄ was calculated as a function of in situ water temperature (T, °C), following (Wanninkhof, 1992):

(4) $$Sc=1\text{,}897.8-114.28\times T+3.2902\times T^2-0.039061\times T^3.$$

To ensure robustness, five empirical models were used to estimate the gas transfer velocity (k) and its standardized from k₆₀₀ (i.e., normalized to Sc number of 600 at 20 °C).

The coefficient k was calculated as (Clough et al., 2007):

(5) $$k=\sqrt{{\displaystyle \frac{DV}{h}}}+2.78e^{-6}{u}_{10}^2{\left({\displaystyle \frac{S_c}{600}}\right)}^{0.5}$$where $\sqrt{{\displaystyle \frac{DV}{h}}}$ is the water current term, which was calculated using the river water velocity (V; m s^–1), average river depth (h; m), and a diffusion coefficient for methane in the water (D; m² s^–1) (Jähne, Heinz & Dietrich, 1987). $2.78e^{-6}{u}_{10}^2{\left({\displaystyle \frac{S_c}{600}}\right)}^{0.5}$ is a wind term. 2.78e^–6 is a conversion factor (cm h^–1 to m s^–1), α is a constant (0.31), u₁₀ is the wind speed at a height of 10 m above the river, and Sc is the Schmidt number for methane (Wanninkhof, 1992).

The wind speed data (u10) was obtained from nearby weather stations located approximately 9–13 km from the sampling sites, which are available from the Automated Meteorological Data Acquisition System (AMeDAS), operated by the Japan Meteorological Agency (JMA) (https://www.jma.go.jp/jp/amedas/). The k₆₀₀ were calculated using the following empirical models (Alin et al., 2011; Raymond et al., 2012):

(6) $$-\mathrm{R}\mathrm{a}\mathrm{y}01:k_{600}=5\text{,}037\times (V\times S{)}^{0.89}\times h^{0.54}$$

(7) $$-\mathrm{R}\mathrm{a}\mathrm{y}02:k_{600}=5\text{,}937\times (1-2.54\times F{r}^2)\times (V\times S{)}^{0.89}\times h^{0.58}$$

(8) $$-\mathrm{R}\mathrm{a}\mathrm{y}05:k_{600}=2\text{,}841\times V\times S+2.02$$

(9) $$-\mathrm{A}\mathrm{l}\mathrm{i}\mathrm{n}03:k_{600}=3.84\times {10}^{-5}+9.72\times {10}^{-5}\times V$$where S is the average slope (unitless) between bridges derived from DEM data obtained from Geospatial Information Authority of Japan (https://www.gsi.go.jp/). The slopes between Nagamine and Shinfukurai bridges, and Joso and Yamato bridges were 0.000915 and 0.001395, respectively. Fr is the Froude number (Fr = V/(gh)^0.5). These models were chosen because they were preferably recommended for streams and small rivers (Raymond et al., 2012). The final gas transfer velocity k (cm h^–1) was calculated as:

(10) $$k=k_{600}\times {\left({\displaystyle \frac{S_c}{600}}\right)}^{-0.5}.$$

The resulting k values from each model were: Clough: 8.74 ± 3.39, Ray01: 23.0 ± 21.8, Ray02: 26.3 ± 24.1, Ray05: 12.7 ± 4.6, Alin03: 24.6 ± 9.5 (mean ± SD; units: cm h^–1). The diffusive methane flux was calculated separately using each model. The standard deviation of the five flux estimates was then used as an indicator of model-based uncertainty. This approach provides a quantitative assessment of variability among widely used gas transfer models, as recommended for small river systems.

The obtained methane concentrations, oxidation rates and outgassing fluxes were combined into a simple mass-balance model to estimate the internal methane loading (Fig. 1). Here, concentration changes of methane between upstream and downstream (Nagamine~Shinfukurai bridges and Joso~Yamato bridges) were modeled based on: (a) the external loading of methane via such as tributary and groundwater inflow (F_external; mol m^–3 d^–1), (b) the internal methane loading mainly produced in river bed (F_internal; mol m^–2 d^–1), (c) methane oxidation in water column (F_oxidation; mol m^–3 d^–1), and (d) methane emission to atmosphere (F_outgassing; mol m^–2 d^–1) as:

(11) $$\begin{array}{rl}{C}_{downstream}& -C_{upstream}={\int }_{t1}^{t2}{F}_{external}dt+{\displaystyle \frac{A}{V}{\int }_{t1}^{t2}{F}_{internal}dt-{\int }_{t1}^{t2}{F}_{oxidation}dt}\\ & -{\displaystyle \frac{A}{V}{\int }_{t1}^{t2}{F}_{outgassing}dt}\end{array}$$where C is dissolved methane concentration (nmol L^–1), t is streaming time (h) calculated by streaming distance (I; m) divided by velocity (u; m s^–1), A is riverbed area (m²) and V is water volume between the bridges (upstream and downstream) (Fig. 1). The F_internal refers to the net rate of methane supply to the river water, including methane production and oxidation processes in the riverbed. If the reach is sufficiently short such that the discharge at the upstream and downstream ends are nearly identical, F_external can be assumed to be zero. This assumption allows F_internal to be estimated using a mass balance approach. Since the flow velocity at Shinfukurai Bridge on August 23, 2023, could not be measured due to the wind at the time of the observation, the calculation assumed that the flow velocity at Shinfukurai Bridge was the same as that of the upstream Nagamine Bridge. In the present study, F_internal was calculated assuming that CH₄ was uniformly loaded from the riverbed within the section. The uncertainty of the internal CH₄ loading estimate (F_internal) was calculated by propagating the independent errors associated with CH₄ oxidation rates (F_oxidation) and diffusive fluxes to the atmosphere (F_outgassing). Specifically, the standard error of F_internal was derived using the root-sum-square method, assuming the errors in F_oxidation and F_outgassing were uncorrelated. To ensure the validity of F_internal estimation, data were excluded when the relative difference in discharge between upstream and downstream sites exceeded 30%, as this may indicate substantial external inflows or exceed expected measurement uncertainty. For all other cases, we assumed that the discharge remained approximately constant within the reach, considering typical measurement uncertainty.

Statistical analysis

Partial least squares (PLS) regression was used to identify significant drivers (independent X variables) of internal methane loading (F_internal, dependent Y variable). PLS was chosen because of its robustness against multicollinearity among X variables and its insensitivity to deviations from normality (Wold, Sjöström & Eriksson, 2001). The X variables included in the analysis were water velocity, average depth, discharge, water temperature, DO, PO₄³⁻, NH₄⁺, NO₃⁻, DOC, chlorophyll a, and BP. Before running the PLS regression, data were mean-centered and scaled to unit variance to mitigate the effects of differing measurement units and ensure numerical stability. To assess the importance of the independent X variables in the model, the variable influence on projection (VIP) score was used. Independent X variables with VIP > 1.2 were considered significant. PLS regression was conducted using JMP 14.3.0 (SAS Institute Inc., Cary, NC, USA).

Environmental variables

The average water depth in the upstream reach (Nagamine and Shinfukurai bridges) ranged from 2.8 to 3.5 m between April and August, exceeding the depth in the downstream reach (1.1 to 2.1 m) during the same period (Fig. 3A). This difference is attributed to the operation of the Fukuoka weir for irrigation, which increased water depth upstream. Similarly, flow velocity was slower in the upstream reach during the irrigation period (April to August, Fig. 3B). On September 26, 2022, both water depth and flow velocity were unusually high due to flooding. The discharge for each month was almost the same between the two bridges in each upstream and downstream reach (Fig. 3C). In certain months (June 2022 in upper reach and September 2022 and May 2023 in lower reach), the discharge difference between upstream and downstream sites exceeded the 30% threshold. For these months and reaches, F_internal values were not calculated. In all other months, discharge differences were within the assumed uncertainty range, and F_internal was estimated accordingly.

Water temperature varied from 1.1 °C (January) to 29.1 °C (August; Fig. 4A). Dissolved oxygen (DO) concentration was lower in the upstream reach during the irrigation period, while NH₄⁺ and PO₄^3– concentrations were higher upstream during the same period (Figs. 4B and S2A, S2B). Specific conductivity was higher during the non-irrigation period (maximum: 231 µS cm⁻¹) and lower during irrigation (minimum: 163 µS cm^–1) (Fig. S2A). Seasonal trends in nitrate (NO₃⁻) concentrations showed higher values in winter and lower values in summer (Fig. 4C). Dissolved organic carbon (DOC) concentrations ranged from 1.2 to 2.3 mgC L^–1 (Fig. S3A). Bacterial production (BP) was lowest in winter and increased in summer, with no clear differences between sites (Fig. S3B). Chlorophyll-a (Chl) concentrations were generally higher in the downstream reach during irrigation (Fig. S3C). The complete dataset, including all raw measurements of environmental parameters and methane-related variables, is provided in File S1.

Methane dynamics

Dissolved methane concentrations fluctuated between 237 ± 19 and 1,271 ± 6 nM (Fig. 5A). Within a given month, methane concentrations tended to decrease downstream. However, higher concentrations were sometimes observed during the irrigation period at Shin-fukurai Bridge compared to Nagamine Bridge in the upper reach. In the lower reach, methane concentrations at Joso Bridge were consistently higher than at Yamato Bridge in both irrigation and non-irrigation periods. The specific methane oxidation rate ranged from –0.00165 ± 0.00046 to 0.0328 ± 0.00105 h^–1, with relatively high values observed in downstream reaches such as Joso and Yamato during the summer irrigation period. Areal and volumetric methane oxidation rates peaked at 66.4 ± 18.5 µmol m^–2 h^–1 and 23.2 ± 6.5 nmol L^–1 h^–1, respectively (Fig. 5B). Notably, during the irrigation period, areal methane oxidation rates were higher in the upstream section (Nagamine and Shin-fukurai: 16.2 ± 18.1 µmol m^–2 h^–1) than in the downstream section (Joso and Yamato: 5.97 ± 4.68 µmol m^–2 h^–1), likely due to greater water depth in the upstream reaches. In contrast, from October to march, the areal oxidation rates showed little spatial variation, with an overall average of 1.1 ± 1.3 µmol m^–2 h^–1. Methane diffusive emission fluxes to the atmosphere ranged from 199 ± 149 to 31.6 ± 9.6 µmol m^–2 h^–1 (Fig. 5C). During the irrigation season, diffusive fluxes increased in both upstream (Nagamine and Shin-fukurai: 88.9 ± 26.7µmol m^–2 h^–1) and downstream reaches (Joso and Yamato: 107 ± 44 µmol m^–2 h^–1). From October to March, the values remained low throughout the river (mean: 53.0 ± 12.7 µmol m^–2 h^–1).

Internal methane loading (F_internal) varied from –33.1 ± 20.7 to 160 ± 46 µmol m^–2 h^–1, with a mean of 61.3 ± 50.8 µmol m^–2 h^–1 (Fig. 6). In both the upper and lower reaches, F_internal tended to increase during the summer and decrease during the winter, indicating a consistent seasonal pattern. No significant difference in F_internal was found between the upper reach (66.0 ± 51.6 µmol m^–2 h^–1) and the lower reach (56.2 ± 51.4 µmol m^–2 h^–1) throughout the year (Welch’s t-test, n = 27, p > 0.05). However, during the latter half of the irrigation period (July to August), F_internal was significantly higher in the upper reach (149 ± 11 µmol m^–2 h^–1)than in the lower reach (92.4 ± 23.8 µmol m^–2 h^–1) (Welch’s t-test, n = 7, p = 0.0112).

In the PLS regression analysis, 58.2% of the variation in F_internal was explained by the explanatory variables (R² = 0.582, n = 27, p < 0.001). Average water depth, water temperature, DO and NO₃^– were selected as explanatory variables with VIP > 1.2 for the dependent variable of F_internal (Table 1). The estimated coefficients of DO and NO₃^– were negative, whereas those of average water depth and water temperature were positive.