Comment: In the paper “Supporting hierarchical soil biogeochemical modeling: Version 2 of the Biogeochemical Transport and Reaction model (BeTR-v2)””, authors Tang, Riley, and Zhu develop a new version of BeTR—a model development framework enabling investigation of different levels of complexity, process representations, and numerical method implementations. The new version (BeTR-v2) implements new numerical algorithms, is said to be a more efficient software, and can be run independent from host models. To demonstrate, the soil biogeochemistry model of ELMv1-ECA was implemented in BeTR-v2. The numerical solution was compared to analytical solutions, and simulations were performed at multiple scales (single layer, 1D, and global), comparing ELM versus ELMBeTR model configurations. Global simulations were benchmarked using ILAMB. Overall, this is a nice paper presenting tools and concepts that will be of interest to the GMD readership.

Response: We appreciate the reviewer’s positive assessment of our work. Based on the raised concerns, we have made changes accordingly in the submitted revision.

Comment: The most interesting result was that the different numerical implementation of ELMv1-ECA in BeTR-v2 led to substantially different predictions. Re-calibrating some of the key parameters was not sufficient to reduce these differences. I agree with the authors that this is an important source of uncertainty that is often not considered in the biogeochemical modeling community. I think this point could be better highlighted in the paper, though. For example, numerical methods as a source of uncertainty is mentioned in the middle of a list in the introduction, but not much is done to highlight or elaborate on this point (even though it becomes one of the main take-home points in the results and conclusion).

Response: Following the reviewer’s suggestions, we added at the end of the abstract “We contend that earth system models should strive to minimize this uncertainty by applying better numerical solvers.” In the introduction, when discussing this type of numerical uncertainty, we now write “In particular, when the differential equations of a model are approximated with inappropriate numerical solvers, the model may obtain answers that better match observations for wrong reasons because calibration may inappropriately make up for deficiencies in the model’s governing equations (i.e., type-I error that gets right answers with poor model formulations). This problem can result in incorrect inference of causality and interactions between processes. For instance, Tang et al. (2015) found that the simulated evapotranspiration agreed better with observations when the coupled equations for soil and root water exchange were purposely solved incorrectly in a sequential manner than when they were solved correctly as tightly coupled. Alternatively, if calibration cannot make up the deficiency caused by the inappropriate numerical method, one may assert that a right model formulation is wrong (i.e., type-II error that gets wrong answers with good model formulations). For example, when the 1D diffusion equation is solved with central difference in both time and space, the numerical solution actually approximates a wave equation instead, and this deficiency cannot be fixed by calibration. Both types of inference error will contribute to the uncertainty of climate-biogeochemistry feedback simulated by ESMs.”

Comment: I also think several aspects of section 2 could be clarified. This is a model description paper, so developments in this new version should be at least briefly described (even if they’ve been detailed elsewhere). For example, how is the code “more efficient” and could more details be provided regarding the new numerical algorithms? Please briefly describe the “multiple-flux-co-limiting solver" at first mention. I also suggest to clarify Table 1. This is important for keeping model configurations straight, but was not intuitive. For example, is the best description of the BeTR configurations really “Based on src/Applications/soil-farm/v1eca”. The reader is left to mine the text to understand what this means. Also, the BeTR configurations used the multiple-flux-co-limiting solver for belowground processes too, right?

Respond: We revised the text by adding brief introductions of new algorithms adopted in BeTR-v2, and also clarified that the BeTR configurations used the multiple-flux-co-limiting solver for belowground processes by default. Specifically, we now write in the last paragraph of section 2.1 “Gaseous and aqueous diffusion are solved together using the dual-phase algorithm (that assumes equilibrium between gaseous and aqueous phases) with the implicit time stepping method (Tang and Riley, 2014), which is equally accurate but simpler than the treatment in BeTR-v1 that requires calculating locations of wetting fronts in the soil. Solid phase diffusion is also solved implicitly. Aqueous advection is solved using the mass-conserving semi-Lagrangian approach (Manson and Wallis, 2000), which is more accurate (by reducing numerical dispersion) than the upstream scheme used in BeTR-v1. Biogeochemical reactions are solved using the multiple-flux-co-limiting algorithm (Tang and Riley, 2016), which considers the production and consumption fluxes concurrently, so that there is no delay between nutrient mineralization and its competition by consumption fluxes within a time step, a critical feature to resolve the nutrient limitation dynamics (Tang and Riley, 2018). To ensure numerical accuracy, within each modeling time step of ELM (which is 30 minute), each solver uses the adaptive time stepping that exits when either the relative difference between solutions based on coarse time step and halved time step is less than 0.1% or when the minimum time step (30 seconds) is reached.”

            We also revised Table 1 based on both reviewers’ suggestions:

