Reply on RC2

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.


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: 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 andRiley (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 2 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 andRiley, 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.