Here we use biogeochemical modules of the Bottom RedOx Model BROM (Yakushev et al., 2017) coupled with a vertical 2 Dimensional Benthic-Pelagic Model 2DBP ( Figure 1 ), that considers processes occurring in the water column, benthic boundary layer (BBL), and sediments together (Yakushev et al., 2020).

The BROM-C module for carbon transformation describes processes in the carbonate system allowing calculation of pH and carbonates saturation states, and rates of formation and dissolution of carbonates ( Figure 2 ). It also parameterizes methanogenesis and both aerobic and anaerobic methane oxidation. The detailed description is given in (Yakushev et al., 2017). To describe the background biogeochemical processes of seasonal formation and decay of organic matter we used the biogeochemical module OxyDep considering transformation of nitrogen between inorganic nitrogen (NUT), phytoplankton (PHY), heterotrophs (HET), dissolved organic matter (DOM), particulate organic matter (POM) and changes in dissolved oxygen content (OXY) (Berezina et al., 2021) ( Figure 2 ). The gas bubble fate module Yakushev et al. (2021) parameterizes bubbles rising and dissolution (see Supporting Material).

These modules were integrated into an existing modular platform [Framework for Aquatic Biogeochemical Modelling, FABM; (Bruggeman and Bolding, 2014)], that couples together the transport driver 2DBP, BROM’s biogeochemical blocks, and the bubble fate module.

The superiority of this model lies in its ability to simulate the transport of CH₄ in both gaseous and dissolved forms, its dissolution and oxidation, as well as the transformation of carbon-containing compounds through chemical processes (e.g., carbonate system equilibration) and biogeochemical processes (e.g., those related to the synthesis and decay of organic matter).

In our studies, the transport model 2DBP represents a 20 km transect with a seep at the center. The transect is positioned in the central part of the Laptev Sea in a region with seep fields characterized by massive CH₄ ebullition [Polygon P1, from (Shakhova et al., 2015)]. 2DBP was forced by model outputs for the seasonal variability of temperature, salinity, ice cover changes, vertical turbulence, irradiance and current velocity from a ROMS-20 km model (Shchepetkin and McWilliams, 2005) applied for the Arctic Ocean (Ommundsen et al., 2008).

The horizontal spatial resolution of 2DBP was 25 m in the center of the transect and increased toward the peripherical parts to 500 m. Vertical resolution of 2DBP was the same as of ROMS (40 grid points for 50 m depth) in the water column and decreased in the limits of 50 cm BBL from 13 to 6 cm and then increased in the upper sediments from 0.5 mm to 2 cm.

After 200 spin-up years with repeated seasonal variations in temperature, salinity, ice cover, vertical turbulence, irradiance, and current velocity, a stable seasonal pattern of biogeochemical characteristics was achieved. This was used as an initial condition for the seeping experiment, with seepage occurring in a 25 m x 25 m area in the center of the transect ( Figure 3 ). Methane release in bubbles was parameterized at 50 cm above the sediment at the upper boundary of the BBL. The flux of methane through the sediment in the release area was not considered, because in reality the bubbles are released through narrow channels but are not “filtered” through the sediment (Sapart et al., 2017; Shakhova et al., 2015; Wild et al., 2022, 2023).

Existing estimates of methane release are in the range 30–170 g m⁻²d⁻¹ or 2500–15000 mM m⁻²d⁻¹ (Chernykh et al., 2023; Shakhova et al., 2019, 2015). Here we made numerical experiments for a moderate seeping (MS) scenario of 3600 mM m⁻²d⁻¹ CH₄ and an intensive seeping scenario of 20000 mM m⁻²d⁻¹ CH₄, consistent with predicted (Shakhova et al., 2010a) and observed flux estimates (see also references in Table 1 ). We assumed that the seepage produces a constant number of bubbles with the same initial diameter, which are transported and dissolved according to the single-bubble dissolution model described in the Supplementary Material .

In the model we parameterize aerobic methane oxidation, a microbially-mediated process, that is generally described with an equation (Chan et al., 2019):