Table 1. Summary of the configurations for the four global simulations.

Minor comments:

Comment: P1 L14-16: This sentence is kind of vague with regard to the algorithms and structural improvements.

Respond: Now the sentence is revised as “Here, we describe the new version, BeTR-v2, which adopts more robust numerical solvers for multiphase diffusion and advection, and coupling between biogeochemical reactions, and improves code modularization over BeTR-v1”.

Comment: P3 L19: change “share” to “sharing”

Response: Done.

Comment: P7 L1-3: The structure of this sentence is confusing. The analytical solutions should follow the colon, not equation 3.

Response: We revised the sentence as “As was done for BeTR-v1, two analytical solutions with different boundary conditions are employed to benchmark the numerical accuracy of the BeTR-v2 reactive transport solver that solves the following equation:”.

Comment: P9 L15: Elaborate. What previous findings?

Response: We adjusted the sentence as “When using the same parameters, the reordering required by ecacnpcaused significant differences in the simulated carbon and nutrient cycling compared to ELMv1-ECA (which confirms our previous findings in Tang and Riley (2018)), and such differences were inferred not correctable by calibration.” So that it clearly says that our previous finding is “reordering causes significant differences in model simulated carbon and nutrient cycling”.

Comment: P10 L18-19: “Comparing ELMv1-ECA and ELMv1-ECA”?

Response: We corrected it with “By comparing ELMv1-ECA-V and ELMv1-ECA”.

Comment: Figure 1: X-axes look like the title of the panel below.

Response: We looked for other approaches, like enlarging the space among panels, or putting the labels into the panel, neither of these methods make the illustration better. Therefore, we decided to stick with the current approach.

Comment: P14: Fig 2c is not referenced in the text.

Response: Now it is properly referenced in the sentence “Accordingly, the soil CO₂ concentration builds up continuously, with a seasonal cycle that has its maximum in July and minimum in March (Figure 2c).”.

Comment: Figure 3: Maybe make ELMv1-ECA a thicker line so that the reader can easily tell which model it is hidden behind.

Response: We made ELMv1-ECA a thicker line and used dash-dot line for ELMv1-ECA-V. Now the figures are more readable.

Comment: Table 2: Similar issue as above. I think one of these should be ELMv1-BeTR-ECA.

Response: Sorry for this typo, now it is corrected by identifying ELMv1-BeTR-ECA on the right side.

Comment: This paper describes an updated version of the Biogeochemical Transport and Reaction Model (BeTR-v2) including updated algorithms for reactive transport and numerical coupling with vegetation and hydrological processes. Simulations are conducted with the standalone version of the model and compared to analytical benchmarks, and a version of the model coupled to the E3SM Land Model is used to conduct and evaluate global simulations with alternate numerical implementations of soil biogeochemistry and plant soil coupling. The simulation results compare well with analytical benchmarks. Coupled land model simulations resulted in different carbon, nitrogen, and phosphorus cycle outcomes for the different numerical implementations.

Overall, the manuscript is well written and provides a clear description of the model developments, the simulations that were performed, and the results. There are a few typographical errors and some areas where clarity could be improved.

Response: We appreciate the reviewer’s positive assessment of our work.  

Comment: Page 3, Line 3: I would word this “sharing of common process representations…”

Response: Done.

Comment: Page 3, Line 19-20: …enable efficient code and knowledge sharING … improvements that have BEEN brought…

Response: Done.

Comment: Page 4, line 10: SINCE significant code rewriting…

Response: This phrase is now replaced with “that since significant code re-writing and data restructuring have occurred after ELM branched out from CLM4.5”.

