Modeling water column gas transformation, migration and atmospheric flux from seafloor seepage

Understanding the fate of gas seeping from the seafloor is crucial for assessing the environmental impacts of both natural and anthropogenic seep systems, such as CH₄ cold seeps, leaking gas wells, and future carbon capture projects. We present a comprehensive modeling framework that integrates physical, chemical, and biological processes to estimate the 3-dimensional water column dissolved gas concentration field and 2-dimensional atmospheric flux field resulting from seafloor seeps. The framework consists of two main components: (1) a gas-phase model that calculates free gas dissolution and direct atmospheric release at the seep site, and (2) a concentration model that combines particle dispersion modeling with an adaptive-bandwidth kernel density estimator and customizable process modules. Applying the framework to a natural CH₄ seep at 200 m depth offshore northwestern Norway (20 May–20 June 2018), we found that dissolved methane was advected northeastward along the coast, spreading across shelves, reefs, and into fjord systems. Within days, the vertical CH₄ concentration profile was near inverted, with near-surface maxima, facilitating atmospheric exchange. Diffusive emissions covered large areas (>10⁵ km²) and was almost 3 times the local free gas flux. Around 0.7 % of dissolved CH₄ reached the atmosphere during a 4 week period, microbial oxidation removed around 65 %, while ∼ 34 % remained in the water column. Uncertainties caused by a range of model framework elements remain substantial, e.g. can estimates of microbial oxidation removal change from 65 % to as low as 5.5 % or as high as 91.4 % depending on rate coefficient assumptions. Our framework provides a globally applicable tool that integrates free and dissolved gas dynamics and accommodates advanced hydrodynamic modeling. Its ability to explicitly resolve spatiotemporal fields enables the inclusion of complex physical and biogeochemical process modules and supports not only the quantification of atmospheric fluxes but also applications that require explicit field representations, such as assessing impacts on local ecosystems.

1 Introduction

Estimates of the contribution of seafloor gas seepage to atmospheric emissions and its impact on ocean environments are highly uncertain due to limited data and understanding of gas transformation and transport mechanisms in the water column. Estimation of total atmospheric gas emissions from seep areas (e.g. Myhre et al., 2016) rely largely on either ship measurements or large-scale atmospheric inversion models. The former of these approaches only gives information on the local flux and requires some sort of up-scaling, while the latter is unable to estimate dispersed sources and/or weaker point sources precisely due to its rough scale and inability to completely decouple atmospheric sources from sinks (Thompson and Stohl, 2014). Quantifying dissolved gas in the water column usually involves measuring dissolved gas via water samples (e.g. Silyakova et al., 2020) or using in situ sensors (e.g. Gentz et al., 2014) which can be time-consuming and often result in poor data coverage. New modeling tools for constraining the environmental impacts of current and future seabed gas seepage from both natural and man made sources are therefore needed.

Gas released at the seabed can enter the atmosphere directly as free gas (bubbles) or via diffusive equilibrium of dissolved gas that has reached the sea surface. To estimate the total atmospheric emissions from a seabed seep and its dissolved distribution in the water column, one must be able to model both pathways simultaneously. Gas content in bubbles is constantly changing due to dissolution (gases in the bubble dissolve in the liquid) and exsolution (gases already dissolved in the liquid enter the bubble) driven by partial pressure gradients across the bubble rim. Additionally, chemical and biological processes can modify local dissolved gas content. Estimating the gas distribution in the water column and total atmospheric flux therefore requires a flexible framework which can integrate processes governing the gas phase dynamics and the hydrodynamics, accommodate atmospheric exchange, and other phenomena that modify water column gas content. Previous modelling efforts have typically focused on single gas phase frameworks including only selected processes (e.g. McGinnis et al., 2006; Graves et al., 2015; Silyakova et al., 2020), however, key steps towards modeling the complete system have been made recently in Dissanayake et al. (2023) and Nordam et al. (2025). We aim to further expand on these studies from a methodological perspective and provide a pilot framework which can integrate all key processes governing free and dissolved transport and transformation of seeped gas and give a realistic estimate of the time varying 3-dimensional (3D) water column concentration field and 2-dimensional (2D) atmospheric release field.

Our approach integrates a gas phase model with a hydrodynamic model using particle dispersion modeling (similar to Dissanayake et al., 2023). It estimates the 3D distribution of gas in the water column and the total (free and diffusive) atmospheric 2D gas release resulting from observed seabed seepage. This approach offer flexible inclusion of atmospheric flux and chemical and biological process modules affecting dissolved gas content in the water column. Explicit concentrations (molecules per volume) are obtained using kernel density estimation. Atmospheric dissolved flux estimates are obtained using a bulk model and atmospheric free gas flux via a gas phase model. We tested the framework by quantifying direct and diffusive atmospheric fluxes as well as 3D dissolved gas distribution between 20 May and 20 June 2018 for a methane (CH₄) seep area offshore Northwestern Norway.

2 Method

Our goals are two-fold: (i) Calculate the combined total amount of seep-derived gas that reaches the atmosphere – both direct free gas release and ventilation of dissolved gas, and (ii) Estimate the impact of seeped gas on the scalar dissolved gas concentration field, i.e., we seek the anomaly ${\mathrm{\Phi }}^{\prime }(x,y,z,t)=\mathrm{\Phi }(x,y,z,t)-{\mathrm{\Phi }}_{\mathrm{0}}(x,y,z,t)$ caused by the seeps, where $\mathrm{\Phi }(x,y,z,t)$ denotes the total concentration and ${\mathrm{\Phi }}_{\mathrm{0}}(x,y,z,t)$ is a background concentration.

The outlined goals are achieved by adapting and integrating existing and new models, of which output from one model serves as a final result or feeds another model. We first use seabed gas volume flux data and a two-phase gas model to calculate the gas dissolution rates and direct atmospheric gas (bubble) release. Dissolved injection rate output from the gas phase model then feeds a concentration model that combines an existing dispersion modeling framework with an adaptive kernel density estimator, including an atmospheric flux module and options for water column process modules. Figure 1 shows the complete framework, with input data in the left column, the modeling steps in the center column, and the final results in the right column. Each modeling section is detailed in the corresponding subsection.

Figure 1 Model framework flowchart. The emboldened modeling steps and associated numbers refer to the four subsections of the Methods section. The “D” in the results column refers to spatial dimensions.

2.1 Gas phase modeling

Free atmospheric gas fluxes and dissolved gas profiles injected to the water column are initially modeled for each observed seep using the seabed free gas flux data and the M2PG1 gas phase model (Jansson et al., 2019). M2PG1 provides an integrated solution of dissolved and free gas in a 1-dimensional water column, with sources and sinks at both horizontal and vertical model boundaries. It simultaneously models gas exchange, dissolution, and associated dissolved gas concentration of five gas species (methane (CH₄), Argon (Ar), Carbon dioxide (CO₂), Nitrogen (N₂), and Oxygen (O₂)) across a user-defined initial spectrum of bubble sizes. The bubble size spectrum and gas distribution across this spectrum vary freely across the spatio-temporal model domain. The model includes several bubble shape and rising speed models, microbial oxidation of CH₄ using first order kinetics (Griffiths et al., 1982; Chan et al., 2019), diffusive exchange with the atmosphere, dissolved transport due to vertical turbulent exchange of water masses, as well as loss due to advection across the model boundary (Jansson et al., 2019). While dynamic solutions are permissible in M2PG1, we have opted for a steady state solution in our modeling framework.

The input parameters include seabed gas flux, bubble characteristics (size distribution, rising speed, dirtiness, flatness), temperature, salinity, microbial CH₄ oxidation rate coefficients (MOx), ambient dissolved gas concentrations (for all five gases), vertical mixing (turbulent) and local ocean currents. Seabed free gas flux data can in theory be obtained by any means available, although hydroacoustics have been used extensively due to its relatively straightforward deployability and large coverage (e.g. Ferré et al., 2020).

In our implementation of M2PG1 we used a new estimation technique to determine the horizontal model domain size in M2PG1. Horizontal domain size was previously chosen ambiguously in M2PG1 (Jansson et al., 2019) and could cause significant exchange rate errors. Our method removes this ambiguity by estimating the horizontal bubble plume extent based on local conditions. Details are provided in Appendix A.

The steady-state output from the M2PG1 simulation provides two key results: (i) direct atmospheric gas flux and (ii) injection rates of dissolved gas to the surrounding water column. The former is a direct output in M2PG1 and the latter, which are key input for the concentration modeling steps (Sect. 2.2–2.4), can be derived from the dissolved gas concentration profiles by calculating the dissolved gas loss q [mol s⁻¹] to the water column at the downstream boundary of each M2PG1 grid cell. The steady-state mass flux assumption gives:

$$\begin{array}{}\text{(1)}& q=A_{\mathcal{M}}^{⟂}U({\mathit{\phi }}_{\mathcal{M}}-{\mathit{\phi }}_{\mathrm{b}}),\end{array}$$

where $A_{\mathcal{M}}^{⟂}$ [m²] is the vertical grid cell area (Appendix A), U [m s⁻¹] the current speed, φ_ℳ [mol m⁻³] the estimated concentration within the grid cell and φ_b [mol m⁻³] the assumed concentration at the upstream boundary.

2.2 Particle dispersion modeling

To estimate the unobservable dissolved gas concentration field anomaly ${\mathrm{\Phi }}^{\prime }(x,y,z,t)$, we must model the advection and spread of the dissolved gas from the seeps. We chose to simulate the transport and dispersion of the gas from the release site using OpenDrift, which is a Lagrangian particle trajectory modeling software (Dagestad et al., 2018). In practice, this means that we distribute (virtually) the released CH₄ over a discrete number of virtual particles, and update the particle positions at discrete times t_n for time-steps $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },N$ according to the output from a hydrodynamic model. Each timestep is separated by a time interval Δt. We then define S[n] as the number of virtual particles seeded at the modeled seep sites at each time step. This generates a total of $Z={\sum }_{n=\mathrm{1}}^{N}S\left[n\right]$ particles indexed by $\mathit{\zeta }=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },Z$. Note that we throughout this manuscript will use square “[⋅]” versus round “(⋅)” brackets to distinguish between discrete and continuous spatiotemporal arguments, respectively.

Once particles are seeded, OpenDrift calculates the trajectory of each particle individually by numerically solving a stochastic differential equation which is consistent with the Lagrangian representation of the advection-diffusion equation (see e.g. Spivakovskaya et al., 2007). The drift in particle position η can be expressed as

$$\begin{array}{}\text{(2)}& \mathrm{d}\mathit{\eta }={\mathit{U}}_{\mathit{\mu }}(\mathit{\eta },t)\mathrm{d}t+\mathbf{B}(\mathit{\eta },t)\mathrm{d}\mathit{W}\left(t\right)\end{array}$$

