We would like to extend our sincere gratitude to the reviewers for their valuable and constructive comments on our manuscript. Their insights have been instrumental in refining and enhancing the quality of our work. We carefully considered each comment and suggestion, and we are pleased to submit this revised version of the manuscript with thorough revisions addressing the raised concerns. The reviewers' feedback has significantly contributed to the clarity and depth of our research, and we believe that their input has ultimately strengthened the rigor and validity of our work. We are confident that this collaborative effort will result in a more impactful contribution to the scientific community. We thank the reviewers once again for their time and expertise, and we look forward to their further evaluation of our manuscript.

Comments from Referee #1:

-    E.g. in order to computationally bridge the gap between pore scale and macro scale, it is crucial to be able to scale up the system size. It is well-documented that OpenFoam is well parallelized and given that reactions are local, so should BRNS. However, the paper seems to lack a description of strong or weak scalability. Could the authors please elaborate on this?

We added more explanation on this topic in the revised manuscript on lines 309-313.

-    In at least one of the test cases 3.1--3.3, it would be useful with a discussion on how grid resolution affects the convergence to the analytical solution.

Convergence of the numerical solution for different grid sizes are addressed for test case 3.2 (lines 406-410).

-    Are osmotic forces included in the model? There is no force term in Eq. (2) due to gradients in the chemical potential, however it could be accommodated by a re-interpretation of the pressure P. Would P3D-BRNS be able to model e.g. droplet motion due to a concentration gradient and solubility difference in the two phases?

Any other forces that are not mentioned in the manuscript are not considered. However, our framework is pretty modular which makes it easy to add/consider other forces in the future. We added a statement regarding this issue on lines 581-583.

-    section 3.1: How was the mesh generated? Is it skewed or staircase-like at the diagonal side? How is the sharp top corner handled? Some more information here would be useful.

Explanation regarding the mesh generation is added on line 336.

-    Caption of Fig. 1 also inserted at line 115.

We removed this (line 119).

-    line 359: incomplete sentence

The sentence has been revised.

-    line 412: simluations --> simulations

Changed as suggested (line 443).

-    line 487: growth --> grow

Changed as suggested (line 521).

We would like to extend our sincere gratitude to the reviewers for their valuable and constructive comments on our manuscript. Their insights have been instrumental in refining and enhancing the quality of our work. We carefully considered each comment and suggestion, and we are pleased to submit this revised version of the manuscript with thorough revisions addressing the raised concerns. The reviewers' feedback has significantly contributed to the clarity and depth of our research, and we believe that their input has ultimately strengthened the rigor and validity of our work. We are confident that this collaborative effort will result in a more impactful contribution to the scientific community. We thank the reviewers once again for their time and expertise, and we look forward to their further evaluation of our manuscript.

Comments from Referee #2:

-    Firstly, the presentation of the model is, in my opinion, inadequate. I would consider myself an "informed reader", but found it very difficult to understand many parts of the model description. I concede that  due to the complexity of the model, it is not an easy task to present it in an easily accessible manner. I have written down many questions for the authors regarding the description of assumptions, equations, limitations etc. which I hope the authors see as constructive criticism of the presentation of the work.

We thank the reviewer for his constructive comments. We think the revised manuscript now provides a clearer description of our model.

-    Secondly, the authors have not satisfactorily separated previous work form new contributions.  Consider the fact that the presented model is a coupling between InterFoam for flow and BRNS for reactive transport. It is clear that the authors have not simply plugged these two together to produce their simulator. What is not clear is what they have done. Example: A simple Internet search of InterFoam lead me to this website:

https://openfoamwiki.net/index.php/InterFoam

which is basically most of Section 2.1.  Yes, not all, but most of the section. While I understand the authors' wish to describe the physical model "from scratch", I am afraid the manuscript reads like a repetition of previous work. That is one thing. And the other thing is that, as mentioned already, the authors do not really explain many of the basic principles well. I cannot say what the best way to solve this problem is, but it should be addressed. At the very least, the authors should CLEARLY highlight what THEIR contribution to the model is, and CLEARLY state which part of the presented model is simply adopted from previous work. The part that has been adopted can then be described in a brief manner, while the part that is new should be described rigorously.

We added a paragraph explaining which parts of the existing OpenFOAM package we are making use of (lines 259-263). We do think that a brief overview on the underlying equations and basic assumptions as provided in Sections 2 should be given in the manuscript. We tried to clarify our statements to improve our explanations, but for concepts taken from the literature we rather refer to the original references for an in-depth explanation.

-    Line 57: "... as well as nonlinear relation of local biomass concentration and the bulk nutrient concentration."

There is no relation between these two, whether linear or non-linear. How biomass concentration and bulk-nutrient concentration are related to each other for a particular problem can only be determined by the solution of a boundary value problem. I don't understand what the authors want to express here.