Comment: Equation 1: Cg is used in the third right-hand-side term (with Ds), and I think it should the Cs instead

Response: The typo is now corrected.

Comment: Page 5, line 18-19: A bit more description of the solver method would be helpful so readers can get a basic understanding without reading a different paper. Also, does the time stepping method account for truncation errors at longer time steps? Is some adaptive time stepping included for cases where the model time step is too long to resolve fast biogeochemical processes (maybe not important in the simulations presented here but potentially important in some applications such as explicit tracking of oxygen concentrations)?

Response: Please also see our response to a similar comment from Reviewer #1, above. For our description of the numerical methods, we revised the text as “Gaseous and aqueous diffusion are solved together using the dual-phase algorithm (that assumes equilibrium between gaseous and aqueous phases) with the implicit time stepping method (Tang and Riley, 2014), which is equally accurate but simpler than the treatment in BeTR-v1 that requires calculating locations of wetting fronts in the soil. Solid phase diffusion is also solved implicitly. Aqueous advection is solved using the mass-conserving semi-Lagrangian approach (Manson and Wallis, 2000), which is more accurate (by reducing numerical dispersion) than the upstream scheme used in BeTR-v1. Biogeochemical reactions are solved using the multiple-flux-co-limiting algorithm (Tang and Riley, 2016), which considers the production and consumption fluxes concurrently, so that there is no delay between nutrient mineralization and its competition by consumption fluxes within a time step, a critical feature to resolve the nutrient limitation dynamics (Tang and Riley, 2018). To ensure numerical accuracy, within each modeling time step of ELM (which is 30 minute), each solver uses the adaptive time stepping that exits when either the relative difference between solutions based on coarse time step and halved time step is less than 0.1% or when the minimum time step (30 seconds) is reached.”

Comment: Section 2.4 and Table 1: I had a hard time keeping track of what the differences were between the different simulations. The short descriptions in Table 1 are not very informative because they refer to specific code directories rather than numerical methods, and include several different contrasting numerical approaches described in only one table column. I would suggest adding more columns to the table to clearly differentiate the features of the different implementations. Separate columns could include plant-soil competition solver, plant allocation solver, and parameterization which all varied across different simulations. I would avoid referring to specific code directories where possible and instead refer to the differences in underlying methods, which is more universal. In the text description (page 10), the use of italicized “ecacnp” in some places and the names of the implementations (e.g., ELMv1-ECA) in others is confusing and seems specific to this code base rather than a general description of numerical approaches. I would suggest using only one terminology, or else including the “ecacnp” terminology in Table 1 so it’s easier to keep track of the different terms.

Response: By following these suggestions, and those from Reviewer #1, we revised the table by adding more entries to describe the differences among model configurations:

Table 1. Summary of the configurations for the four global simulations.

Comment: Page 10, line 18-19: “Comparing ELMv1-ECA and ELMv1-ECA” - these are both the same. Should one be different?

Response: It is corrected as “By comparing ELMv1-ECA-V and ELMv1-ECA”.

Comment: Figure 2: There was not an explanation of how column integrated heterotrophic respiration, soil surface CO2 flux, and CO2 infiltration rate were calculated and what exactly they represent. I assume the surface flux takes transport of gaseous and dissolved CO2 into account whereas integrated HR is instantaneous production?

Response: To address this issue, we revised the caption of Figure 2 as “(a) Heterotrophic CO₂ flux simulated by the 10-cm thick single layer model; (b) column integrated heterotrophic flux (by summing up contributions from all layers in the soil column), soil surface CO₂ flux (from capillary exchange and diffusion considering equilibrium between gaseous and aqueous phases) and CO₂ infiltration flux; (c) evolution of soil CO₂concentration corresponding to panel (b). ”

Comment: Figure 3: Why was accelerated spinup used here instead of the normal spinup or historical simulation?

Response: We made this choice because, for accelerated spinup, all models started with the same initial conditions. However, when exiting accelerated spinup, there will be a significant increase of the soil organic matter pools (as accelerated spinup means to shorten the time needed for the model to reach equilibrium) through rescaling (based on simulation dependent rescaling factors) of the soil organic matter pools, which further amplifying the model difference (see further discussions in Koven et al. (2013), where the accelerated spinup method was developed). Therefore, using accelerated spinup reduces the aliasing impact resulting from the rescaling factors.