where U_μ(η,t) represents displacement produced by the underlying (mean) velocity field and the second term represents displacement from random, diffusive processes and is composed of a diffusivity matrix B(η,t) and increments of a Wiener process dW(t). The advective term (U_μ(η,t)) is determined by velocity fields obtained from the hydrodynamic model. OpenDrift represents the diffusivity B(η,t) as a diagonal matrix with a horizontal and a vertical diffusivity. If available, these diffusivities can be directly read from the hydrodynamic model output. Otherwise, OpenDrift can also estimate the diffusivity coefficients using one of several built-in parametrizations. Finally, OpenDrift returns individual (traceable) positions η_ζ[n] for each seeded particle at each time-step they spend in the model domain.

2.2.1 Particle mass

To associate the particle distribution with dissolved gas content, we latch a particle mass Γ_ζ to each seeded particle, which explicitly corresponds to the number of moles each particle represents (this mass has no influence on the particle buoyancy). Each particle is thus interpreted as a virtual single-point representation of some local spatial distribution of Γ_ζ moles of dissolved gas molecules.

The initial mass, i.e. mass at release, of an arbitrary particle ζ is scaled such that the total released particle mass from all modeled seeps combined at timestep n approximates the total number of moles of gas dissolved in the water column during the time interval Δt centered on t_n. In practice, we distribute the integrated sum of modeled (using the gas phase model) injected gas molecules from $t_n-\mathrm{\Delta }t/\mathrm{2}$ to $t_n+\mathrm{\Delta }t/\mathrm{2}$ evenly over the seeded particles. The mass of particle Γ_ζ seeded at time-step n is then obtained by

$$\begin{array}{}\text{(3)}& {\mathrm{\Gamma }}_{\mathit{\zeta }}\left[n\right]={\displaystyle \frac{\mathrm{\Delta }t\sum _{p=\mathrm{0}}^P{\mathrm{{\rm Y}}}_p\left[n\right]}{S\left[n\right]}}\end{array}$$

where ${\mathrm{{\rm Y}}}_{\mathrm{1}}\left[n\right],{\mathrm{{\rm Y}}}_{\mathrm{2}}\left[n\right],\mathrm{\dots },{\mathrm{{\rm Y}}}_P\left[n\right]$ [mol s⁻¹] are total injected dissolved gas from all P modeled seeps. Approximation to the modeled dissolved gas release profiles at each modeled seep is achieved by seeding different amount of particles at different depths. Particle masses are then subsequently individually adjusted at each time-step to simulate processes affecting gas content. Each particle thus has a successively constructed mass time-series, where the current mass Γ_ζ[n] is determined by the previous mass Γ_ζ[n−1] and selected mass modification functions.

2.2.2 Particle count

Our framework must be able to model an extensive 3-dimensional domain (e.g. larger ocean regions), making computational complexity a challenge. Both computation time and estimation quality increase with the number of active particles present in the domain. This makes it crucial to be able to strike a decent compromise between the two, which typically involves removal of particles that have been present in the domain for a certain number of time-steps. Total particle count Λ[n] in the domain can be expressed as

$$\begin{array}{}\text{(4)}& \mathrm{\Lambda }\left[n\right]=\mathrm{\Lambda }[n-\mathrm{1}]+S\left[n\right]-L\left[n\right]-\mathrm{\wp }\left[n\right],\end{array}$$

where L[n] is the number of particles leaving the modelled geographical domain, and ℘[n] the number of removed particles. A constant particle count is obtained when $S\left[n\right]\sim L\left[n\right]+\mathrm{\wp }\left[n\right]$. Since ℘[n] represents non-physical loss of gas, the model simulation would ideally run with a spin-up time that ensures S[n]∼L[n]. Unfortunately, this typically results in unreasonable computation times and/or spin-up periods, making particle removal necessary. To limit errors caused by removed particles, we apply a function that redistributes mass from all removed particles to nearby non-removed particles. The redistribution is weighted according to the inverse distance from the removed particle within a user defined distance limit d_max, giving a non-removed particle τ an added mass of

$$\begin{array}{}\text{(5)}& {\mathit{\gamma }}_{\mathit{\tau }}={\mathrm{\Gamma }}_{\mathit{\theta }}{\displaystyle \frac{\left|\right|{\mathit{\eta }}_{\mathit{\tau }}-{\mathit{\eta }}_{\mathit{\theta }}|{|}_{\mathrm{2}}^{-\mathrm{1}}}{{\sum }_{\mathit{\zeta }\in \mathcal{T}}\left|\right|{\mathit{\eta }}_{\mathit{\zeta }}-{\mathit{\eta }}_{\mathit{\theta }}|{|}_{\mathrm{2}}^{-\mathrm{1}}}}\end{array}$$

from the removed particle θ. Here, τ∈𝒯 and 𝒯 is the set of non-removed particles with indices ζ satisfying $\Vert {\mathit{\eta }}_{\mathit{\zeta }}-{\mathit{\eta }}_{\mathit{\theta }}{\Vert }_{\mathrm{2}}\le d_{max}$, and $\Vert \cdot {\Vert }_{\mathrm{2}}$ denotes the Euclidean norm. This solution changes the problem of non-physical loss of dissolved gas to one of non-physical redistribution. This can affect model results by shifting particle mass towards the seed location, since the density of particles are in general higher closer to the release point. However, we consider this artifact less problematic than mass simply disappearing.

2.3 Grid Projected Adaptive-bandwidth Kernel Density Estimator

Having an explicit relationship between dissolved gas content (of seep origin) and particle mass Γ_ζ allows us to infer gas concentrations by evaluating the particle mass per unit volume, which we refer to as the particle density. Let us assume that the particles $\mathit{\zeta }=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },Z$ are scaled/weighted samples (using their mass Γ_ζ) from an unknown, smooth, underlying particle density field $\mathit{\varphi }(x,y,z,t)$ which approximates the seep-induced gas concentration anomaly field ${\mathrm{\Phi }}^{\prime }(x,y,z,t)$. Estimation of ${\mathrm{\Phi }}^{\prime }(x,y,z,t)$ can then be done via the estimate $\widehat{\mathit{\varphi }}$ of $\mathit{\varphi }(x,y,z,t)$ using the particle data set.

To get $\widehat{\mathit{\varphi }}$, we employ a discrete spatiotemporal grid $[i,j,k,n]$, where $i=\mathrm{1},\mathrm{2},\mathrm{\dots },I,j=\mathrm{1},\mathrm{2},\mathrm{\dots },J,k=\mathrm{1},\mathrm{2},\mathrm{\dots },K,n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },N$, and $I,J,K,N$ denote the number of grid cells in east, north, vertical and temporal dimensions, respectively. Grid cell center positions are given by $[x_i,y_j,z_k,t_n]$, with horizontal resolution Δλ (in both directions), vertical resolution Δz, and temporal resolution Δt. We then bin all mass in the temporal and vertical domains and obtain separate estimates $\widehat{\mathit{\varphi }}[i,j]$ of ϕ(x_i,y_j) for each resulting depth layer k and time-step n to form the final estimate $\widehat{\mathit{\varphi }}[i,j,k,n]$. Obtaining $\widehat{\mathit{\varphi }}[i,j,k,n]$ thus translates to solving a series of 2-dimensional density estimation problems (see e.g. Silverman, 1986).

Due to the extensive model domain and the need to obtain one estimate for every depth layer and time-step (K×N estimates), our density estimator needs to be fast and allow reliable density estimates from limited particle counts. It must also handle regions with low and high concentrations and concentration gradients as well as complex boundaries like fjords and islands. A commonly used density estimator in similar contexts is the histogram estimator, which unfortunately has several well-known limitations in these applications (the histogram estimator and its drawbacks are detailed in Appendix B). Previous studies on concentration estimation from particle dispersion model data have shown that Kernel Density Estimators (KDEs) can offer far superior information exploitation than the histogram estimator and overcome many of its drawbacks (see e.g. De Haan, 1999; Vitali et al., 2006; Björnham et al., 2015; Barbero et al., 2024; Yang et al., 2026). One remaining challenge in our specific application, however, is the lack of available KDEs tailored to coastal ocean regions that appropriately adapt to spatial density variability (adaptive bandwidth) and complex boundary geometries (bathymetry). We therefore formulated a new 2-dimensional adaptive-bandwidth KDE to provide our density estimates.

Kernel density estimation is a standard non-parametric way to estimate the density of a random variable using kernel functions (Silverman, 1986). This offers density estimates that are differentiable, grid cell size independent and generally more realistic than histogram estimators, without lower density limitations. In our case, a kernel density estimate involves placing a symmetric, smooth, and weighted kernel function at each particle position. By summing up the kernel contributions, the density field $\widehat{\mathit{\varphi }}$ at position r₀ located within the volume V can be estimated via the general kernel density estimator formula

$$\begin{array}{}\text{(6)}& \widehat{\mathit{\varphi }}\left({\mathit{r}}_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{1}}{V}}\sum _{\mathit{\zeta }=\mathrm{1}}^Z{\mathrm{\Gamma }}_{\mathit{\zeta }}{K}_h\left(\left|\right|{\mathit{r}}_{\mathrm{0}}-{\mathit{\eta }}_{\mathit{\zeta }}|{|}_{\mathrm{2}}\right).\end{array}$$

Here, $K_h\left(\mathit{\xi }\right)\equiv (\mathrm{1}/h)K(\mathit{\xi }/h)$, where K(ξ) is a non-negative, normalized, and symmetric kernel function, h a bandwidth (smoothing) parameter, and r₀ the estimate position.

It is well established that the choice of kernel shape K(ξ) is of less importance, as long as it adheres to the kernel function requirements. We define the base kernel K(ξ) as a standardized 2-dimensional Gaussian, i.e. $K\left(\mathit{\xi }\right)=\exp \left(-{\mathit{\xi }}^{\mathrm{2}}/\mathrm{2}\right)/\mathrm{2}\mathit{\pi }$.

Selecting an appropriate bandwidth h is crucial, as a poor choice can cause large errors (De Haan, 1999), particularly due to over-smoothing (Larsen et al., 2002). Several methods exist for selecting h by evaluating the statistical properties of the collected data, but they typically rely on strict assumptions on the underlying field. For heterogeneous fields, such as ours, where statistical properties vary across the domain, local adaptation of h is necessary to give realistic density estimates in both high and low particle count areas. Furthermore, the presence of complex boundaries in the form of bathymetry and coastlines introduces additional challenges, both for providing valid estimates and for computational complexity. To handle these challenges, we have proposed a KDE that is bathymetry bounded and estimates a locally adapted kernel bandwidth h using an expanded version of Silverman's rule (Silverman, 1986) which accommodates correlated, weighted data. Computational complexity is constrained via grid-projection and pre-computation of kernels. Testing and validation of the estimator were done using synthetic simulations (see Appendix C).

2.3.1 Grid projection and pre-computed kernels

