A very good piece of work, I am struggling to find issues with the paper and work that has been conducted.

First, in the interest of transparency I would like to mention that one of the authors is at present working in the same department as me. I informed the editor of this, and he said he would take this information into account when assessing the review.

Overall, I think the paper is interesting, and the basic approach outlined is sound. I think this would be a usfeful addition to the literature, where many studies have only focused on the immediate bubble transport, and been more vague about the fate of the dissolved methane (which is after all the greater part).

My main issues with the paper are first that I think the choice of rate coefficient for microbial oxidation could be discussed further. According to appendix E, half-lives found in the literature vary from 5 days to two years, and then to just consider a single rate coefficient which corresponds to a half-life of around 30 days seems a bit simple. I'm not saying you necessarily need to do more simulations, as I think the point of the paper is more to demonstrate a method than to give definite numbers, but I think this could be discussed further.

Second, the paper is quite long, and dedicates a lot of time to describing a partially new method for calculating concentrations. I am not completely convinced by the authors' arguments that this was necessary, and even if it was (necessary), I think much of Section 2.3 could go in the appendix, as it distracts from the main point. (In my opinion, the main point of the paper is that it's a good idea to combine a bubble model with a 3D transport-fate model for dissolved methane. The rest is just details.)

Finally, I have a large-ish number of more detailed comments below, mostly minor suggestions.

I recommend that the paper be published, subject to revisions.

Line 7: In what sense was the modelling "successful"? Did you compare against measurements or other estimate, or do you just mean that the model(s) ran successfully and produced reasonable output?

Line 22: "Atmospheric measurements are currently the only approach..." Surely the current study is based on the assumption that also modelling is a possible approach to estimating atmospheric emissions from seep areas? 

Line 33: Alien gases?

Line 41: What is "in situ data"?