It is widely assumed that the oxidation of CH₄ in the water column follows first-order kinetics (Kitidis et al., 2010; Wåhlström and Meier, 2014), and there exist estimates of the rates of aerobic methane oxidation in the sea water determined using different techniques for different concentrations of methane and for different oceanic regions. Most of the estimated CH₄ oxidation rates fall into the range 0.001–10 nM d⁻¹, and CH₄ oxidation activity increases in ocean environments with high CH₄ concentrations (Mau et al., 2013). The CH₄ oxidation fluxes for the outer ESAS regions measured with a C₃H₄ radiotracer following the procedure of (Sandbeck and Reeburgh, 1989) were low, about 0.1 ± 0.13 nM d⁻¹, corresponding to rates 0.001–0.003 d⁻¹ (Shakhova et al., 2015). Higher rate estimates were obtained by (Mau et al., 2013) for Arctic water in a Spitzbergen fjord: 0.01–0.08 d⁻¹ (using ³H-CH₄) and 0.001–0.1 d⁻¹ (using ¹⁴C-CH₄). Bussmann et al. (2017) found that CH₄ oxidation rates varied between riverine (0.011 d⁻¹), mixed (0.006 d⁻¹), and polar (0.028 d⁻¹) waters of the Laptev Sea. Incubation studies in the deep Gulf of Mexico obtained sometimes higher CH₄ oxidation rates: 0.2 ± 0.1 d⁻¹ (Chan et al., 2019), 0.0001–0.200 d⁻¹ (Kessler et al., 2011), 0.0001–0.425 d⁻¹ (Crespo-Medina et al., 2014; Weinstein et al., 2016) estimated rates 0.012–0.13 d⁻¹ at different depths of the Hudson Canyon.

In this model, the specific oxidation rate of CH₄ by O₂ was set to k_{CH 4} =0.01 d⁻¹, which lies in the midrange of estimated oxidation rates. Fluxes through the sea surface for the gases (oxygen, carbon dioxide, methane) were parametrized as bubble-mediated air-sea gas exchange (see Supplementary Material 2 ). The influence of ice was considered as a restriction on the flux scaling with its fractional coverage of the sea surface.