To improve computational times, we have implemented a grid-projected estimator (Sole-Mari et al., 2019). This involves obtaining a preliminary density $\stackrel{\mathrm{̃}}{\mathit{\varphi }}[i,j]$ using the histogram estimator via Eq. (B1), i.e. calculate the accumulated particle mass of all particles within each grid cell. All mass then belongs to a discrete grid, where any difference in position r between two locations of interest [i,j] and [i₀,j₀] can be expressed as $(i-i_{\mathrm{0}},j-j_{\mathrm{0}})\mathrm{\Delta }\mathit{\lambda }$.

Furthermore, we pre-compute a set of normalized kernels with fixed, discrete bandwidths given by:

$$\begin{array}{}\text{(7)}& h\left[\mathit{\omega }\right]={\displaystyle \frac{\mathit{\omega }\mathrm{\Delta }\mathit{\lambda }}{\mathrm{3}}},\mathrm{where}\mathit{\omega }=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots }\mathrm{\Omega }.\end{array}$$

Each initial non-discrete bandwidth estimate h^′ from the data-driven bandwidth algorithm is then mapped to the nearest candidate in h[ω]. Kernel support is set to ±ωΔλ beyond which the kernel contribution is set to zero. Leaked kernel mass is added back across the kernel domain using the Kernel function.

These simplifications make complexity scale with particle-containing cells instead of particles. It also allows for fast vectorized operations which drastically reduce computation time while giving negligible errors for large grids (Sole-Mari et al., 2019).

2.3.2 Data-driven adaptive bandwidth selector

The conventional Silverman's rule of thumb (Silverman, 1986) selects the optimal bandwidth h that minimizes the integrated mean square error under the assumption of Gaussian distributed data with variance σ². For a d-dimensional Gaussian kernel, one obtains

$$\begin{array}{}\text{(8)}& h={\left({\displaystyle \frac{\mathrm{4}}{d+\mathrm{2}}}\right)}^{\frac{\mathrm{1}}{d+\mathrm{4}}}{N}^{-\frac{\mathrm{1}}{d+\mathrm{4}}}\mathit{\sigma },\end{array}$$

where N is the number of data samples, i.e., the number of particles. For d=2, the expression simplifies to

$$\begin{array}{}\text{(9)}& h=N^{-\frac{\mathrm{1}}{\mathrm{6}}}\mathit{\sigma }.\end{array}$$

Here we modify Silverman's rule to yield reasonable estimates for our multi-modal, correlated, non-homogeneous, weighted data set by: (i) adapting h locally for a square shaped horizontal “adaptation” grid of size P×P surrounding each particle containing grid cell where we assume near normal, unimodal distribution, (ii) estimating the effective (uncorrelated) sample size N_eff and (iii) implementing bias corrections to a σ-estimator for weighted data-sets. The size P is determined by an integral length scale estimate (as outlined below) of the entire I×J 2-D grid. For an arbitrary cell [i,j], the adaptation grid is defined by letting l and m be discrete indices on the grid such that 1≤l, m≤P with step size Δλ in both directions. The number of particles contained in the local grid (prior to grid projection) is denoted as N_g, and the pre-computed histogram density $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$ serves as the underlying data for obtaining estimates of the local h. We will now describe the procedure of obtaining estimates of N_eff and σ to get the local h for an arbitrary adaptation grid.

Spatial correlations present in environmental data decrease the effective degrees of freedom in the sample set and we must therefore estimate and use the effective sample size N_eff to obtain a reasonable h for the local grid (e.g. Larsen et al., 2002). To estimate the effective number of spatially uncorrelated samples, a measure of correlation length is employed. We use the so-called integral length scale known from turbulence theory and statistical physics, see e.g., Yaglom (1987), Frisch (1995), and Pécseli (2000) as our objective measure of the correlation length. The correlation length L_c in terms of the integral length scale is formally defined by (Frisch, 1995)

$$\begin{array}{}\text{(10)}& L_{\mathrm{c}}={\displaystyle \frac{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left|R\right(\mathit{\varrho }\left)\right|\mathrm{d}\mathit{\varrho }}{R\left(\mathrm{0}\right)}},\end{array}$$

where $R\left(\mathit{\varrho }\right)=\mathbb{E}\left\{\mathit{\varphi }\left(x\right)\mathit{\varphi }(x+\mathit{\varrho })\right\}$ is the spatial autocorrelation function (ACF) of a continuous univariate spatial random process ϕ(x), ϱ is a spatial lag coordinate, and 𝔼{⋅} is the statistical expectation operator.

We now proceed by defining a local integral length scale L_c for the binned adaptation window. First, we estimate the local one-dimensional ACF of $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$ along each row l and column m of the square grid by using the standard unbiased ACF estimator (e.g. Percival and Walden, 1993)

$$\begin{array}{}\text{(11)}& {\displaystyle }{\widehat{R}}_l^{\mathrm{row}}\left[\mathit{\lambda }\right]={\displaystyle \frac{\mathrm{1}}{P-\mathit{\lambda }}}\sum _{m=\mathrm{1}}^{P-\mathit{\lambda }}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m+\mathit{\lambda }}\mathrm{for}l=\mathrm{1},\mathrm{2},\mathrm{\dots },P\text{(12)}& {\displaystyle }{\widehat{R}}_m^{\mathrm{col}}\left[\mathit{\lambda }\right]={\displaystyle \frac{\mathrm{1}}{P-\mathit{\lambda }}}\sum _{l=\mathrm{1}}^{P-\mathit{\lambda }}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l+\mathit{\lambda },m}\mathrm{for}m=\mathrm{1},\mathrm{2},\mathrm{\dots },P\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}$ is the binned particle density in grid cell [l,m] and $\mathit{\lambda }=\mathrm{0},\mathrm{1},\mathrm{\dots },P-\mathrm{1}$ is the discrete horizontal spatial lag index. Assuming local spatial homogeneity, we then arithmetically average the ACF estimates over all P rows and columns, respectively, to yield two one-dimensional ACF estimates as

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\widehat{R}}^{\mathrm{row}}\left[\mathit{\lambda }\right]={\displaystyle \frac{\mathrm{1}}{P}}\sum _{l=\mathrm{1}}^P{\widehat{R}}_l^{\mathrm{row}}\left[\mathit{\lambda }\right]\text{and}\\ & {\widehat{R}}^{\mathrm{col}}\left[\mathit{\lambda }\right]={\displaystyle \frac{\mathrm{1}}{P}}\sum _{m=\mathrm{1}}^P{\widehat{R}}_m^{\mathrm{col}}\left[\mathit{\lambda }\right].\end{array}\end{array}$$

We now assume local spatial isotropy and let the arithmetic average of the two perpendicular ACFs serve as a representative single ACF for the adaptation window

$$\begin{array}{}\text{(14)}& \widehat{R}\left[\mathit{\lambda }\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{R}}^{\mathrm{row}}\left[\mathit{\lambda }\right]+{\widehat{R}}^{\mathrm{col}}\left[\mathit{\lambda }\right]\right).\end{array}$$

Using the estimated ACF, we can finally estimate the local one-dimensional integral length scale ${\widehat{L}}_{\mathrm{c}}$ for the adaptation window by discretizing Eq. (10) as

$$\begin{array}{}\text{(15)}& {\widehat{L}}_{\mathrm{c}}={\displaystyle \frac{\sum _{\mathit{\lambda }=\mathrm{0}}^{P-\mathrm{1}}\left|\widehat{R}\right[\mathit{\lambda }\left]\right|\mathrm{\Delta }\mathit{\lambda }}{\widehat{R}\left[\mathrm{0}\right]}}.\end{array}$$

We now express the correlation length in terms of the associated number of samples as $N_{\mathrm{c}}={\widehat{L}}_{\mathrm{c}}/\mathrm{\Delta }\mathit{\lambda }$. It is easy to show that $\mathrm{1}\le N_{\mathrm{c}}\le P$. We then define the number of effectively uncorrelated particles N_eff as

$$\begin{array}{}\text{(16)}& N_{\mathrm{eff}}={\displaystyle \frac{N_{\mathrm{g}}}{N_{\mathrm{c}}}},\end{array}$$

and it directly follows that $N_{\mathrm{g}}/P\le N_{\mathrm{eff}}\le N_{\mathrm{g}}$. The interpretation of N_eff is straightforward: if all particles are spatially uncorrelated, then N_eff=N_g, and if all particles are fully correlated (e.g., if they are all trapped in a coherent structure), then N_eff attains its lower limit $N_{\mathrm{eff}}=N_{\mathrm{g}}/P$.

To obtain an estimate $\widehat{{\mathit{\sigma }}^{\mathrm{2}}}$ of the variance σ² in the two-dimensional binned data, we need to account for the loss of degrees of freedom due to shortening of the residual vector (Bessel's correction) and weighting as well as increased variance due to the binning process itself. The estimate of variance for the binned particle density ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}$, can then be expressed as

$$\begin{array}{}\text{(17)}& \widehat{{\mathit{\sigma }}^{\mathrm{2}}}=\left({\displaystyle \frac{{\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}\Vert {\mathit{r}}_{l,m}-\mathit{\mu }{\Vert }_{\mathrm{2}}^{\mathrm{2}}}{{\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}}}\right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathcal{B}}}\right)+{\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}},\end{array}$$

where r_l,m denotes the grid cell center point position vectors, ${\sum }_{l,m}\equiv \sum _{l=\mathrm{1}}^P\sum _{m=\mathrm{1}}^P$, and

$$\begin{array}{}\text{(18)}& \mathit{\mu }={\displaystyle \frac{{\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}{\mathit{r}}_{l,m}}{{\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}}}\end{array}$$

is the weighted mean position vector, and $\left(\frac{\mathrm{1}}{\mathrm{1}-\mathcal{B}}\right)$, where

$$\begin{array}{}\text{(19)}& \mathcal{B}={\displaystyle \frac{{\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}^{\mathrm{2}}}{{\left({\sum }_{l,m}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{l,m}\right)}^{\mathrm{2}}}},\end{array}$$

is a bias correction term that accounts for Bessel's correction and the reduced degrees of freedom due to uneven sample weights (Kish, 1965, pp. 86–88). The variance increase due to the binning process (Sheppard's correction, see e.g. Vardeman, 2005) is included through the correction term ${\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}=\mathrm{\Delta }{\mathit{\lambda }}^{\mathrm{2}}/\mathrm{12}$.

The final local bandwidth estimate (for each particle-containing grid cell) then follows from Eq. (9):

$$\begin{array}{}\text{(20)}& \widehat{h}=N_{\mathrm{eff}}^{-\frac{\mathrm{1}}{\mathrm{6}}}\widehat{\mathit{\sigma }}.\end{array}$$

2.3.3 Boundary solution

We establish a boundary solution for the density estimator by interpolating bathymetry data onto the model grid across all predefined depth layers, creating a matrix of “permissible” and “impermissible” cells for gas. The boundary control is imposed at the kernel estimation stage before summation, by directly modifying kernels whose support contains impermissible cells (see Fig. 2). While being computationally intensive, this greatly simplifies the boundary control and entirely omits the difficulties of finding a reliable boundary solution that handles the complex bathymetry and physical processes appropriately.