Lines 90 and onwards: Personally, I find the notation with the square brackets a bit odd, though I suppose it is a matter of taste. And why use subscripts for particle number and square brackets for time? How about subscripts for particle and superscript for time? Or two subscripts? (There's lots of nice inspiration to be had from the world of general relativity: https://en.wikipedia.org/wiki/Christoffel_symbols#Definition_in_Euclidean_space )

Lines 95-100: Strictly speaking, this is not a correct description of the equation of motion for particles in OpenDrift. It solves an SDE, not an ODE, and there is no such thing as the diffusive velocity vector in the scheme used in OpenDrift. (If you try to calculate the diffusive velocity, you will find that it goes to infinity as dt ->0). The correct equation can be found in numerous papers, see for example Eq. (2) in Spivakovskaya et al. (2007): https://link.springer.com/article/10.1007/s10236-007-0102-9

Lines 107-115: I think these expressions are not quite consistent. On line 107 you say that the mass of particle \zeta, at time n, \Gamma_\zeta[n], is given by the previous mass and some functions. In Eq. (3), you say that \Gamma_\zeta[n] is defined by the release. At this stage, I am also confused about how you represent different release rates for different sites, and how you deal with the vertical variation in dissolution from the bubbles, if all particles are seeded with the same amount of mass, but presumably that will be described later?

Line 122: Who refers to this as "particle death"? Very dramatic term, what's wrong with "particle removal"?

Line 130: Isn't there something wrong in Eq. (5)? Surely if \gamma_s is supposed to be based on the distance between particle s and particle \theta, then the position of particle s would be expected to appear on the right-hand side of the equation?

Line 148: If I understand correctly, you said that you found histogram unsuitable, but you are essentially using a histogram in the vertical direction. Why not use 3D KDE?

Line 150: Just a comment, but I think KDEpy is pretty efficient. Might be easier (and maybe faster?) than making a new implementation.

Line 153: Whether the density estimate is differentiable or not depends on the choice of kernel.

line 157: I don't think the variable V is defined.

Line 162: I wouldn't worry about the kernel being consistent with diffusive transport. And in any case, if you wanted the kernel to be consistent with diffusion, wouldn't you have to let the bandwidth grow based on the diffusivity? And you also later truncate the kernel.

Lines 199, 202, and elsewhere: Two different notations: N_{eff} and N_\mathrm{eff}

Line 260: It's probably not important in the grand scheme of things, but redistributing the mass instead of mirroring fail in some simple test-cases. For example for uniform distribution on a bounded 1D domain, it will fail to reproduce a uniform distribution, but rather giving too low density near the boundaries.

Line 276: Wouldn't U_a more commonly be the wind at 10 m height?

Lines 292-296: Do I understand correctly that any particle loses mass based on the sum of the flux \beta over the entire horizontal grid? So even a particle in an area with no wind will lose mass based on the average loss?

Line 318: Polynomial fit to visual observations of what?

Line 333 and Figure 6: What type of average was used?

Line 336: Just out of curiosity: Does the deposited concentration actually decrease exponentially, or does it just look approximately exponential? Obviously it cannot be a _true_ exponential if it depends on environmental variables like temperature, but even for a theoretical case of a perfectly homogeneous water column, is it actually exponential?

Line 343: What are the units of the atmospheric fluxes?

Line 348: 1,2,3,N, ? Missing some dots? But also, no need to repeat the concept of the indexing, but rather mention how many particles were used.

Lines 352 and 353: No need to repeat terrain-following

Line 357: "therefore" ... strictly speaking it does not follow logically from the fact that the model is used by the authorities for acute happenings, that it is also the best model when running a modelling study several years later. (Note: I'm not saying that NorKyst isn't the best model, only that it doesn't follow logically from the statement)

Line 392: Missing units on the rate coefficient.

Line 393: How is the "mass modification term" used? And shouldn't it be negative?

Line 405: Missing closing bracket. Also, I believe mol/L (molar) is a more common unit than mol/m3.

Lines 415-427: I'm curious about the thickness of the vertical layers, and particularly the top layer, in the model used to calculate the flux from dissolved to the atmosphere. If this layer is to thin, the layer might be depleted on a faster timescale than mixing from the lower layers, since that mixing is modelled by particles that live in a different world, with a different timestep. On the other hand, if the layer is too thick, it will allow the escape of methane from too deep in the water column. In a hybdrid model like this, where the escape to the atmosphere and the transport of the particles are modelled by different approaches, this seems to me to be a point worth investigating and discussing.

Line 423: Reads like the microbial flux is driven by wind speed.

Line 424-425: Do I read correctly that the loss due to particles leaving the domain is of the same order of magnitude as the loss to microbial oxidation? You have used a fairly high biodegradation rate, which is perfectly fair I think, but some studies claim to have found much slower oxidation rates (see for ekample https://pubs.acs.org/doi/full/10.1021/acs.est.7b02732). How would you use your model to study the case of slower biodegradation, without the loss due to particles leaving the domain being completely dominant?

 

Lines 456-457: I would have thought that also the vertical eddy diffusivity used in the transport model is quite important. You say that the mass transfer coefficient is stated to have an uncertainty of 20%, my guess would be that the vertical eddy diffusivity from the ocean model, and the effect of that on the transport of methane towards the surface, has a (much) larger uncertainty than 20%. Any comments?

Lines 490-502: I'm not sure I understand how the M2PG1 model works originally, or the implications of choosing a horizontal cell size. How can the original model take ambient concentrations into account for dissolution of methane from the bubbles, and for mass transfer of dissolved methane to the atmosphere, without calculating a concentration? And it also seems to me that the dissolution mostly happens near the sea floor, where the plume is narrow, while the mass transfer to the atmosphere only happens at the surface, where the plume is wide, so maybe the cell size should vary with depth? Finally, I wonder how sensitive the results are to this choice. Does much of the mass-transfer from the atmosphere happen in the M2PI1 model, or is that only calculated based on the concentration fields from the particles? And to what degree does the ambient concentration in the water column hinder the dissolution of methane from the bubbles?

Line 533: Should be "piece-wise constant function"

Line 536: I don't undestand why the histogram estimator is "highly sensitive" to the position of the origin. I get that the results can look different if you shift the origin, but "highly sensitive"?

Line 550 and 561: What do the different numbers of particles (1900000 and 10^6) refer to?

Line 551: What does \mathcal{Z}^{100 \times 100} mean? To me it suggests some esotheric 10000 dimensional discrete space. Perhaps you mean i,j \in [1, 2, ... 100]?

Line 553: Isn't there something odd in Eq. (C1)? The sum is over $t$, but there is no $t$ in the expression. Also, the units do not match. You mention a timestep later, but it should be included here to obtain displacement from velocity, and the random vector should be scaled by sqrt(dt) to match the units of sqrt(D).

Lines 565-567: Again, I would suggest including KDEpy (https://kdepy.readthedocs.io/en/latest/) in the comparison. I'm not sure how it will compare, but in my limited experience I've found it quite fast. Also, I'm not sure if Silverman's rule is a very relevant option, I believe it is well known to be less than optimal for multi-modal distributions.

Line 573: "theoretically ideal". Only a minor detail, but I'm not sure if this would be theoretically ideal (depends on the theory, I guess). In a particle model, where the particles are transported by the diffusivity, this increases the variance of the overall distribution. But the kernel bandwidth also increases the variance of the overall distribution. I seem to rememeber that it is easy to show (at least of all particles have the same kernel and bandwidth) that the variance of the distribution obtain by KDE is the sum of the variance of the particle positions, and the variance of the kernels. So I wonder if this "theoretically ideal" bandwidth wouldn't in fact double-count the diffusion? You could test it against an analytical solution of the diffusion equation in 1D.

Line 590: Some words missing at the end here?

Figure C2: The histogram example uses very small cells (10000 cells and only 1000 particles), it would probably be closer to the "ground truth" with larger cells.

Line 634: I suggest a manual proofread as well, there are still a few typos here and there (including on this line).

Just one thing I forgot to mention in the review comments, that you could perhaps also discuss: What are the advantages of 3D modelling compared to 1D modelling? You mention some of the drawbacks (considerable computational effort, loss of mass through particles exiting the domain), so it would be interesting to discuss if it is worth the effort.

Oh, and one more small point: I think the units on the colorbar in Fig. 9 are wrong, the numbers seem to be in mol/L not mol/m3?

The authors present a numerical modeling setup for simulating the chemical, physical, and biological processes associated with natural gas seeps, dissolved gas in the water column, and its eventual release to the atmosphere. The model has been applied to simulate seeps on the Norwegian continental shelf. Please find my specific comments and questions on the manuscript below:

Lines 128: “The redistribution is weighted according to the inverse  distance from the dying particle within a user defined distance limit dmax,:  How does the diffusive transport of particles that is adjusted by this re-distribution of mass and affect the model results?

Page 5: last line: “This solution changes the problem of non-physical loss of dissolved gas to one of non-physical redistribution, which is generally considerably less problematic.” 

How do you say this is less problematic? Based on any evidence or proof?

Lines 35-40 “Our aim is to provide a framework which can integrate all key processes governing free and dissolved transport and transformation of 

seeped gas to provide a full 3-d concentration field in the water column and total atmospheric release estimates.” There has been previous modeling studies in literature focusing on these aspects. See the studies https://pubs.acs.org/doi/full/10.1021/acs.estlett.3c00493, https://www.nature.com/articles/s41467-024-53780-7, https://sintef.brage.unit.no/sintef-xmlui/handle/11250/2730544 , https://pubs.acs.org/doi/full/10.1021/acs.est.5c03297

Line 345 “modeling step, were calculated using the dissolved gas profiles” Are these modeled or measured?

Line 505: In the model presented the vertical binning size is dependent on the bubble rising speeds wo. The bubble size which controls the rising speed changes with pressure changes and the mass transfer at different depths. Hence the bubble size distribution (BSD) changes at different depths. BSD at what depth level was used when you decided the vertical bin size and why?

Line 550 : “OpenDrift by seeding N = 1900000 particles“. This is a very large number of particles that were used in the simulations. How did you decide (what is the basis) the number of particles to be used? and how much computer resources needed for these simulations and time taken. It would be good to give an idea of this to the readers. 

Line 620: In this  section where you describe the oxidation rates you are presenting several ranges. For example rates of CH4 oxidation 

Varying  from 10−8 to 10−2 nM s−1,) and ((kox) range from 0.02·10−6 to 1.74·10−6 s−1,. You should state which numbers were used in your study in the Norwegian waters and justify the reasons for choosing these numbers as there is a large variation.

Section 4 Conclusions 

I expect large variability in the biodegradation rates will affect the results presented to a large extent. I believe should be discussed in the manuscript and included in the conclusion. Please see the recent studies presented in https://sintef.brage.unit.no/sintef-xmlui/handle/11250/2730544 and https://pubs.acs.org/doi/full/10.1021/acs.est.5c03297