The simulated seasonal variability after 200 years spin-up reflects the basic features of biogeochemical seasonality in the water column in this region ( Figure 3 ). An intensive phytoplankton bloom starts after the ice melts in August followed by zooplankton growth and increasing concentrations of particulate and dissolved organic matter (POM, DOM). Maximum concentrations of phytoplankton (Phy) reach 0.7 µM N at depth 25 m; this suggests a maximum Chl a concentration of ca. 2 mg/m³ (assuming a Redfield molar ratio of C:N = 106/16 and a typical low-light-adapted C:Chl ratio of 30 gC/gChla; (Cloern et al., 1995), which is consistent with observed values of 1–2 mg/m³ for the Laptev Sea in September 2008 (Polyakova et al., 2021). Zooplankton (Het) increases to ∼0.2 µM N, roughly uniform over the 50 m column depth; this corresponds to ∼2 g dw/m² [assuming Redfield C:N and carbon 42% of dry weight (Beers, 1966)] which is consistent with values 0.1-1.5 g dw/m² observed on the Laptev Sea shelf during autumn 1993 (Kosobokova et al., 1997). Concentrations of POM and DOM increase after the bloom to 0.1 µM N and 3.5 µM N respectively. Concentrations of NUT (i.e. nitrate and nitrite) are exhausted during the phytoplankton bloom and restored to maximum concentrations of ∼7 µM N in winter. The water column is well oxygenated. Dissolved oxygen has highest concentrations in winter (about 380 µM) and lowest in the bottom layer in summer after the bloom, when oxygen is consumed by decomposition of organic matter (310 µM). These features generally corresponds to the observed changed of the distributions in the Laptev Sea in the ice free period from August to November analyzed in (Xie et al., 2023). Simulated baseline changes of pH (total scale) show reasonable seasonal variability for pH, it varies from 7.9 in January to 7.4 in August. pH decreases during the summer period due to influx of low salinity and low alkalinity riverine water, and increases weakly during the phytoplankton bloom. Aragonite saturation state varies from 2.05 in January to 0.8 in August. DIC is around 2300 µM C during spring and summer and starts to decrease at the end of August to 2100 µM in October, due to phytoplankton uptake, then gradually increases to 2200 µM at the end of November, due to DOM decomposition ( Figure 3 ). Generally, the surface water of this region is strongly affected by river discharge; the surface layer with high salinity and alkalinity is replaced by a low alkalinity plume water during the warm season, and alkalinity increases again in November ( Figure 3 ). Note however that the seepage area modeled in this study lies outside the area of extreme acidification impacted by export of eroded carbon and very strong river influence (Semiletov et al., 2016).

To summarize the statistical characteristics for the intercomparison between the model and the ensemble means values of the observed data collected in 2015–2020 (Xie et al., 2023), we use the cost function based method as described in (Eilola et al., 2011). The cost function C is computed for the model as follows:

Here, the mean bias (M − D) of the model means (M) relative to the mean of observations (D) is normalized by the standard deviation (SD) of the observations. The cost function for the bias of the ensemble mean is calculated similarly. The cost function value is determined separately for each variable (salinity, temperature, oxygen, nutrients, pH and alkalinity) at each station and depth. The following ranges are used for interpreting C: 0 ≤ C < 1 (good), 1 ≤ C < 2 (reasonable), and 2 ≤ C (poor).

Overall, nutrients, alkalinity, and pH are modeled well with good accuracy (NUT = 0.69, Alk = 0.92, pH = 0.51). Oxygen and salinity have reasonable accuracy (Oxy = 1.16, S = 1.86), suggesting some room for improvement in modeling these variables. The poor accuracy of the temperature model (T = 2.68) can be potentially linked to the observations being taken during anomalously warm years 2019 and 2020, in the Laptev Sea. These years were characterized by higher than average temperatures, which might not have been adequately represented in the historical data used to develop the model.

Vertical distributions of the modeled and observed characteristics are shown in Supplementary Figure S2 in Supplementary Material .

Starting from the model’s date 01.01.2012, bubble methane ebullition was parameterized in the center of the transect with a constant intensity of the mentioned above either 3600 mM CH₄ m⁻²d⁻¹ [moderate seeping (MS) scenario] or 20000 mM CH₄ m⁻²d⁻¹ [intensive seeping (IS) scenario].

Depending on the initial methane flux and bubble radii, bubble dissolution takes place at different depths with different intensity. Calculated vertical distributions above the seep of methane in bubbles, bubble radii, rising rate, its vertical flux and dissolved methane are shown in Table 1 . We see that the shapes of these distributions are sensitive to initial bubble diameter. As bubbles ascend from the bottom, their radii decrease due to dissolution but increase in response to the reduced pressure in the upper layers. The model shows that larger bubbles rise faster and dissolve more slowly, causing their radii to grow as they ascend ( Figure 4 ). Consequently, smaller bubbles can be completely dissolved within the water column ( Table 1 ). Such an effect was also observed in Arctic thermokarst lakes (Semiletov et al., 1996). Note that the modeled dissolved methane concentration is higher in the bottom water that is constantly subjected to intensive bubble dissolution ( Table 1 , right column). These concentrations are higher than published observations (see references in Table 1 ), but lower than recent estimates obtained for an intensive seep in the Laptev Sea (Semiletov, 2024, pers. comm.). The methane oxidation rate used here, k_{CH 4} =0.01 ⁻¹, is closer to the minimum estimates (see above) and its increase leads to a decrease of the bottom concentration of dissolved methane.

The modeled ebullition starting 01.01.2012 leads to increased dissolved methane concentrations and changes in the carbonate system ( Figures 4 , 5 ). Aerobic oxidation of the dissolved methane causes production of CO₂, a corresponding increase of DIC, and decrease in pH and dissolved oxygen ( Figures 4 , 5 ).

Figure 4 shows small-timescale changes of CH₄, DIC, pH, OXY and Ω _Ar above the seeping point, that can be explained by lateral currents parameterized in the model (mainly of tidal scale), that shift water back and forth water in horizontal direction. The processes of methane oxidation take a certain time, so the anomalies in the fields of oxygen, DIC and pH can be displaced horizontally from the point with maximum bubble concentration ( Figure 6 ) resulting in an export of “methane-oxidation effect” to the adjacent waters.

The region in the vicinity of the seep is in a larger degree affected by seeping and demonstrates long-term changes in the carbonate system. Figure 4 shows that these changes take several years (about 7 years in this MS model assumptions, i.e. intensity of seeping, considered volume of water, bottom depth etc.).

CO₂ produced from the oxidizing methane contributes to an increase in DIC, that can diffuse into the sediments in the region near the seep. This increase of DIC (with no change in alkalinity) results in a decrease of aragonite saturation in the bottom boundary layer and in the upper sediments. Therefore, the model predicts dissolution of calcium carbonate (aragonite form) in the upper layer of the sediments ( Figures 5 , 7 ). Potential dissolution of calcium carbonate present in the sediments buffers changes in alkalinity in the pore water and therefore in the bottom water that keeps the water column seasonal variability features the same after the start of seeping.

In the bottom layer, the most intensive changes should be detected at about 1 km distance from the seep ( Figures 6 , 7 ), but in the water column the carbonate system parameters can be affected up to 10 km from the seep ( Figure 7 ) with the MS model assumptions.