Figure 2 Sample kernel with support $i,j=\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{10}$, containing impermissible cells in its kernel support. Density is indicated by blue shading and impermissible cells (Bathymetry/land) and blocked cells where mass from the shown kernel is not permitted to access are shown as black and yellow colored cells, respectively. Cells identified as being in the line of sight between the kernel center (green dot) and cell [1,8], according to Bresenham's line algorithm, are grayed out (although there is nothing particular about this line).

Impermissible cells are treated as impenetrable obstacles, reflecting that dissolved gas cannot cross land or shallow bathymetric boundaries. Any density within, or “blocked” by impermissible cells, is considered misplaced. A cell is defined as blocked if it lacks a clear line of sight to the kernel center. We determine line of sight using Bresenham's line algorithm (Bresenham, 1965). This is an efficient incremental algorithm relying solely on integer arithmetics that identify grid cells located between an origin cell (x₀,y₀) and a target cell (x₁,y₁) (Fig. 2). The algorithm is implemented on a normalized grid with unit cell lengths and initialized by first defining the direction, or step coefficients s_x and s_y in the x and y directions, respectively:

$$\begin{array}{}\text{(21)}& s_x=\left\{\begin{array}{ll}\mathrm{1},& \text{if }{x}_{\mathrm{0}}<x_{\mathrm{1}}\\ -\mathrm{1},& \text{if }{x}_{\mathrm{0}}>x_{\mathrm{1}}\end{array}\right.,s_y=\left\{\begin{array}{ll}\mathrm{1},& \text{if }{y}_{\mathrm{0}}<y_{\mathrm{1}},\\ -\mathrm{1},& \text{if }{y}_{\mathrm{0}}>y_{\mathrm{1}}.\end{array}\right.,\end{array}$$

and an error term $\mathit{\epsilon }=\mathrm{d}x-\mathrm{d}y$, where $\mathrm{d}x=|x_{\mathrm{1}}-x_{\mathrm{0}}|$ and $\mathrm{d}y=|y_{\mathrm{1}}-y_{\mathrm{0}}|$. The algorithm then iteratively updates (x₀,y₀), using ε to determine whether to step in x or y direction (s_x, s_y are not updated) via the following criteria:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \mathrm{2}\mathit{\epsilon }>-\mathrm{d}y:x_{\mathrm{0}}\Rightarrow x_{\mathrm{0}}+s_x\mathrm{and}\\ & \mathrm{2}\mathit{\epsilon }<\mathrm{d}x:y_{\mathrm{0}}\Rightarrow y_{\mathrm{0}}+s_y.\end{array}\end{array}$$

At each time step, (x₀,y₀) is added to the list of grid cells, thereby iteratively forming the line of sight to the target cell. Density in blocked or impermissible cells are redistributed to permissible cells according to the kernel function.

2.4 Atmospheric flux and mass modification functions

Changes in dissolved gas content due to processes within or at the boundaries of the water column are included by modifying the particle masses. Total mass change of a particle ζ is estimated at each time-step n using predefined mass modification functions that couple particle properties to the gridded field processes. Here, our modeling framework is relatively flexible, and can even accommodate models where it is necessary to keep track of higher order parameters, such as microbe stocks (e.g. a Monod model). This is made possible since we can model any parameter explicitly across the domain. We will only describe mass modification due to dissolved atmospheric exchange of gas here, but a mass modification function for microbial oxidation of CH₄ is presented in the application section (Sect. 3).

Atmospheric flux can be implemented following any theory using the surface layer concentration as input data. Here we propose a simple solution by applying the bulk equation from Wanninkhof (2014):

$$\begin{array}{}\text{(23)}& F=\mathit{\kappa }\left({\mathrm{\Phi }}_{\mathrm{atm}}-{\mathrm{\Phi }}_{\mathrm{sw}}\right),\end{array}$$

where F [mol m⁻² s⁻¹] is the gas flux across the sea-air interface, κ [m s⁻¹] is the gas transfer velocity, and Φ_atm and Φ_sw [mol m⁻³] are atmospheric and surface water concentrations, respectively. The gas transfer velocity κ can be expressed as

$$\begin{array}{}\text{(24)}& \mathit{\kappa }(U_{\mathrm{a}},T_{\mathrm{s}})=C_{\mathit{\kappa }}{U}_{\mathrm{a}}^{\mathrm{2}}{\left({\displaystyle \frac{\mathit{\text{Sc}}\left(T_{\mathrm{s}}\right)}{\mathrm{660}}}\right)}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

where Sc(T_s) and C_κ are empirically derived constants and U_a is the wind speed 10 m above the sea surface. The Schmidt number Sc(T_s) is the empirically derived, gas-specific temperature-dependent ratio between sea water kinematic water viscosity and the diffusion coefficient of the gas. The C_κ coefficient lumps together a set of various processes that govern sea/air exchange and has been determined for CO₂ and a wind speed range of $\mathrm{4}<U_{\mathrm{a}}<\mathrm{15}$ m s⁻¹ using inverse modeling for global estimates. Validity for other gases and wind ranges is not fully known.

Let the gridded estimate of the 2D spatiotemporal atmospheric flux field $\mathit{\beta }(x,y,t)$ be denoted $\widehat{\mathit{\beta }}[i,j,n]$, using the same horizontal and temporal grid cells as $\widehat{\mathit{\varphi }}[i,j,k,n]$. We then assume an initial equilibrium between the atmospheric concentration and background surface concentration, which is disturbed by the (modeled) seep-derived dissolved gas. The difference between surface water and atmospheric concentration in Eq. (23) is then simply the surface layer (k=0) concentration $\widehat{\mathit{\varphi }}[i,j,\mathrm{0},n]$. To obtain an estimate of the gas transfer coefficient κ, we project re-analysis atmospheric 10 m above sea level wind speed and sea surface temperature data onto all $i,j,$ and ns, delivering the gridded gas transfer coefficient field estimate $\widehat{\mathit{\kappa }}[i,j,n]$ using Eq. (24). The gridded atmospheric dissolved flux field estimate is then given by

$$\begin{array}{}\text{(25)}& \widehat{\mathit{\beta }}[i,j,n]=\widehat{\mathit{\kappa }}[i,j,n]\widehat{\mathit{\varphi }}[i,j,\mathrm{0},n]\mathrm{\Delta }{\mathit{\lambda }}^{\mathrm{2}}\mathrm{\Delta }t,\end{array}$$

where $\widehat{\mathit{\beta }}[i,j,n]$ is the integrated atmospheric flux from grid cell $[i,j,n]$.

Loss of gas due to atmospheric flux is implemented by modifying the mass of all particles present in the surface layer. To ensure efficient computation and mass conservation, we assume that the entire contribution to the atmospheric flux from a surface layer particle occurs within the grid cell where that particle resides, disregarding the effects of mass distribution through the density kernels. Errors associated with this assumption are expected to be small, since wind and temperature fields and consequently, gas transfer velocities are generally smooth on typical kernel bandwidth scales. It is also mass conserving, because atmospheric flux varies linearly with dissolved gas concentration (Eq. 25). Furthermore, since grid cell concentration depends linearly on the total cell gas content (i.e., the sum of all particle masses in that cell) and the gridded gas transfer velocity, relative flux contributions from particles can be estimated using products of particle masses and cell specific gas transfer velocities. The mass loss due to atmospheric exchange for a surface layer particle α at time-step n can then be expressed as

$$\begin{array}{}\text{(26)}& {\mathit{\gamma }}_{\mathit{\alpha }}\left[n\right]={\displaystyle \frac{{\mathrm{\Gamma }}_{\mathit{\alpha }}\left[n\right]\widehat{\mathit{\kappa }}\left[c\right(\mathit{\alpha }),n]}{{\sum }_{\mathit{\zeta }\in \mathcal{A}}{\mathrm{\Gamma }}_{\mathit{\zeta }}\left[n\right]\widehat{\mathit{\kappa }}\left[c\right(\mathit{\zeta }),n]}}\sum _{i=\mathrm{1}}^I\sum _{j=\mathrm{1}}^J\widehat{\mathit{\beta }}[i,j,n],\end{array}$$

where α∈𝒜 and 𝒜 denotes the set of all surface-layer particles, and c(α) denotes the indices i,j where particle α resides.

3 Application

We applied the modeling framework to a well documented natural CH₄ seep site offshore northwestern Norway located in the Hola trough (Fig. 3), where coral reefs and CH₄ seeps coexist (Chand et al., 2008). These seeps were investigated not only to assess the mechanisms governing CH₄ fluxes to the atmosphere, but also to evaluate their potential impact on cold water coral ecosystems (Sert et al., 2025; Argentino et al., 2025). A thorough description of the data, site characteristics, environmental conditions, and seabed flux estimates are presented in Ferré et al. (2024). In short, the observed seeps are weak, and our focus is therefore on examining system dynamics and fractional distribution of gas, rather than on quantifying environmental impacts or contributions to the atmospheric CH₄ budget.

Figure 3 Bathymetric map of the application area and location of Conductivity Temperature Depth (CTD) station, observed seep-associated flares indicated by yellow and pink dots during the 20–22 May 2018 survey. Seeding locations (where particles in the particle trajectory model is released), estimated as the flux weighted average position of the incorporated seeps, are indicated by the yellow and pink triangle (see Sect. 3.2). Coloring reflects which seeding location each seep observation is pooled into.

We modeled the resulting direct and diffusive atmospheric gas release, as well as 3D concentration from 45 observed CH₄ seeps for the period between 20 May and 20 June 2018. A 1-month period was chosen since it captures a relatively wide range of periodic variability in both ocean and atmospheric circulation patterns and yields relatively modest computation times. The OpenDrift simulation required 2–3 d on a supercomputer and the concentration modeling 5–6 h on a workstation laptop.

3.1 Seep gas phase modeling

Free and dissolved gas profiles and direct free atmospheric gas flux were modeled individually for each of the 45 seeps using M2PG1 (see Sect. 2.1) to steady state, using observed and inferred input data and settings, as outlined in the following sections.

3.1.1 Gas phase modeling input data and settings

Temperature and salinity data were extracted from a Conductivity Temperature Depth (CTD) cast performed in 20 May 2018 (Figs. 3 and 4a). Seabed gas flux for each seep was estimated using single beam echosounder data (Simrad EK-60 scientific SBE splitbeam echosounder) obtained between 20 and 22 May 2018 and are presented in Ferré et al. (2024). All other input parameters had to be inferred as outlined in the following paragraphs since we lack observations.

Figure 4 Input parameters used in the M2PG1 model runs. (a) Conservative temperature and absolute salinity (at the CTD station obtained in 20 May 2018. (b) Bubble size distribution (Veloso et al., 2015) and bubble rising speed (Fan and Tsuchiya, 1990) for different bubble sizes as a function of effective radius $r_{\mathrm{E}}=(a^{\mathrm{2}}b{)}^{\mathrm{1}/\mathrm{3}}$. (c) Bubble flatness and surface area for various bubble sizes as function of effective radius for Jansson et al. (2019), Leblond et al. (2014), and Spherical flatness parametrization. Note that the linear flatness (Jansson flatness) appears non-linear in the figure since its linearity is with major spheroid axis and not effective radius.

M2PG1 requires an initial bubble size distribution, and we used the polynomial fit to visual observations of bubbles as presented in Veloso et al. (2015). Note that since M2PG1 takes into account the non-spherical shape of bubbles, the bubble size distribution is given using the effective radius $r_{\mathrm{E}}=(a^{\mathrm{2}}b{)}^{\mathrm{1}/\mathrm{3}}$, where a and b are the major and minor axis of the spheroid, respectively (Fig. 4b). We used the bubble rising speed model from Fan and Tsuchiya (1990) using their recommendation for bubble contamination, and the linear flatness parametrization from Jansson et al. (2019) (see Fig. 4b and c). We describe and discuss the bubble rising speed model and deformation parametrization selection in Appendix D.

The horizontal domain size was determined using Appendix A and an assumed barotropic current of $\stackrel{\mathrm{‾}}{U}=\mathrm{0.1}$ m s⁻¹, horizontal diffusivity D_h=0.01 m² s⁻¹, rising speed 〈w〉=0.25 m s⁻¹ and a H=200 m deep water column. We estimated σ_w using the bubble rise spectrum (Fig. 5) to σ_w=0.025 m s⁻¹. This resulted in an estimated model area of A_ℳ∼88 m² and grid cell side-lengths 9.4 m. Vertical and temporal resolution does not affect grid cell concentration but must obey the Courant-Friedrichs-Lewy numerical stability condition. Here we use a grid cell height of 1 m to obtain $A_{\mathcal{M}}^{⟂}=\mathrm{9.4}$ m² and a time-step of 0.0625 s.

Figure 5 Probability density for bubble rising speed of bubbles in the bubble plume using the bubble rising speed model from Fan and Tsuchiya (1990) and bubble size distribution from Veloso et al. (2015).

We assumed a constant vertical mixing coefficient of D_v=0.001 m² s⁻¹ and background dissolved gas concentrations were set to the default values from Jansson et al. (2019) at $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{4}}$ mol L⁻¹, $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ mol L⁻¹, $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}$ mol L⁻¹, $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ mol L⁻¹ and $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{8}}$ mol L⁻¹ for N₂, O₂, CO₂, CH₄, and Ar, respectively.