We agree with the reviewer and reworded our statement to clarify that we refer to changes of these concentrations (lines 59-61). 

-    Line 67 and Line 73: "For the continuum scale a ..." and "... the soil at the microscale, the ..." 

Expressions like microscale and continuum scale mean something different for each reader. I suggest the authors either use domain-specific terms like "pore scale" and "Darcy scale" or define how they use scale-related terms at the beginning of the manuscript.

Throughout the manuscript, all instances of “microscale”, and “macroscale” or “continuum scale” are changed to “pore scale” and “Darcy scale”, respectively.

-    The authors' excessive use of parentheses diminishes legibility. I advise against structuring sentences around parentheses. In one particular case on Lines 76--79, the authors have even opted to using nested parentheses which further reduces the reader's ability to comprehend the intended message and which has resulted in the authors forgetting to close the brackets.

This passage has been revised avoiding excessive use of parentheses (lines 80-83). Throughout the revised manuscript, the use of parentheses has been reduced and the use of nested parentheses has been avoided.

-    Lines 88--89:

Why not the Stokes equations for Re<<1? In any case, the NS equations should be plural. 

The Navier-Stokes equations are chosen for a couple of reasons: 1) the velocities we are dealing with satisfy very low Mach numbers, and 2) the Navier-Stokes equations are more general compared to the Stokes equation, that in turn serves well the purpose of having easier time in extending the model in the future (line 117-118). The NS equations are now given in plural (line 93).

-    Eqs. (1) and (2):

This is the incompressible form of the NS equations. Do the authors really assume that the gas phase is incompressible?

The pressure gradients (and Reynolds numbers) considered in our simulations are in the range that changes in gas compressibility are negligible (lines 139-140).

-    Eq. (5):

How does one obtain the individual velocities of the fluids? How does one obtain the velocity of the interface?

Explanation for individual velocities are added on lines 128-130.

-    Line 141:

How are u_alpha und u_beta different from u_i?

Alpha and beta are indices of u_i; an explanation is now provided on line 142.

-    Line 144: “The latter term in equation (7) is active only in the presence of an interface.”

There are three terms in Eq. (7), not only two.

As far as I can tell, all three terms are only active in the presence of an interface. If div u =0 (Eq. 1), then the second term is also zero.

The confusing statement has been removed (L 147).

-    Lines 145--146: “… which computationally helps with maintaining stiffness of the interface.”

It really is not clear what the stiffness of the interface is.

Complementary explanation on interface stiffness/sharpness has been added on lines 148-150.

-    Eq. (8):

What is the max operator acting on? The magnitude of velocity?! Which is the set of values for which the maximum is to be determined?

How does Eq. (8) give the relative velocity between the phases?

u_c is not the relative velocity but referred to as the compressive velocity. All the variables, as mentioned in the text, are global attributes rather than phase attributes. So the individual velocity and hence the relative velocity between phases are put of the scope of all the equations. The explanation for the max operator is added on line 154-155.

-    Lines 150—151: “In simulation scenarios introduced in this paper, c_alpha has been assigned the value of 1, unless otherwise is stated.”

So, if c_alpha =1, what does this mean for Eq. (8)? How does it simplify?

Further explanations on the term c_alpha has been added on lines 156-158.

-    Line 154: “… kappa as mean interface curvature, …”

Same question as with the max operator, what is the mean here performed over? Is it the mean curvature for each interface or for all interfaces in the domain or …?

Mean interface curvature is calculated in each computational grid cell (now clarified on line 161).

-    Line 155: “… interfacial tension force …”

I presume this force should act tangentially to the interface. How is this enforced? sigma and kappa are scalars. grad alpha is a vector which, as far as I can tell, is orthogonal to the interface

“…Interfacial tension force…” model used in our manuscript is referenced in lines 164-165 and and its explanation/derivation is not in the focus of our paper

-    Eq. (19):

u_n should be defined.

A definition of u_n has been added on line 236.

-    Eq. (21):

How is this interface normal vector used in the simulations?

The use of the interface normal vector at the triple point is clarified now on line 245-246.

-    Lines 242—245: “In short, the modifications they proposed and which are used here are: 1) smoothing the indicator function to have a better measure of the interface curvature, 2) sharpening the indicator function for computation of the interfacial tension force, 3) filtering the capillary pressure force parallel to the interface, and 4) filtering capillary fluxes based on the capillary pressure gradient.”

It is next to impossible for most readers to parse what these four steps actually mean in terms of the mathematical formulation presented here.

What does it mean to smoothen and then sharpen the indicator function? What is meant by "filtering" here? How is the capillary pressure force calculated? What is it defined as? What is a capillary flux? How is it obtained mathematically?

We used the enhanced interFoam solver developed by Raeini’s work directly in our model. We settled with pointing out the key factors added by this solver from Raeini and refer the interested readers for more details to the reference provided. The purpose of the filtering is briefly explained on lines 257-258 , and the meaning of smoothing interface is explained before on line 155-156).