Comment: Page 18, line 7-8 and Table 2: I would use the PFT names rather than numbers which are not meaningful to readers to are not closely familiar with this land model

Response: We now spelled explicitly the PFT names in the text as PFT-4 (broad leaf evergreen tropical tree) and PFT-6 (broadleaf deciduous tropical tree). We kept the use of PFT-4 and PFT-6 in the table entries (for formatting purpose as the use of full name will make the table too busy), but do note the PFT names in the table’s caption.

Comment: Page 18, line 11: “other variables” - Explain which variables

Response: We now spelled them explicitly as “(NBP, vegetation carbon, soil carbon, and total ecosystem carbon)”.

Comment: Table 2: Both columns have the same heading “ELMv1-ECA”. One should be ELMv1-BeTR-ECA.

Response: This was a typo brought in when we revised the initial submission before discussion. Now it is corrected by identifying “ELMv1-BeTR-ECA” on the right side.

Comment: Page 19, line 16-18: This sentence feels oddly judgmental. The previous sentence reports better agreement with some benchmarks but there also seems to be worse agreement with others, so it might be more balanced to say that numerical differences can significantly change model outcomes even without changing the underlying differential equations of a model.

Response: We revised the text to address this concern: “This better agreement between ELMv1-BeTR-ECA and some benchmarks suggests that the numerical difference can significantly influence the performance of a supposedly good mathematical representations of ecosystem biogeochemistry.”

Comment: Page 21, line 3: Estimates of P dynamics also changed, not just N.

Response: We revised the sentence as “Therefore, a model calibration using C cycle variables led to very different estimates of N cycle parameters and thereby different nitrogen dynamics and phosphorus dynamics (through N and P co-modulated biogeochemical feedbacks).”

Comment: Table 3: The units of the numbers are never described and it’s not clear whether a higher number means a better or worse fit to the benchmarks.

Response: In the caption of Table 3, we now added “All metrics are normalized to the range from 0 to 1, where greater values indicate better performance.”.

Comment: Page 23, line 12: What is meant specifically by “numerically more robust”?

Response: To clarify this issue, we revised the sentence as “We found that because the multiple-flux-co-limiting numerical solver more tightly couples plant and soil processes during nutrient competition (so that it is numerically more robust than the solver used by default ELMv1-ECA model; Tang and Riley (2018, 2016)).”

Comment: Page 24, line 2: I don’t think there is a basis here to decide whether model parameters are “incorrect” or that a particular numerical coupling is “inappropriate.” This study shows that different numerical approaches can yield different results. Without a clear demonstration that one approach or the other fails relative to some benchmarks I don’t think it can support a declaration that one is right or wrong. A more balanced wording might be that different numerical approaches can significantly change model behavior and that care should be taken to evaluate whether re-parameterization is necessary following numerical changes. This was clearly demonstrated in this study where the land model needed to be recalibrated following a change to the numerical coupling scheme.

Response: We modified the tone of the last sentence by using “inappropriate numerical coupling will potentially result in incorrect model parameters that may affect predictions of carbon cycling variables under a changing climate and increasing atmospheric CO₂ concentrations.” We maintain our opinion that a proper numerical approach is essential to obtain numerical solutions that are consistent with the differential equations being solved. In previous papers (Tang and Riley, 2016, 2018), we learned that the solution strategy used by the default ELMv1-ECA model delays the availability of newly mineralized inorganic nutrient for uptake to the next time step, while the multiple-flux-co-limiting solver synchronizes mineralization and uptake as formulated by the governing equations.Therefore, as time step size decreases, the default ELMv1-ECA model will not converge to its governing equations. If we assume that the governing equations are correct, then the default ELMv1-ECA has to make up for this deficiency by using inappropriate parameters. In turn, we may obtain the right answer for wrong model formulations because inappropriate numeric solutions and model calibration together may mask the insufficiency in model formulations, or assert that a correct model formulation is incorrect because the numerical code has to use inappropriate parameter values.