We provide an overview of microbial oxidation rate coefficient (k_ox) observations presented in literature and their associated uncertainty in Appendix E and Fig. 6. Here, we use the simple average $k_{\mathrm{ox}}=\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{7}}$ s⁻¹ of the full compiled dataset in Table E1 which include cold seep environments, hydrothermal vents, and human-made releases.

Figure 6 Maximum CH₄ oxidation rate coefficients (k_ox) obtained from datasets found in literature and detailed in Appendix E and Table E1. The x-axis is logarithmic, meaning that the bars cover different ranges, i.e. the histogram bars are narrower at smaller scales. Vertical dashed lines indicate simple average and median of the values in the table. For the application offshore Northwestern Norway, we used the average of all compiled values.

3.1.2 Gas phase modeling results

Most of the CH₄ gas is dissolved in the water column, with concentration appearing to decrease near exponentially towards the sea surface (Fig. 7a). Hourly seabed gas flow rate was ∼97 mol of which 93 % dissolved below 100 m depth and only ∼0.28 % reaching the atmosphere. Integrated atmospheric release from free gas over the 1-month period was 183.1 mol. Free CH₄ gas content closely follows the total free gas content throughout the water column (Fig. 7a) and loss of total free gas volume (bubble shrinkage, collapse, and dissolution) dominates over other gases replacing CH₄ in bubbles. The resulting change in dissolved gas profiles for the four other gases (N₂, O₂, CO₂ and Ar) due to bubble transit, i.e. the transport of gas molecules by entering bubbles, rising, and subsequently dissolving at shallower depths, was therefore negligible (never exceeding 0.1 % of background values). Atmospheric flux from the 45 seeps varied considerably from $<{\mathrm{10}}^{-\mathrm{7}}$ to $\sim {\mathrm{10}}^{-\mathrm{5}}$ [mol s⁻¹] (Fig. 7c), mostly due to large variations in seabed fluxes (Ferré et al., 2024). Dissolved gas injection rates, which are needed as input in the particle dispersion modeling step, were calculated using the modeled (by M2PG1) dissolved gas profiles (not shown) and Eq. (1) and are shown for the 45 seeps in Fig. 8a.

Figure 7 (a) Vertical profiles of total free gas content for the 5 gases (colored axes) and the total free gas (black, lower x-axis) in the water column for all seeps combined. (b) Distribution of gas content on bubble sizes for all seeps and gases combined at different depths (note power in scale). (c) Free atmospheric gas flux at the sea surface for the 45 seeps (logarithmic color scale).

Figure 8 (a) Dissolved CH₄ release rate profiles (q_m) for all observed seeps on a logarithmic scale for the two cluster groups (see Fig. 3 for seep/group locations, group west in yellow and group East in pink), and (b) Resulting accumulated dissolved CH₄ release rate from each groups (lines, upper x-axis) and histogram of released particles at each modeling time-step (hourly). The smaller bottom bar at the western cluster reflects the slightly shallower depth in this area.

3.2 Particle data set

Using OpenDrift, we simulated a particle data set of N particles with associated 4D-positions (3D space and time) for the 1-month period, as described in Sect. 2.2.

The advective and diffusive components for the drift model (Eq. 2) were determined using velocity vectors and diffusivity coefficients (throughout the water column) obtained from the NorKyst-800 hydrodynamic modeling system. NorKyst-800 is based on the Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams, 2005) framework and is a terrain-following, free-surface, primitive equations model with 35 vertical layers and a horizontal resolution of 800 m (Albretsen et al., 2011). The model is eddy-permitting and can resolve major fjord systems and other coastal bathymetric features such as troughs. Vertical turbulent exchange is computed using the General Length Scale closure scheme (Umlauf and Burchard, 2003). OpenDrift has an option to automatically check if diffusivity coefficients are reasonable from a physical perspective, and we used a non-zero fallback value of D=0.2 m² s⁻¹ when OpenDrift deemed the input data unphysical. OpenDrift did not excessively use the fallback value, and from tests where we increased and decreased the value, we found that the choice of fallback value had negligible impact on the results.

We seeded 500 particles every time-step with a particle lifetime of 4 weeks and simulated particle trajectories by updating their positions every 5 min. The effect of vertical mixing on the particle trajectory was modeled using a sub-timestep of 30 s. We configured OpenDrift to store the particle positions every full hour during the simulation.

The initial particle mass was calculated using seabed gas flux data and Eq. (3). To saturate the particle field prior to the study period, the simulation was initiated on 20 April (1-month spin-up time). With this setup, total particle mass in the domain would increase in the first 10–15 d of the study period due to the re-distribution of removed particle mass. However, particle count would remain approximately constant, with Λ ∼ 330 000 particles present at each time-step.

Due to the wide range in seep intensity and relatively closely clustered seep positions (Fig. 7) combined with the limitation of 500 release particles each time-step, we chose to aggregate the seeps into two seeding locations to promote smooth release profiles and reduce round-off effects. Grouping was done based on visual inspection of the seep positions (Fig. 3) and seed locations were calculated using the flux-weighted average position of each group, given by:

$$\begin{array}{}\text{(27)}& \mathit{r}={\displaystyle \frac{{\sum }_{g\in \mathcal{G}}{\mathrm{{\rm Y}}}_g{\mathit{r}}_g}{{\sum }_{g\in \mathcal{G}}{\mathrm{{\rm Y}}}_g}},\end{array}$$

where 𝒢 denotes the set of all seep indices included in the seed position. The seed locations and their associated seeps are indicated by matching colored triangles and dots in Fig. 3. The 500 particles were then distributed according to the injection rate profiles for each seed location with a 1 m vertical resolution, resulting in the profiles/histograms shown in Fig. 8b.

3.3 Concentration and atmospheric flux estimation

Concentration and atmospheric fluxes were estimated on a $[i,j,k,n]$ grid with cell sizes Δλ=800 m, Δz=25 m, and Δt=3600 s, covering the time period between 20 May and 20 June and geographical region between (12.5° E, 68.5° N), (12.0° E, 72.1° N), (21° E, 72° N), and (20.1° E, 68.45° N). Although particle data were technically available outside of this region, we chose this boundary to avoid potential edge effects in the hydrodynamic model and to constrain computation time. Kernel bandwidths were estimated for each cell location in the 4D grid (each [$i,j,k,n$]) using a P×P sized adaptation grid where we determined P using an integral length scale estimate of every 2D (i,j) layer of the particle data. The size of P, typically varying between ∼ 7000 and ∼ 20 000 m, agreed reasonably well with observations and theory on meso-scale eddy sizes in the region (Dugstad et al., 2021). Boundary conditions were implemented as described in Sect. 2.3.3 using [i,j]-interpolated IBCAO v. 4 bathymetry data (Jakobsson et al., 2012). Spatiotemporal gas transfer velocities $\widehat{\mathit{\kappa }}[i,j,n]$ were estimated from grid interpolated ERA V reanalysis wind and sea surface temperature data (Hersbach et al., 2023) which, together with the surface layer concentration estimates $\widehat{\mathit{\varphi }}[i,j,\mathrm{0},n]$, gave the atmospheric flux field estimates $\widehat{\mathit{\beta }}[i,j,n]$ using Eq. (25).

Particle mass was adjusted at each time-step using the mass modification terms (γ) for atmospheric flux (Eq. 26) and redistribution from removed particles (Eq. 5). We also added a mass modification term for microbial oxidation, a crucial process when simulating the evolution of dissolved CH₄ content in the ocean (Appendix E). We used a simple first order kinetics formula (Eq. E1), with the same rate coefficient $k_{\mathrm{ox}}=\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{7}}$ s⁻¹ as in the gas phase model (see Appendix E and Fig. 6 for the determination of k_ox). Microbial oxidation was then included by imposing a mass loss γ_ox=k_oxΔtΓ_ζ[n] at each time-step. In principle, this corresponds to discretization of Eq. (E1) using a standard first order forward finite difference scheme. Mass modification of any particle at any time step could then be calculated by summing up the three applied mass modification terms: (i) mass loss to microbial oxidation, (ii) mass loss to the atmosphere and (iii) mass gained from nearby removed particles.

3.4 Application results

