the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
P3DBRNS v1.0.0: a threedimensional, multiphase, multicomponent, porescale reactive transport modelling package for simulating biogeochemical processes in subsurface environments
Amir Golparvar
Matthias Kästner
Martin Thullner
Download
 Final revised paper (published on 01 Feb 2024)
 Supplement to the final revised paper
 Preprint (discussion started on 03 Jun 2022)
 Supplement to the preprint
Interactive discussion
Status: closed

RC1: 'Comment on gmd202286', Anonymous Referee #1, 08 Dec 2022
This paper describes a new opensource toolbox for simulating reactive biogeochemical processes at the pore scale. It couples OpenFoam, for interface tracking and fluid flow, with the Biogeochemical Reaction Network Simulator (BRNS), for complex reactions. An improved VoF solver is implemented based on the procedure proposed by Raeini et al. (2012). The code is validated with three test cases where analytical solutions are available. Finally, the capability of the model to represent realistic scenarios is tested on a 3D microstructure with coupled twophase flow and bacterial growth following a complex set of reactions.
The manuscript is well written and clearly and concisely documents the model and the implementation. Coupling twophase flow to complex reactions in a complicated pore space is a nontrivial task. However, based on the thoroughly presented benchmark tests it seems this code is able to do it in an efficient manner. Hence I believe that the P3DBRNS toolbox will be of great benefit to the community. Provided the authors can satisfactorily answer the minor remarks below, I will recommend publication in GMD.
Minor remarks:
 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 welldocumented 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?
 In at least one of the test cases 3.13.3, it would be useful with a discussion on how grid resolution affects the convergence to the analytical solution.
 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 accomodated by a reinterpretation of the pressure P. Would P3DBRNS be able to model e.g. droplet motion due to a concentration gradient and solubility difference in the two phases?
 section 3.1: How was the mesh generated? Is it skewed or staircaselike at the diagonal side? How is the sharp top corner handled? Some more information here would be useful.Typos/formatting:
 Caption of Fig. 1 also inserted at line 115.
 line 359: incomplete sentence
 line 412: simluations > simulations
 line 487: growth > growCitation: https://doi.org/10.5194/gmd202286RC1 
AC1: 'Reply on RC1', Amir Golparvar, 31 Aug 2023
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 welldocumented 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 309313.
 In at least one of the test cases 3.13.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 406410).
 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 reinterpretation of the pressure P. Would P3DBRNS 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 581583.
 section 3.1: How was the mesh generated? Is it skewed or staircaselike 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).
Citation: https://doi.org/10.5194/gmd202286AC1

AC1: 'Reply on RC1', Amir Golparvar, 31 Aug 2023

RC2: 'Comment on gmd202286', Anonymous Referee #2, 23 Jun 2023
I was surprised to see how long this manuscript has been under review. The authors do deserve a speedy and fair review. I will do my best to deliver both.
I will say right away that I am sceptical about two things in this manuscript:
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.
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.
Here are questions and comments to specific parts of the manuscript:
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 nonlinear. How biomass concentration and bulknutrient 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.
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 domainspecific terms like "pore scale" and "Darcy scale" or define how they use scalerelated terms at the beginning of the manuscript.
Section 1:
The authors' excessive use of parentheses diminishes legibility. I advise against structuring sentences around parentheses. In one particular case on Lines 7679, 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.
Lines 8889:
Why not the Stokes equations for Re<<1? In any case, the NS equations should be plural.
Eqs. (1) and (2):
This is the incompressible form of the NS equations. Do the authors really assume that the gas phase is incompressible?
Eq. (5):
How does one obtain the individual velocities of the fluids? How does one obtain the velocity of the interface?
Line 141:
How are u_alpha und u_beta different from u_i?
Lind 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.
Lines 145146: “… which computationally helps with maintaining stiffness of the interface.”
It really is not clear what the stiffness of the interface is.
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?
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?
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 …?
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.
Eq. (19):
u_n should be defined.
Eq. (21):
How is this interface normal vector used in the simulations?
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?
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.
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.
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 nonezero at the top and bottom walls.
Figure 5:
The axis labelled "Location (m)" should be labelled "y (m)" or similar.
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 runtime of the simulations.”
How are these different timestep sizes chosen?
Lines 471—473: “Since microbial growth takes place at much larger time scales than the porescale 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 gaswater interface?
Why do the authors not have planktonic cells in the water phase a priori?
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 microCT discretisation used here is 6 micrometres. That means the airwater 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.
Figure 6b:
Are these simply the initial conditions?
Lines 508:
Of the five videos available here, I was only able to play the first one. And it was of poor quality.
Lines 523—526: “ Our modelling framework is properly designed for simulating biogeochemical processes such as carbon nitrogensulfurphosphorus 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.
Line 526: “It comes with the benefit of explicit recognition of the soil structure, …”
What is an "explicit recognition of soil structure"?
Citation: https://doi.org/10.5194/gmd202286RC2 
AC2: 'Reply on RC2', Amir Golparvar, 31 Aug 2023
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 259263). 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 indepth 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 nonlinear. How biomass concentration and bulknutrient 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 5961).
 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 domainspecific terms like "pore scale" and "Darcy scale" or define how they use scalerelated 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 7679, 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 8083). Throughout the revised manuscript, the use of parentheses has been reduced and the use of nested parentheses has been avoided.
 Lines 8889:
Why not the Stokes equations for Re<<1? In any case, the NS equations should be plural.The NavierStokes equations are chosen for a couple of reasons: 1) the velocities we are dealing with satisfy very low Mach numbers, and 2) the NavierStokes 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 117118). 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 139140).
 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 128130.
 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 145146: “… 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 148150.
 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 154155.
 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 156158.
 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 164165 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 245246.
 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 257258 , and the meaning of smoothing interface is explained before on line 155156).
 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 259263. 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 401405.
 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 nonezero at the top and bottom walls.An improved description for test case 3.3 has been added to avoid confusions on the used (channellike) model setup and the considered velocities (L 422425). 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 runtime of the simulations.”
How are these different timestep sizes chosen?Information on timestep size is added on line 468472.
 Lines 471—473: “Since microbial growth takes place at much larger time scales than the porescale 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 gaswater 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 472473. 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 gaswater 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 microCT discretisation used here is 6 micrometres. That means the airwater 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 fluidfluid 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 nitrogensulfurphosphorus 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 gaswater 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 573582).
 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 562563.
Citation: https://doi.org/10.5194/gmd202286AC2

AC2: 'Reply on RC2', Amir Golparvar, 31 Aug 2023