-    Line 254:

How much of the mathematical formulation is simply inherited from InterFoam and how much is contributed by the authors? This is not clear. There is no point in repeating at least in part the description of a module of openFOAM which already exists.

The contribution of our work is now explicitly described on lines 259-263. As mention above, we think a brief introduction of all governing equations and basic assumptions should be provided in the manuscript. 

-    Lines 378—379: “… , strong depletion of the tracer in aqueous phase as well as the concentration jump across the interface are observed.”

A better explanation is needed here.

An improved explanation of this topic is now given on lines 401-405.

-    Line 394: “The bacteria residing in the channel, …”

There has been no mention of a channel up till now? Which channel is meant here?

Line 400: “where 𝒖 is the velocity (which has only a constant 𝑥-component), …”

One has to mention beforehand that flow is not being solved for here but assumed to be given as this or that. Some explanation should be given as to why the velocity can be none-zero at the top and bottom walls.

An improved description for test case 3.3 has been added to avoid confusions on the used (channel-like) model setup and the considered velocities (L 422-425). Note that as mentioned in the manuscript all assumptions regarding the model setup of this test case are taken from the reference study of Cirpka and Valocchi (2007) in order to ensure compatibility of our results with the analytically derived reference. 

-    Figure 5:

The axis labelled "Location (m)" should be labelled "y (m)" or similar.

The axis labels in Figure 5 have been updated as suggested.

-    Lines 436—439: “An important note to make here is with a relatively high influx, advection transport acts as the bottleneck for numerical time steps. Hence, reactions are performed at a quite slower pace (i.e. larger time steps). This separation of processes helps improve the overall run-time of the simulations.”

How are these different time-step sizes chosen?

Information on time-step size is added on line 468-472.

-    Lines 471—473: “Since microbial growth takes place at much larger time scales than the pore-scale transport processes no significant growth takes place during the simulated time period.”

Which time scales are we talking about here? How long does the system take to achieve steady state? How is biomass treated in terms of the distinction between sessile and planktonic cells? Are planktonic cells allowed to be transported in the system? Can bacterial cells also attach to the gas-water interface?

Why do the authors not have planktonic cells in the water phase a priori?

In the presented demonstration example we did not consider planktonic cells which is now clarified on lines 472-473. This was done for the sake of simplicity since the focus of this example is on the reactive transport of the dissolved species. In general, the chemotactic movement and attachment of bacterial cells to the gas-water interface are not considered in our model so far, but the advective transport of planktonic bacterial cells is considered. The criteria for transition from transient to steady state was determined by changes of water content becoming negligible (the exact value of the total simulated time to reach steady state is not available as it was not recorded but roughly it was in order of days).

-    Figure 6a:

The term "water saturation" is misleading. Be that as it may, it is really not clear how there can be a distribution of the water volume fraction as shown. I would expect there to be either water or air. At the interface between the two, one can have alpha values between zero and one, but this interface is miniscule in size, on the order of magnitude of Angstroms. The micro-CT discretisation used here is 6 micrometres. That means the air-water interface should be contained within one layer of computational cells. This doesn't seem to be the case here. Instead, there seems to be a distribution of the volume fraction within each bulk phase. This doesn't make sense to me.

The caption of Figure 6 has been updated. The 2D cutting plane is somewhat arbitrary in the sense that it might (and it can) cut through the 3D domain where the fluid-fluid interface exists. Also, in reality, the interface indeed is small but in simulating it using VOF, the interface can span several computational cells that in turn surely causes diffusion of the numerical errors. One main reason of all smoothing/sharpening steps of the interface explained beforehand is to reduce such effects.

-    Figure 6b:

Are these simply the initial conditions?

We clarified that shown biomass concentrations are nearly reflecting the initial conditions (line 507).

-    Lines 508:

Of the five videos available here, I was only able to play the first one. And it was of poor quality. 

We double checked the referenced videos and all were working fine. Also the originals that are submitted on the github repo are working fine.

-    Lines 523—526: “ Our modelling framework is properly designed for simulating biogeochemical processes such as carbon nitrogen-sulfur-phosphorus cycles in soil as well as mixing and migration of contaminants in both unsaturated soil and water aquifers.”

If the model cannot handle biomass transport processes, then it is not properly designed to simulate these processes.

As mentioned above, the chemotactic movement and attachment of bacterial cells to the gas-water interface are not considered in our model so far, but the advective transport of planktonic bacterial cells is considered. We think that this makes the model a suitable tool for number of systems and the processes therein. We outline a number of features not included in the model so far and that the modular structure of the model provides the opportunity to include them at a later stage (lines 573-582).

-    Line 526: “It comes with the benefit of explicit recognition of the soil structure, …”

What is an "explicit recognition of soil structure"?

A better explanation on this term has been added on lines 562-563.