The results from the gas phase modeling step is shown and described in Sect. 3.1.2 and the dissolved concentrations, atmospheric flux and fate of CH₄ molecules in the modeled domain in the following sections. Animations of the time varying 3D CH₄ concentration field and 2D diffusive release field are shown in Video S1 (https://doi.org/10.5446/69942, Dølven, 2025c) and Video S2 (https://doi.org/10.5446/69941, Dølven, 2025d), respectively. It is important to note that the modeling results presented are subject to a wide range of relatively uncertain assumptions concerning various model coefficients and that the main aim here is to test the modeling framework and investigate the dynamics of the system, rather than conclude about absolute values.

3.4.1 3D CH₄ concentration field

The averaged distribution pattern of CH₄ throughout the study period is strongly affected by generally northeastward currents that transports gas along the coast, following the shelf and shelf break. The gas enters the more open fjord systems, and to a lesser degree inner fjords. North of 70° the CH₄ plume disperses more, branching into a northward plume leaving the coast and a coastal plume that keeps following the coastline. The concentration anomaly is generally small, around 2–4 orders of magnitude lower than typical oceanic CH₄ background concentration values ($\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ mol m⁻³, Fig. 9), due to the weak seabed release.

Figure 9 Nine depth layers of modification (i.e. the estimate $\widehat{\mathit{\varphi }}$ of the anomaly ${\mathrm{\Phi }}^{\prime }(x,y,z,t)$)) to the CH₄ concentration on 1 May as indicated on top of each panel. Typical background concentration in the ocean is $\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ mol m⁻³ for reference. The bathymetric boundary for the different layers are delineated with a grey contourline.

The 3D concentration field is very dynamic due to the energetic regional current regimes, and shows variability on ranging from tidal (∼ 12 h) to fortnightly periods. A visual representation of the temporal variability of the top 9 layers in the water column (down to 200 m depth) is shown in Video S1.

Figure 10 Fractional vertical distribution of CH₄ in the water column 28 d after release. Color scale shows the fraction of total depth integrated CH₄ in the water column.

Most of the CH₄ is displaced upward relatively quickly from the trough and pushed on top of the shelf break (Fig. 9). Vertical distribution of CH₄ is therefore characterized by a quick (a few days) shift from CH₄ being mainly located close to the seafloor at the release site (∼ 150–200 m depth) to shallower depths (Fig. 10). A thorough analysis of mechanisms causing the rapid shift in location of CH₄ is outside the scope of this study, however, upwelling within troughs and along the shelf break is well documented in this region (e.g. Slagstad et al., 1999).

3.4.2 Diffusive CH₄ flux to the atmosphere, microbial oxidation and non-physical redistribution and loss

The time-integrated 2D diffusive atmospheric CH₄ release over the study period is shown in Fig. 11 and the complete time series can be found in Video S2. Within the model domain, most of the CH₄ remains in the water column or is consumed by microbes, with a total of ∼ 0.76 % (∼ 528 mol) being exchanged diffusively with the atmosphere. This diffusive release is roughly three times greater than the local free gas release (0.27 %) and does not account for any diffusive release occurring outside of the model domain. The diffusive flux extends across a broad area spanning several hundred kilometers and shows pronounced temporal variability that is strongly correlated with wind speed (Fig. 13a), although no clear effect of surface water depletion is observed after storm events. Microbial CH₄ consumption exceeds atmospheric flux by a factor one to two orders of magnitude, emphasizing its crucial role in regulating dissolved CH₄ levels (Fig. 13b). Loss from the domain due to particles leaving the model domain is also substantial (∼ 5 to ∼ 50 times the atmospheric loss) and shows clear tidal excursions patterns with periodic variability (Fig. 13c). Non-physical methane loss due to removed isolated particles (i.e. too far away to be re-distribution) is on the same order of magnitude as loss to the atmosphere via diffusive release (Fig. 13d).

Figure 11 Modeled accumulated diffusive release of CH₄ within the model domain from the seeps between 20 May and 20 June 2018.

3.4.3 Fate of released CH₄

We analyzed the vertical redistribution and partitioning of CH₄ among available sinks within the model domain over a four week period (the particle lifetime). This analysis is important not only for evaluating CH₄ molecules potential to reach the atmosphere, but also in cases impacts on water column and/or seafloor ecosystems are of interest. Excluding removed particles and particles leaving the model domain, the accumulated fractional water column CH₄ loss due to atmospheric exchange shows an exponential increase the first couple of days, with a subsequent near linear slope until the end of the four week period. The initial rapid gradient increase in atmospheric loss fraction corresponds to a vertical redistribution of dissolved CH₄, where the concentration maximum shifts from ∼ 200 to ∼ 10 m depth (Fig. 10). After four weeks, around 0.7 % of dissolved CH₄ molecules had been transferred to the atmosphere (Fig. 12a and b). Microbial oxidation dominates over both atmospheric diffusive and free-gas fluxes transforming 65 %, while around 34 % of the CH₄ remains in the water column at the end of the particle lifetime.

Figure 12 (a) Accumulated fractional loss to atmospheric exchange and microbial oxidation and fraction CH₄ that remains in the water column and (b) accumulated fractional loss to atmospheric exchange and microbial oxidation in days after release.

Figure 13 (a) Loss of CH₄ from the water column from atmospheric equilibrium, (b) microbial oxidation, (c) particles leaving the model domain and (d) mass loss due to deactivation of particles that are unable to redistribute its mass.

3.4.4 Interpretation of results

Diffusive exchange of CH₄ exceeds the local free gas release and spreads over a large ocean region, making it almost impossible to detect and quantify using conventional measuring instrument. This also poses a challenge for atmospheric inversion models, since these are better at detecting point sources rather than weak releases over large regions (Thompson and Stohl, 2014). These limitations highlight the uncertainty in quantifying the impact of seabed seepage on the atmospheric CH₄ budget, particular when considering the potential increased seepage in recent decades due to e.g. thawing marine permafrost, hydrate dissociation (e.g. Serov et al., 2015; Ruppel and Kessler, 2017), and anthropogenic disturbances of the seafloor (e.g. drilling).

Although here, the estimated total atmospheric flux is small, the impact of more extensive seepage (such as e.g. Mau et al., 2017) could be significant and at the same time difficult to observe and/or trace.

Even though the dissolved gas spreads out over cold water coral reef areas, the CO₂ generated by microbial oxidation is likely too small to have any measurable effect on the local ocean environment and cold water corals. This primarily reflects the weak seabed fluxes. For more intense and/or localized seepage, e.g. a leaking gas well, this might not be the case. It is also worth noting that the influence of seabed gas seepage on cold-water coral ecosystems remains a sparsely explored field of research.

An additional major caveat, both regarding atmospheric fluxes and potential impact on the ocean environment, is the uncertainty caused by the microbial oxidation rate coefficient k_ox assumption in methane oxidation rates (MOx), which vary by several orders of magnitude and correspond to half-lives for dissolved methane (or methane turnover) ranging from 5 d to nearly 2 years (Fig. 6, Table E1). To examine sensitivity to k_ox, we conducted a limited coefficient-sweep experiment, rerunning the framework using the lowest and highest reported “cold seep” rate coefficients from Table E1: $\mathrm{0.02}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹ (low) and $\mathrm{0.98}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹ (high) (Gründger et al., 2021), as well as two intermediate values.

The low and high rate coefficient runs resulted in a ∼ 34 % increase (705 mol) and ∼ 41 % decrease (309 mol), respectively, in atmospheric emissions during over the model domain and particle lifetime. The impact on the final fate of dissolved CH₄ molecules were also considerable (Table 1). After 4 weeks, ∼ 5.5 % of dissolved CH₄ were consumed and 1.16 % released to the atmosphere in the low rate coefficient run, compared to ∼ 91 % consumed and 0.36 % released in the high rate coefficient run. These results highlights the importance of selecting a suitable MOx rate coefficient and illustrates the huge span of results one can obtain when coefficients in the modeling framework are poorly constrained. It is also important to note that since MOx is a biologically mediated process, it can vary substantially on relatively small spatiotemporal scales (Valentine et al., 2001; Ruff et al., 2015) depending on a wide range of factors, including CH₄ concentration, water temperature, salinity (Steinle et al., 2015; Osudar et al., 2015), nutrient availability (Knief, 2015), and the presence of trace elements (Hanson and Hanson, 1996). Thus, assuming a constant rate coefficient is in itself a potentially problematic simplification, since it will most likely vary considerably across the model domain. One would, for instance, expect MOx to decrease with the distance from seep area due to rapid dilution of methane and varying environmental stress for the CH₄ oxidizing microbes. To reduce the uncertainty concerning MOx rates, future studies must therefore not only constrain rate coefficients, but also improve our ability to model MOx dynamics within modeling frameworks. Including more complex MOx parametrizations is possible in our framework since we allow explicit modeling of higher order fields at each time-step and location in the modeling domain.

Table 1 Fate of CH₄ 4 weeks after it is dissolved in water using a sweep of different microbial oxidation rate coefficients (k_ox).

Another important source for uncertainty in the modeled fate of CH₄ arise from the atmospheric bulk model. The atmospheric gas transfer coefficient function (Eq. 24) was derived based on, and designed for global CO₂ estimates (Wanninkhof, 2014) and has an estimated uncertainty of ∼20 %, even for its intended use. We must expect considerably higher uncertainties due to our application in a local coastal region (as opposed to an ocean basin), where wind speeds may exceed the validity range, and since the specific coefficient C_κ has been determined solely for CO₂.

Uncertainties in the eddy diffusivity, vertical transport and distribution are also expected to be large. The choice of grid cell thickness can also modify the end result. If the grid cells are too thin, and temporal resolution too coarse, there is a risk of depletion of the surface layer between the model output time-steps. On the other hand, if the grid cells are too thick, one would incorporate CH₄ from depths where exchange with the atmosphere is unrealistic, thereby violating the assumptions of the atmospheric exchange bulk model. One can evaluate whether the surface thickness is sufficiently thick by comparing typical values for Eq. (26) with the typical mass of surface layer particles and ensure that the atmospheric loss is considerably smaller than the surface layer gas content (i.e. that γ_α[n]≪Γ_τ[n] always).

4 Conclusions

We implemented and applied a new framework for modeling the impact of seabed gas seepage on spatiotemporal water column concentrations and atmospheric gas exchange with the ocean. The application uncovered and highlighted important aspects of the dynamics regarding the fate of seeped gas from the seabed, such as a highly distributed diffusive release which considerably exceeds local free atmospheric gas fluxes.

Estimation uncertainties arise from a relatively wide range of sources which should be addressed in future studies. In particular, mass loss due to microbial oxidation pose a significant challenge since rate coefficients are shown to exhibit large variability which cause considerable differences in the modeled atmospheric fluxes and concentrations. Current parameterizations of mass loss due to atmospheric ventilation are also simple and developed for large scale ocean regions, not coastal areas. Our results are also sensitive to poorly resolved diffusivity coefficients, which can greatly affect the water column distribution of dissolved gas. A steady state assumption for the seepage itself might also be a problematic assumption in areas where seepage is known to vary strongly over time (e.g. Ferré et al., 2020).

From a pure modeling perspective, non-physical re-distribution of dissolved gas is also a source for potential error and when estimating the final fate of gas, the limitation in model domain size makes inference about final state difficult. In cases where the aim is to estimate the final fate of released gas, one could argue that it is more suitable to use a 1D approach, e.g. as suggested in Nordam et al. (2025). Using a steady state solution of the gas phase model is also a drawback and might currently cause significant errors, especially for intense seep sites where background gas concentrations can be significantly altered.

On the other hand, the framework is flexible and reasonably fast and makes it possible to employ complex, locally adapted existing hydrodynamic models and can include advanced process modules, thereby capturing not only idealized processes but also complex hydrodynamic and chemical/biological phenomena. It has a wide range of potential applications, not only for monitoring known gas seeps, but also for risk assessments concerning future potential increased seepage due to e.g. hydrate thawing (James et al., 2016), leaking gas wells, and integrity of subsea legacy carbon storage reservoirs (e.g. Torsæter et al., 2024) and other leaking industrial installations. Certain studies also requires a 3D spatiotemporal concentration field, e.g. when studying the potential effect of seepage on biological processes in a specific area. Aspects of the framework (e.g. the kernel density estimator) can also be complementary to established frameworks for post-processing and analysis of ocean particle dispersion data and contaminant spreading in the ocean in general (e.g. the ChemicalDrift module, Aghito et al., 2023).

Aside from improving process and atmospheric exchange modules, future developments should consider a full on-line coupling between concentration model and gas phase model as well as a proper validation study to ensure realistic results.

Appendix A: M2PG1 Model Grid Cell Dimensions

Here we propose a solution for determining a reasonable horizontal model grid cell size assumption for M2PG1, as no established method currently exist. Since M2PG1 assumes horizontally invariant concentrations within the predefined model domain, the choice of the horizontal model domain size directly affects the concentration within the model grid cells and, in turn, gas transfer and dissolution. Defining the horizontal dimensions of the model grid cells must therefore be done with care and should reflect the horizontal extent of the modeled bubble plume to obtain realistic results. We determine the horizontal and vertical gas phase model grid cell area, respectively denoted A_ℳ and $A_{\mathcal{M}}^{⟂}$, by modeling the 2-dimensional spread of the bubble cloud.

We assume that the seeps are point sources and that the bubbles drift with a barotropic current with mean speed $\stackrel{\mathrm{‾}}{U}$ and random velocity fluctuations governed by a horizontal diffusivity D_h. In this framework, horizontal bubble spread is caused by (i) differences in accumulated horizontal displacement resulting from varying rising speeds of bubbles with different sizes (slow/fast bubbles spend more/less time in the velocity field) and (ii) turbulent (random) effects in the horizontal flow, modeled as diffusion. The horizontal extent of the bubble plume increases towards the sea surface and we use the estimated spread at half of the total water column depth H to minimize estimation errors at the surface/bottom.

The spread due to differences in rising speed can be estimated using the probability density P for bubble rising speeds w_o in the bubble cloud. The distribution of rising speeds for bubbles in the bubble cloud can be described by the discrete probability P[w_o], and can be derived from the chosen initial discrete bubble size distribution (BSD) and bubble rising speed model. This is done by estimating the bubble rising speed of all bubbles in the BSD and re-bin the results, using the fractional weights from the BSD, according to bubble sizes. We obtain the weighted distribution average and standard deviation as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& \langle w\rangle ={\displaystyle \frac{{\sum }_o{w}_{o}P\left[w_o\right]}{{\sum }_{o}P\left[w_o\right]}}\text{and}\\ & {\mathit{\sigma }}_{\mathrm{w}}=\sqrt{{\displaystyle \frac{{\sum }_o{\left(w_o-\langle w\rangle \right)}^{\mathrm{2}}P\left[w_o\right]}{{\sum }_{o}P\left[w_o\right]}}},\end{array}\end{array}$$

where w_o are discrete rising speeds, P[w_o] associated probabilities, 〈w〉 the weighted average, and σ_w weighted the standard deviation (see Fig. 5 for an example). Note that the BSD is expected to change with height above the seafloor (which also changes P[w_o]). For the purpose of this calculation, however, we assume the BSD remains unchanged.

Along-flow spread Δx_rs can then be expressed as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& \mathrm{\Delta }{x}_{\mathrm{rs}}=\stackrel{\mathrm{‾}}{U}\mathrm{\Delta }{t}_{max},\mathrm{where}\\ & \mathrm{\Delta }{t}_{max}={\displaystyle \frac{H}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{1}}{\langle w\rangle -{\mathit{\sigma }}_{\mathrm{w}}}}-{\displaystyle \frac{\mathrm{1}}{\langle w\rangle +{\mathit{\sigma }}_{\mathrm{w}}}}\right].\end{array}\end{array}$$

Horizontal displacement due to current diffusivity acts in both along-flow (x) and cross-flow (y) direction and can be expressed by the 2D Gaussian solution to the diffusion equation for a point source,

$$\begin{array}{}\text{(A3)}& p(x,y,t)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{4}\mathit{\pi }{D}_{\mathrm{h}}t}}}{e}^{-(x^{\mathrm{2}}+y^{\mathrm{2}})/\mathrm{4}{D}_{\mathrm{h}}t},\end{array}$$

where $p(x,y,t)$ is the normalized count of bubbles at position (x,y) and time t and 4D_ht is the variance of the spread in both directions. We constrain diffusive spread using twice the standard deviation 2σ_D of the distribution at $H=H/\mathrm{2}$ given by

$$\begin{array}{}\text{(A4)}& \mathrm{\Delta }{x}_{\mathrm{D}}=\mathrm{\Delta }{y}_{\mathrm{D}}=\mathrm{2}\sqrt{\mathrm{2}{D}_{\mathrm{h}}}\cdot \sqrt{\mathrm{0.5}{t}_{\mathrm{H}}}\mathrm{where}{t}_{\mathrm{H}}={\displaystyle \frac{H}{\langle w\rangle }}.\end{array}$$

and

$$\begin{array}{}\text{(A5)}& A_{\mathcal{M}}=(\mathrm{\Delta }{x}_{rs}+\mathrm{\Delta }{x}_D)\mathrm{\Delta }{y}_D\end{array}$$

giving horizontal grid cell side lengths of $\sqrt{A_{\mathcal{M}}}$ (since M2PG1 uses square cells).

An estimate of the vertical grid cell area, which is needed to estimate dissolved gas injection profiles is easily obtained and defined as

$$\begin{array}{}\text{(A6)}& A_{\mathcal{M}}^{⟂}=\sqrt{A_{\mathcal{M}}}\mathrm{\Delta }{z}_{\mathcal{M}},\end{array}$$

where z_ℳ is the vertical grid cell size.

Appendix B: The histogram estimator

A commonly used density estimator based on data from particle dispersion models is the histogram estimator. The histogram estimator for the concentration estimate $\widehat{\mathit{\varphi }}$ at position r₀ using a predefined grid with grid cell volume V can be expressed as

$$\begin{array}{}\text{(B1)}& \widehat{\mathit{\varphi }}\left({\mathit{r}}_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{1}}{V}}\sum _{\mathit{\zeta }=\mathrm{1}}^Z{\mathrm{\Gamma }}_{\mathit{\zeta }}K({\mathit{r}}_{\mathrm{0}},{\mathit{\eta }}_{\mathit{\zeta }}),\end{array}$$

where Γ_ζ and η_ζ represent the mass and positions (respectively) of particles, and

$$\begin{array}{}\text{(B2)}& K({\mathit{r}}_{\mathrm{0}},{\mathit{\eta }}_{\mathit{\zeta }})=\left\{\begin{array}{ll}\mathrm{1}& \text{when }{\mathit{\eta }}_{\mathit{\zeta }}\text{ shares the same}\\ & \text{grid cell as }{\mathit{r}}_{\mathrm{0}},\\ \mathrm{0}& \text{otherwise}.\end{array}\right.\end{array}$$

Using the histogram estimator implies modeling a smooth, continuously distributed property with a discontinuous, quantized, and piece-wise constant function, which introduces several drawbacks with this estimator-property pairing. Firstly, the estimator is highly dependent on the choice of grid cell size: fine grids result in noisy and unrealistic estimates in regions with medium to low particle counts, while coarse grids lead to significant loss of information in areas with high particle counts. Secondly, the histogram estimator can be sensitive to the chosen position of the origin. In addition, the minimum concentration estimate is limited to one particle per grid cell, which can significantly influence e.g. atmospheric flux estimates (for instance if that concentration exceeds the atmospheric background concentration).

Some of these issues can be mitigated by adjusting the grid cell size, however, the problems prevail in highly heterogeneous domains containing regions with low particle saturation (unless one adobts an unstructured grid). Seeding more particles is always a remedy, however, we are still left with inefficient use of the particle position data and potentially unfeasible computational complexity (see Sect. 2.2.2). We have therefore formulated an adaptive bandwidth, 2D grid-projected Kernel Density Estimator (KDE) specifically for OpenDrift output data to calculate the concentration field.

Appendix C: Density estimator testing and validation

The adaptive kernel density estimator was developed, tested, and compared with other estimators using a numerical toy model that generates data resembling typical OpenDrift data. The toy model gives full control of all parameters and allows to efficiently test various scenarios due to the low computational cost of each run. Here we compare the adaptive bandwidth KDE against other estimators as explained below.

C1 Toy model and test simulation

The synthetic data was designed to mimic output data from OpenDrift by seeding particles at seed location r₀ and calculating their position at time step T (at time TΔt) as

$$\begin{array}{}\text{(C1)}& {\mathit{r}}_T={\mathit{r}}_{\mathrm{0}}+\sum _{\mathit{\vartheta }=\mathrm{0}}^{T-\mathrm{1}}\left({\mathbf{U}}_{\mathit{\vartheta }}\mathrm{\Delta }t+{\mathit{\xi }}_{\mathit{\vartheta }}\sqrt{\mathrm{2}D\mathrm{\Delta }t}\right),\end{array}$$

where ${\mathbf{U}}_{\mathit{\vartheta }}=[u_{\mathit{\vartheta }},v_{\mathit{\vartheta }}]$ represents a spatially uniform, time varying velocity field, $\mathit{\vartheta }=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },T$ and Δt is the modeling time-step, ${\mathit{\xi }}_{\mathit{\vartheta }}=[{\mathit{\xi }}_{\mathit{\vartheta },x},{\mathit{\xi }}_{\mathit{\vartheta },y}]$ where each component is sampled from a standard normal distribution 𝒩(0,1) and D represents a diffusivity coefficient (assuming isotropy). The velocity U_ϑ was calculated as

$$\begin{array}{}\text{(C2)}& {\mathbf{U}}_{\mathit{\vartheta }}={\displaystyle \frac{{\mathbf{U}}_{\mathrm{0}}+{\mathbf{\Xi }}_{\mathit{\vartheta }}}{\Vert {\mathbf{U}}_{\mathrm{0}}+{\mathbf{\Xi }}_{\mathit{\vartheta }}\Vert }}\Vert {\mathbf{U}}_{\mathrm{0}}\Vert \end{array}$$

where ${\mathbf{U}}_{\mathrm{0}}=(u_{\mathrm{0}},v_{\mathrm{0}})$ gives the initial velocity, ${\mathbf{\Xi }}_{\mathit{\vartheta }}=\left(\sin \right[\frac{\mathit{\vartheta }}{\mathrm{50}}]u_{\mathrm{0}},v_{\mathrm{0}})$ and $\left|\right|\cdot \left|\right|$ gives the euclidean norm. The normalization with $\Vert {\mathbf{U}}_{\mathrm{0}}\Vert$ is necessary to ensure conservation of mass in the field. We choose $U_{\mathrm{0}}=(\mathrm{0.1},\mathrm{0})$ and D=0.14 and released a total of 2×10⁶ particles from a point source at ${\mathit{r}}_{\mathrm{0}}=(\mathrm{10},\mathrm{10})$ over the course of 400 timesteps. The histogram density estimate (Eq. B1) of the full simulation was considered the simulation “Ground-truth” and is shown in Fig. C1. The computation time for generating the test data was 8.3 s with a Intel Core Ultra 9 185H processor.

Figure C1 Left: Synthetically generated particle dispersion model data for 2×10⁶ particles (purple) and randomly picked “test data” of 2000 particles (green dots) in a domain with a simple “impermissible” elliptic boundary. Right: Histogram estimate of the full 2×10⁶ particle dataset, representing the “ground truth” in the test scenario for the adaptive kernel density estimator.

Figure C2 Density estimates using the test dataset (N=2000, green dots in Fig. C1) using, from left to right, a histogram estimator, a Silverman (non-adaptive) bandwidth estimator from the scipy.kde package, a Time-dependent bandwidth kernel density estimator (TKDE) with boundary control, and the Adaptive bandwidth kernel density estimator (AKDE) developed here. The upper panel figures show the density estimates and the lower panel figures the residuals from the ground truth estimate shown in Fig. C1. The impermissible region (land/bathymetry) is shaded in grey.

C2 Testing and evaluation of different estimators

We implemented and tested four estimators (i) The histogram estimator, (ii) An Time-dependent bandwidth estimator, (iii) The Silverman bandwidth estimator from the gaussian_kde function from the scipy.kde python package, and (iv) The adaptive bandwidth estimator used in the present study. All estimators were tested on a data set where we picked every 1000th particle from the full data set of 2×10⁶ resulting in a total of 2000 particles for density estimation. All estimates were done grid-projected as described in Sect. 2.3.1 and for the final time-step only and all particles had a mass of 1.

For the histogram estimator estimate, we used Eq. (B1) and for the time-varying bandwidth estimator, we defined the bandwidth as

$$\begin{array}{}\text{(C3)}& h_{\mathrm{tv}}=\sqrt{\mathrm{4}D{t}_{\mathit{\vartheta }}}\end{array}$$

which is the theoretically ideal bandwidth for the time-varying estimate. In a real-world scenario, the diffusion coefficient varies and we cannot estimate the correct diffusion coefficient unless information about the local diffusivity is given from the hydrodynamic model. Although the bandwidth function could be suitable when such information is available, complex bathymetry may introduce challenges as discussed below. For the estimate using the provided scipy.kde.gaussian_kde, we used default settings and the bw_method='silverman' setting (SciPy Community, 2024). We refer to the package documentation for details (URL in the reference list). For the adaptive bandwidth estimator, we followed the algorithm described in Sect. 2.3.2 and estimated h locally for each particle-containing grid cell. For the “in-house” coded estimators (all but scipy.kde.gaussian) we included the boundary control explained in Sect. 2.3.3.

A comparison of the four estimators and a residual analysis plot are shown in Fig. C2. In addition to a visual comparison, we evaluate how the max values in the field aligns and do a simple R² statistic. The Histogram estimator gives a noisy result, with a very high max value of 9000 compared to 5672 for the ground truth and a low R²=0.53. The the non-adaptive Silverman results in an unrealistically smooth estimate, with a very low maximum value of 927 and low R²=0.59. While the time varying bandwidth estimator works relatively well in the open “unbounded” part of the domain, it over-smooths when encountering the boundary – highlighting a problem with time varying bandwidth in bounded domains where the stochastic process is limited by physical obstacles. Nonetheless, it performs better than the non-adaptive Silverman and Histogram estimators, achieving an R²=0.71. The adaptive bandwidth estimator is in general slightly over-smooth, however, it significantly outperforms the three other estimators with a maximum value of 4531, which is the closest to the ground truth of 5672 and has a high R²=0.90 (Fig. C2).

The total computation time for doing all the KDE estimates (including the kernel adaptation) were less than 1 second with a Intel Core Ultra 9 185H processor, and the adaptive kernel density estimator was only slightly slower than the gaussian_kde function from the scipy.kde package (both $\le {\mathrm{10}}^{-\mathrm{3}}$). A simple comparative performance study as well as script for further testing and evaluation of the adaptive kernel density estimator and the estimators used for comparison is available in (Dølven and Hanssen, 2025) and corresponding GitHub repository (https://github.com/KnutOlaD/akd_estimator/releases/tag/v1.2.0, last access: 12 November 2025).

Appendix D: Rising speed model and flatness parametrization

It is well-known that the terminal bubble rise velocity U_b varies non-linearly with the bubble size (e.g., Fan and Tsuchiya, 1990; Leifer and Patro, 2002). In addition, it depends on fluid and gas parameters, and the degree of contamination. Fan and Tsuchiya (1990, Eq. 2.11) showed that the terminal bubble rise velocity can be written as

$$\begin{array}{}\text{(D1)}& U_{\mathrm{b}}={\left(U_{\mathrm{b}\mathrm{1}}^{-c}+U_{\mathrm{b}\mathrm{2}}^{-c}\right)}^{-\mathrm{1}/c},\end{array}$$

where U_b1 dominates for small bubbles, and U_b2 dominates for large bubbles. The dimensionless parameter c (called n in Fan and Tsuchiya, 1990 and d in Leifer and Patro, 2002) is a measure of the degree of contamination, which directly affects the surface tension of the bubbles. By fitting Eq. (D1) to multiple experimental data sets, Fan and Tsuchiya (1990) found that $\mathrm{0.8}\le c\le \mathrm{1.6}$, where the lower limit corresponds to contaminated bubbles, and the upper limits corresponds to clean bubbles. We follow their recommendation, and apply c=1.2 for our assumed moderately contaminated conditions, see Fig. 4b and Jansson et al. (2019) for further details.

Bubble deformation is an important factor in bubble dissolution and exchange rates of gas since it changes the surface area to volume ratio of the bubbles. Deformation can be characterized by a dimensionless flatness ratio, defined as $f\equiv a/b$. In addition to spherical flatness, two parametrization options are available in M2PG1: Leblond flatness (Leblond et al., 2014), where $f=\mathrm{0.45}+\mathrm{1.4}\ln (a/b_{\mathrm{ref}})$ for a>1.48 and Jansson flatness (Jansson et al., 2019), where $f=\mathrm{1}+\mathrm{0.3064}(a/b_{\mathrm{ref}})$ (coined here, referred to as “linear flatness” in Jansson et al. (2019) and b_ref=1 mm. While Jansson flatness parametrization has support for bubbles where a<1.48 mm, it lacks an empirical basis other than a fair agreement with Leblond for relatively small bubbles. The divergence between the two models at larger bubble sizes (Fig. 4c) can lead to misrepresentations when modeling distributions that are skewed towards larger bubbles. Nonetheless, we use Jansson flatness parametrization since here our observations indicate that most of the gas is confined to smaller bubble sizes (Ferré et al., 2024).

Appendix E: Brief review of existing estimates of MOx rate coefficients

Oxidation of CH₄ to carbon dioxide is achieved by several groups of aerobic methanotrophs and reaction rates vary substantially, depending on existing microbial consortia, stoichiometry of the involving nutrients, and overall succession of the methanotrophs (Hanson and Hanson, 1996). Nonetheless, reaction rate measurements with radiotracer assays highlight first-order reaction kinetics and a CH₄ decay rate following

$$\begin{array}{}\text{(E1)}& {\displaystyle \frac{\mathrm{d}\mathit{\varphi }\left(t\right)}{\mathrm{d}t}}=-k_{\mathrm{ox}}\mathit{\varphi }\left(t\right)\Rightarrow \mathit{\varphi }\left(t\right)={\mathit{\varphi }}_{\mathrm{0}}{e}^{-k_{\mathrm{ox}}t}\end{array}$$

where k_ox is the reaction rate and ϕ₀=ϕ(0) the initial concentration. Most measurements of microbial CH₄ oxidation (MOx) in marine environments are focused on locations where CH₄ concentrations exceed the background levels. The rates of CH₄ oxidation in suboxic zones, hydrothermal vents, and cold seeps exhibit substantial variability, spanning several orders of magnitude from 10⁻⁸ to 10⁻² nM s⁻¹, primarily due to spatiotemporal fluctuations in CH₄ concentrations. In contrast, half-life or CH₄ oxidation rate constants (k_ox) are independent of CH₄ concentration and provide a more accurate representation of the water column's MOx capacity. The rate coefficients (k_ox) range from $\mathrm{0.02}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{1.74}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹, corresponding to halving times of approximately five days to two years. However, CH₄ can remain stable for decades in oxygen-limited environments where aerobic CH₄ oxidation is inhibited.

Ward et al. (1987)Ward et al. (1989)Pack et al. (2015)de Angelis et al. (1993)Valentine et al. (2010)Steinle et al. (2016)Sansone and Martens (1978)Mau et al. (2012)Steinle et al. (2017)Mau et al. (2020)Weinstein et al. (2016)Uhlig et al. (2018)Sert et al. (2023)de Groot et al. (2024)Gründger et al. (2021)Gründger et al. (2021)Gründger et al. (2021)Gentz et al. (2014)Mau et al. (2017)Mau et al. (2017)Sert et al. (2020)Mau et al. (2013)

Table E1 Methane Oxidation Rate Coefficients (k_ox) in units of 10⁻⁶ s⁻¹ (µHz) from various studies. We have obtained the maximum k_ox reported in the studies unless ranges are given. Half-lives are calculated by solving for ϕ(t_0.5)=0.5ϕ₀ in Eq. (E1), i.e. $t_{\mathrm{0.5}}=\ln \left(\mathrm{2}\right)/k_{\mathrm{ox}}$. In Gründker et al. (2021) only May data was included from 2016 and the difference in turnover time between 2016 and 2017 is because the maximum rate coefficient was $\mathrm{1.85}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻¹ in 2016 and $\mathrm{2.01}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻¹ in 2017, but this difference is rounded off in the table.