Supporting hierarchical soil biogeochemical modeling: version 2 of the Biogeochemical Transport and Reaction model (BeTR-v2)

Reliable soil biogeochemical modeling is a prerequisite for credible projections of climate change and associated ecosystem feedbacks. This recognition has called for frameworks that can support flexible and efficient development and application of new or alternative soil biogeochemical modules in Earth system models (ESMs). The the Biogeochemical Transport and Reaction model version 1 (BeTR-v1) code (i.e., CLM4-BeTR) is one such framework designed to accelerate the development and integration of new soil biogeochemistry formulations into ESMs and to analyze structural uncertainty in ESM simulations. With a generic reactive transport capability, BeTR-v1 can represent multiphase (e.g., gaseous, aqueous, and solid), multi-tracer (e.g., nitrate and organic carbon), and multi-organism (e.g., plants, bacteria, and fungi) dynamics. Here, we describe the new version, Biogeochemical Transport and Reaction model version 2 (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. BeTR-v2 better supports different mathematical formulations in a hierarchical manner by allowing the resultant model be run for a single topsoil layer or a vertically resolved soil column, and it allows the model to be fully coupled with the land component of the Energy Exascale Earth System Model (E3SM). We demonstrate the capability of BeTR-v2 with benchmark cases and example soil biogeochemical (BGC) implementations. By taking advantage of BeTR-v2's generic structure integrated in E3SM, we then found that calibration could not resolve biases introduced by different numerical coupling strategies of plant–soil biogeochemistry. These results highlight the importance of numerically robust implementation of soil biogeochemistry and coupling with hydrology, thermal dynamics, and plants – capabilities that the open-source BeTR-v2 provides. We contend that Earth system models should strive to minimize this uncertainty by applying better numerical solvers.

1 Introduction

Soil biogeochemical (BGC) modeling is essential for predicting terrestrial ecosystem dynamics in natural and managed lands and is an important component of Earth system models (ESMs) that perform projections and hindcasts of ecosystem–climate feedbacks (e.g., Golaz et al., 2019, Hurrell et al., 2013). Over the years, many soil BGC models have been developed, each of which has its merits and deficiencies (e.g., Zaehle et al., 2014). Nonetheless, when the carbon dynamics simulated by many different ecosystem models (including those integrated within ESMs) are synthesized, large uncertainties prevail (e.g., Davies-Barnard et al., 2020; Todd-Brown et al., 2014). Several sources, including parametric, structural, numerical, and boundary conditions (e.g., forcing and initial conditions), have been cited to explain the large modeling uncertainty (Ahlstrom et al., 2013; Huntzinger et al., 2017; Tang and Riley, 2018). Accordingly, various strategies have been proposed to alleviate or constrain these types of uncertainty (Luo et al., 2012; Riley et al., 2021; Wieder et al., 2015; Zhu et al., 2019).

From the many uncertainty quantification studies, one urgent need identified is for a software platform to support systematic comparisons of how different soil BGC formulations influence simulated soil BGC feedbacks within a single ecosystem model, thereby facilitating more process-specific hypotheses testing. To address this infrastructural gap, the Biogeochemical Transport and Reaction model version 1 (BeTR-v1) was developed as a submodule of the Community Land Model version 4 (CLM4; Tang et al., 2013). The design philosophy of BeTR is that, assuming all other factors are equal, one can investigate (1) the various levels of mechanistic complexity and model formulations even for the same number of biogeochemical processes being represented, (2) the same mathematical formulation implemented with different numerical methods (e.g., Tang and Riley, 2018), and (3) the sharing of common process representations across diverse ecosystem types. For point (1), some researchers have argued that more explicit representation of processes are needed to resolve complex soil BGC dynamics (Berardi et al., 2020; Riley et al., 2021; Tang and Riley, 2020; Wieder et al., 2015). For instance, active transport of various substrates in soils is needed if microbial dynamics are to be explicitly considered (e.g., Ahrens et al., 2015; Dwivedi et al., 2017; Grant, 2013; Riley et al., 2021). Point (2) is less discussed in land BGC modeling but has been clearly called out in computational biology and computational geophysical fluids (e.g., Gross et al., 2018; Petrovskii and Petrovskaya, 2012). Specifically, 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., a type I error that gets the 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., a type II error that gets the wrong answers with good model formulations). For example, when the 1-D 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. For point (3), aqueous environments, such as wetlands, groundwater, rivers, lakes, and oceans, share many biogeochemical processes with soils, but different representations are often used (e.g., early diagenesis in marine sediments vs. terrestrial soil organic matter decomposition; Koven et al., 2013; Munhoven, 2021). Therefore, developing a generic modeling infrastructure could benefit the mechanistic consistency in integrated Earth system BGC modeling (e.g., Fisher and Koven, 2020).

Since the publication of BeTR-v1, we made a number of new developments in numerical algorithms for transport and BGC process coupling (Tang and Riley, 2016, 2014). We also made better use of the object-oriented programming feature of Fortran 2003 to enable more efficient code sharing among the Fortran modules. All of these are integrated into the Biogeochemical Transport and Reaction model version 2 (BeTR-v2): a platform to accelerate soil BGC model development and enable efficient code and knowledge sharing among ESM BGC components. Below, we describe the improvements that have brought into BeTR-v2 since BeTR-v1 in detail, give some examples based on its integration with the land module of the Energy Exascale Earth System Model (ELM; Burrows et al., 2020), and apply the model to evaluate whether parameter calibration can resolve the uncertainty associated with different numerical implementations.

2 Model developments

BeTR-v2 adopts the same governing equations as in BeTR-v1 but implements improved numerical algorithms (developed after the development of BeTR-v1) for both tracer transport and biogeochemical reaction coupling (to be described in Sect. 2.1). Compared to BeTR-v1, BeTR-v2 is also better modularized and structured software-wise. Importantly, BeTR-v2 is a standalone capability to support more efficient model development and soil BGC coupling to alternative host models (E3SM in this case), while BeTR-v1 is embedded within CLM4.5 and cannot be run independently. It is noted that since significant code rewriting and data restructuring have occurred after ELM branched out from CLM4.5, a comparison between BeTR-v1 and BeTR-v2 is not feasible.

2.1 Governing equations and numerical implementation

The same as BeTR-v1, BeTR-v2 solves the following one-dimensional reactive-transport equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\theta }}_{\mathrm{w}}{C}_{\mathrm{w}}\right.& +\left.{\mathit{\epsilon }}_g{C}_{\mathrm{g}}+\left(\mathrm{1}-{\mathit{\theta }}_{\mathrm{w}}-{\mathit{\epsilon }}_g\right)C_s\right)\\ & ={\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\theta }}_{\mathrm{w}}{\mathit{\tau }}_{\mathrm{w}}{D}_{\mathrm{w}}{\displaystyle \frac{\partial C_{\mathrm{w}}}{\partial z}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\epsilon }}_g{\mathit{\tau }}_{\mathrm{g}}{D}_{\mathrm{g}}{\displaystyle \frac{\partial C_{\mathrm{g}}}{\partial z}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\left(\mathrm{1}-{\mathit{\theta }}_{\mathrm{w}}-{\mathit{\epsilon }}_g\right)D_{\mathrm{s}}{\displaystyle \frac{\partial C_s}{\partial z}}\right)\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(q_{\mathrm{w}}{\displaystyle \frac{\partial C_{\mathrm{w}}}{\partial z}}\right)-E_{\mathrm{bub}}-T_r+R_{\mathrm{bgc}},\end{array}\end{array}$$

where θ_w is volumetric water content (m³ m⁻³), C_w is aqueous concentration (mol m⁻³), ε_g is air-filled porosity (m³ m⁻³), C_g is gaseous concentration (mol m⁻³), C_s is adsorbed concentration (mol m⁻³), τ_w is aqueous tortuosity (unitless), τ_g is gaseous tortuosity (unitless), D_w is aqueous diffusivity (m² s⁻¹), and D_g is gaseous diffusivity (m² s⁻¹). D_s is solid-phase diffusivity (m² s⁻¹), as an approximation for processes like bioturbation (Jarvis et al., 2010) and cryoturbation (Koven et al., 2009), which do not have mechanistically based formulations. q_w is aqueous advection (m s⁻¹), E_bub is ebullition flux (mol m⁻³ s⁻¹), T_r is transport via transpiration flux (mol m⁻³ s⁻¹; including aerenchyma-enabled gas transport), and R_bgc is biogeochemical reaction rate (mol m⁻³ s⁻¹).

Depending on how z (m) is interpreted, the transport can be regarded as either vertical or horizontal, provided that the advection velocity q_w is defined accordingly.

Besides the one-dimensional solution, BeTR-v2 also solves standalone single-layer models (a new feature that BeTR-v1 does not have), whose governing equation is as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\theta }}_{\mathrm{w}}{C}_{\mathrm{w}}\right.& +\left.{\mathit{\epsilon }}_g{C}_{\mathrm{g}}+\left(\mathrm{1}-{\mathit{\theta }}_{\mathrm{w}}-{\mathit{\epsilon }}_g\right)C_s\right)\\ & =R_{\mathrm{bgc}}+{\displaystyle \frac{g_{\mathrm{as}}}{\mathrm{\Delta }z}}\left(C_{\mathrm{atm}}-C_{\mathrm{g}}\right),\end{array}\end{array}$$

where g_as is the conductance (m s⁻¹) for the gaseous exchange of the tracer between soil and air, C_atm is the atmospheric gas concentration (mol m⁻³), and Δz (m) is the single-layer thickness (taken as 10 cm in the examples below). Simulations from Eq. (2) can be used for comparisons with incubation experiments.

To achieve numerically robust solutions and to adapt to the hydrological dynamics in ELM (which solves surface runoff, variably saturated flow, and subsurface drainage sequentially), Eq. (1) is solved with the operator splitting approach (Simpson and Landman, 2007). Specifically, tracer loss through surface runoff is solved first, then biogeochemical reactions are solved, followed by advection, diffusion, ebullition (using the hydrostatic approximation from Tang et al., 2010), and subsurface drainage. 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, which 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 min) each solver uses the adaptive time stepping that exits when either the relative difference between solutions based on the coarse time step and halved time step is less than 0.1 % or when the minimum time step (30 s) is reached.

2.2 Code structure

The open-source code BeTR-v2 can be found at https://github.com/BeTR-biogeochemistry-modeling/sbetr/tree/jinyuntang/ELM-BeTR-ECA.r (last access: 21 February 2022) (and it is also archived at https://doi.org/10.5281/zenodo.5526854), where each folder has a readme.md file explaining the purpose of each sub-directory or each source code file. For interested users, the major code to look into is under the directories “sbetr/src/betr” (core betr code), “sbetr/src/Applications” (customized BGC modules), “sbetr/src/driver/” (driver files for 1-D models and APIs to couple BeTR with land models like ELM), and “sbetr/src/jarmodel/” (driver files for single-layer models). In the directory “src/Applications/soil-farm/”, each folder (except “bgcfarm_util”) is a soil BGC implementation. For this paper, the BGC implementation ELMv1-BeTR-ECA under the directory “sbetr/src/Applications/soil-farm/v1eca/” is used to test the global simulation capability of BeTR, while “ecacnp” under “sbetr/src/Applications/soil-farm/ecacnp/” is used to demonstrate the standalone simulation capability. More details about why two implementations are used for different purposes will be explained in Sect. 2.4. Instructions on how to build and run BeTR on typical platforms can be found at https://github.com/BeTR-biogeochemistry-modeling/sbetr (last access: 21 February 2022).

2.3 Model benchmark and verification

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:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D{\displaystyle \frac{\partial C}{\partial z}}\right)-u{\displaystyle \frac{\partial C}{\partial z}},\end{array}$$

where D (m² s⁻¹) and u (m s⁻¹) are diffusivity and advection velocity, respectively.

The two analytical solutions are as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}C& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{erfc}\left({\displaystyle \frac{z-ut}{\mathrm{2}\sqrt{\mathrm{D}t}}}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\exp \left({\displaystyle \frac{uz}{D}}\right)\mathrm{erfc}\left({\displaystyle \frac{z+ut}{\mathrm{2}\sqrt{\mathrm{D}t}}}\right)\\ & +\left[\mathrm{1}+{\displaystyle \frac{u}{\mathrm{2}D}}\left(\mathrm{2}L-z+ut\right)\right]\\ & \times \exp \left({\displaystyle \frac{uL}{D}}\right)\mathrm{erfc}\left({\displaystyle \frac{\mathrm{2}L-z+ut}{\mathrm{2}\sqrt{\mathrm{D}t}}}\right)\\ & -\sqrt{{\displaystyle \frac{u^{\mathrm{2}}t}{\mathit{\pi }D}}}\exp \left[{\displaystyle \frac{uL}{D}}-{\displaystyle \frac{{\left(\mathrm{2}L-z+ut\right)}^{\mathrm{2}}}{\mathrm{4}\mathrm{D}t}}\right],\end{array}\end{array}$$

for a pulse tracer imposed at the top of a 1-D soil column of length L (Kumar et al., 2009), where erfc(x) is the complementary error function of x, and

$$\begin{array}{}\text{(5)}& \begin{array}{rl}C& =C_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{\mathrm{2}}{A}_i\exp \left(-{\displaystyle \frac{u}{\mathrm{2}D}}-{\displaystyle \frac{\sqrt{\mathrm{2}}z}{\mathrm{4}D}}\sqrt{u^{\mathrm{2}}+\sqrt{u^{\mathrm{4}}+\mathrm{16}{D}^{\mathrm{2}}{\mathit{\omega }}_i^{\mathrm{2}}}}\right)\\ & \times \sin \left({\mathit{\omega }}_{i}t-{\displaystyle \frac{\sqrt{\mathrm{2}}{\mathit{\omega }}_{i}z}{\sqrt{u^{\mathrm{2}}+\sqrt{u^{\mathrm{4}}+\mathrm{16}{D}^{\mathrm{2}}{\mathit{\omega }}_i^{\mathrm{2}}}}}}\right),\end{array}\end{array}$$

for a wave-like tracer imposed at the top of a 1-D soil column, where $C_{\mathrm{0}}=\mathrm{12}/\mathrm{23}$ mol of tracer per cubic meter, $A_{\mathrm{1}}=\mathrm{9}/\mathrm{23}$ mol of tracer per cubic meter, $A_{\mathrm{2}}=\mathrm{2}/\mathrm{23}$ mol of tracer per cubic meter, ${\mathit{\omega }}_{\mathrm{1}}=\mathrm{2}\mathit{\pi }/\mathrm{365}$ d⁻¹, and ω₂=2π d⁻¹. For both cases, $u={\mathrm{10}}^{-\mathrm{7}}$ m s⁻¹, $D=\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹, and L=35.17 m, and the lower boundary condition is as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial C}{\partial t}}+u{\displaystyle \frac{\partial C}{\partial z}}=\mathrm{0}.\end{array}$$

These two tests have been integrated as part of the standard test suite for BeTR-v2. After doing “make test” by following the instructions, their output can be found under “regression-tests/tests/standalone/” as analytical-adr-b1 and analytical-adr-b2, respectively.

2.4 Example soil BGC implementation

We re-implemented the soil BGC from ELMv1-ECA (Zhu et al., 2019) to illustrate the capability of BeTR-v2. By re-implementation we mean to solve the ELMv1-ECA soil BGC in the form of Eq. (1) in BeTR-v2 based on input fluxes of organic matter and nutrients from ELM and return fluxes of carbon and nutrient fluxes to ELM (after calculations in BeTR-v2). The ELMv1-ECA soil BGC model consists of a century-like soil carbon decomposition cascade based on Koven et al. (2013), but competition of inorganic nitrogen (ammonium and nitrate) and phosphorus between plants and (the implicitly represented) microbial organisms is calculated using equilibrium chemistry approximation kinetics (Tang and Riley, 2013; Zhu et al., 2016). Differing from other nutrient competition paradigms, particularly the popular relative demand approach that is implemented in many existing BGC models (Riley et al., 2018), the ECA approach explicitly considers the kinetic traits of the involved plants, microbes, and soil sorption surfaces and is more capable of resolving the diverse dynamical processes involved in plant–soil nutrient competition (Medvigy et al., 2019; Tang and Riley, 2021; Zhu et al., 2017). ELMv1-ECA also allows plant organs to vary their carbon to nutrient ratios to fluctuate within empirically inferred ranges. The performance of ELMv1-ECA was evaluated comprehensively by using the International Land Model Benchmarking (ILAMB) system (Collier et al., 2018; Zhu et al., 2019).

In BeTR-v2, a soil BGC formulation is implemented in two-steps: (1) a single-layer implementation of the given soil BGC process, which is then populated through BeTR-v2's reactive transport code to form (2) a 1-D vertically resolved implementation of the soil BGC process. Take ecacnp for example, the single-layer implementation is in the folder “ecacnp/ecacnp1layer/”, which is extended to 1-D soil column using codes in the folder “ecacnp/ecacnpNlayer/”, and corresponding model parameters are set using codes in the folder “ecacnp/ecacnpPara/”. The single-layer code is called by the driver file in “sbetr/src/jarmodel/driver” to conduct single-layer simulations. In contrast, BeTR-v1 did not separate the single-layer models from the vertically resolved models. The soil BGC formulations in the folder “src/Applications/soil-farm/ecacnp” adopt the mathematical formulation of the soil BGC module from ELMv1-ECA but concurrently solve the incoming litter input, soil carbon decomposition, nitrification–denitrification, plant–soil nutrient competition, etc., when updating the soil BGC state variables. The treatment in ecacnp enables us to cleanly separate vegetation and soil BGC in ELM but requires a significant amount of reordering of the subroutines and two additional parameters for plant BGC; i.e., besides the introduction of two parameters on the turnover of plant nitrogen and phosphorus storage pools, ecacnp has to compute soil BGC and plant–soil nutrient competition after all the phenological updates (including plant carbon and nutrient allocation for growth and vegetation mortality induced by fire and other natural factors), whereas ELMv1-ECA computes plant–soil nutrient competition before updating the phenological variables. When using the same parameters, the reordering required by ecacnp caused 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 to be not correctable by calibration. Therefore, we did another implementation of the soil BGC of ELMv1-ECA in BeTR-v2 by following the exact subroutine execution order of ELMv1-ECA and named it ECAv1-BeTR-ECA0 (the corresponding code is in src/Applications/soil-farm/v1eca). ELM-BeTR-ECA0 can run with the same parameters as those in ELMv1-ECA, thus minimizing the need to change parameters to obtain acceptable simulations. Since we found that ELM-BeTR-ECA0 performed significantly differently to ELMv1-ECA, we include ELMv1-BeTR-ECA by calibrating ELM-BeTR-ECA0 against ELMv1-ECA to evaluate to what extent the structural difference can be reduced by adjusting model parameters.

Overall, ELMv1-BeTR-ECA0 differs from ELMv1-ECA in the following two aspects: (1) in belowground BGC, ELMv1-BeTR-ECA0 solves the plant–soil nutrient competition using the multiple-flux co-limiting solver such that nutrient mineralization and immobilization are concurrently coupled (Tang and Riley, 2016), whereas ELMv1-ECA solves nutrient mineralization and uptake asynchronously using the explicit Euler scheme, where newly mineralized nutrients at the current time step become available only in the next time step. (2) For aboveground BGC, ELMv1-BeTR-ECA0 solves plant carbon and nutrient allocation for growth using the multiple-flux co-limiting solver to avoid sporadic negative litter flux inputs to the soil BGC, whereas ELMv1-ECA updates these state variable sequentially using the explicit Euler scheme and only corrects negative variables (which occur sporadically) at the end of each simulation step. ELMv1-BeTR-ECA differs from ELMv1-BeTR-ECA0 with values of the model parameters. Besides these three configurations, we also included ELMv1-ECA-V that differs from ELMv1-ECA by using the multiple-flux co-limiting solver for plant carbon and nutrient allocation. A brief description of the four model configurations is given in Table 1.

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

All in all, by using ecacnp soil BGC for standalone simulations (which can be created by following the instructions in the front page of the sbetr Git repository), we show how BeTR-v2 can support an easy extension from a single-layer model to a 1-D vertically resolved model. By comparing ELMv1-ECA-V and ELMv1-ECA, we verify the correct implementation of the multiple-flux co-limiting solver. Using ELMv1-BeTR-ECA0 for ELM-coupled simulations, we demonstrate that different numerical implementations for the same BGC formulation results in significantly different model behavior. More importantly, by using ELMv1-BeTR-ECA we then show that such model behavior differences cannot be resolved by recalibrating model parameters.

We drove the standalone 1-D model simulation with 365 d half-hourly time step output from an earlier ELM simulation Zhu et al. (2019). Model input variables include vertically resolved litter fluxes (of carbon, nitrogen, and phosphorus), soil temperature and moisture, infiltration flux, and atmospheric pressure and temperature. The example simulations are from the CA-Gro Ameriflux site: Groundhog River, Ontario, with boreal mixed-wood forest (see https://ameriflux.lbl.gov/sites/siteinfo/CA-Gro for more descriptions of the site, last access: 21 February 2022). The single-layer simulation used inputs for the first soil layer. For all global simulations, we used 200 year accelerated spinup, 600 year regular spinup (Koven et al., 2013; Zhu et al., 2019), and a transient simulation from 1850–2010. All forcing data are the same from Zhu et al. (2019), and all simulations are 1.9 latitude by 2.5^∘ longitude spatial resolution.

3 Results

3.1 Benchmark of the BeTR-v2 transport code

We found that the BeTR-v2 advection–diffusion solver solutions agreed well with the analytical solutions for both benchmarking cases (Fig. 1). For the pulse tracer solution (left column in Fig. 1), the depth-averaged root-mean-square error (RMSE) and R² of linear fitting are 0.013 mol m⁻³ and 0.999, respectively, and for the wave-like solution (right column in Fig. 1), these two metrics are 0.0050 mol m⁻³ and 0.999, respectively. As time increases, we observed some increase in the RMSE for both benchmark cases. However, the overall magnitudes are still small (Fig. 1e, f). Overall, we infer that using the exponential spatial discretization from the ELM (which is similar to that in the CLM, Oleson et al., 2013) for the BeTR-v2 solver led to very accurate reproduction of the analytical solutions.

Figure 1 Comparison of the analytical solutions with the numerical solutions derived from BeTR-v2. Panels (a) and (b) are comparisons of simulated and analytical tracer concentrations sampled at the end of 6 different days. In panels (c) and (d), blue dots are plotted using analytical solutions on the x axis and numerical solutions on the y axis; red lines are linear regressions of these data. Panels (e) and (f) are the depth-averaged root-mean-square error (RMSE) through time.

3.2 Single-layer and 1-D column model simulations

In the example standalone simulations using the ecacnp soil BGC, the single-layer and the vertically resolved 10-layer models demonstrate very similar temporal patterns for heterotrophic CO₂ production (Fig. 2a and b). For the 1-D vertically resolved model simulation, surface CO₂ flux is smaller than the heterotrophic CO₂ respiration plus infiltration CO₂ flux throughout the 30 year period. Accordingly, the soil CO₂ concentration builds up continuously, with a seasonal cycle that has its maximum in July and minimum in March (Fig. 2c). We did not run the 1-D model to equilibrium because the example case here is only used for demonstration and to check if the model behavior agrees with our physical intuition. Additionally, this example does not account for root respiration. Therefore, we do not expect the soil CO₂ to be as high as 10 000 ppm as it is often observed (Hirano et al., 2003); however, in global simulations, we do find soil CO₂ below 30 cm could reach as high as a few percent. In a previous study with BeTR-v1 (Tang et al., 2013), we showed that some water-dissolved soil CO₂ can be lost through the hydrological outflow. However, this hydrological loss is small in the above simulation because the forcing data does not include lateral runoff (i.e., lateral runoff is set to zero for standalone simulations for simplicity) and the leaching is negligible. Overall, we conclude that BeTR-v2 consistently represents vertically resolved soil BGC reactions and transport.

Figure 2 (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). The single-layer model used the forcing data for the first layer of the 10-layer model.

3.3 Global simulations

We first ran three global model configurations for a 200 year accelerated spinup using the same parameters (see Koven et al., 2013, for information on how the accelerated spinup is designed). Because there was very little competition among the plant organs for carbon and nutrients, the spinup simulations led to almost identical model outputs between ELMv1-ECA and ELMv1-ECA-V (Fig. 3), suggesting that resolved plant carbon and nutrient allocation using the multiple-flux co-limiting solver has negligible influence and that this solver is properly implemented. ELMv1-BeTR-ECA0 (cyan lines) and ELMv1-ECA (blue lines) produced very different latitudinal patterns of primary flux variables, such as GPP (gross primary productivity; Fig. 3a), ER (ecosystem respiration; Fig. 3c), and therefore NBP (net biome productivity; Fig. 3e), particularly in the tropics. Noticeable differences among fluxes are also found in the middle to high latitudes. Accordingly, the vegetation, soil, and total ecosystem carbon pools differ substantially between ELMv1-ECA and ELM-BeTR-ECA0 (Fig. 3b, d, f) in the tropics and northern middle to high latitudes. We found that the differences are caused by different nutrient uptake rates simulated by the two models. In particular, with the multiple-flux co-limiting solver, both plants and soil carbon processes in ELM-BeTR-ECA0 are overall less nitrogen limited given the same nutrient uptake parameters, but as the simulation proceeds in time, less soil carbon is accumulated and soils become less fertile, consistent with our previous findings (Tang and Riley, 2018).

Figure 3 Comparison of simulated flux and state variables at the end of the 200th year during the accelerated spinup for four model configurations (see Table 1). ELMv1-ECA-V differs from ELMv1-ECA as it uses the multiple-flux co-limiting solver for allocating carbon and nutrients for plant growth. ELMv1-BeTR-ECA0 differs from ELMv1-ECA-V by further solving belowground nutrient competition using the multiple-flux co-limiting solver. ELMv1-BeTR-ECA differs from ELMv1-BeTR-ECA0 by using re-calibrated parameters for two tropical plant functional types: broadleaf evergreen tropical trees and broadleaf deciduous tropical trees. We note that the soil carbon here is smaller than vegetation carbon because factors of accelerated decomposition have not been applied.

Since the simulated differences between ELMv1-BeTR-ECA0 and ELMv1-ECA are caused by different numerical implementations and are greatest in the tropics, we calibrated ELMv1-BeTR-ECA0 with respect to ELMv1-ECA for plant parameters important for plant–soil coupling for two tropical plant functional types (Table 2) and to analyze to what extent these differences can be minimized. We did not use the other plant function types because including them in the calibration is unlikely to change the conclusion of our analysis, as we will see later, while this decision significantly reduced the carbon footprint of this computationally intensive study. We conducted the parameter calibration using predicted total column soil C, vegetation C, GPP, ER, and leaf area index from ELMv1-ECA as constraints. We applied the grid-search approach (e.g., Bergstra and Bengio, 2012) to iteratively minimize model differences until no significant further improvement could be obtained. We ended up with 10 iterations wherein each iteration consists of 50 simulations.

Table 2 Comparison of calibrated (ELMv1-BeTR-ECA) and default (ELMv1-ECA) parameters. PFT-4 is a broadleaf evergreen tropical tree and PFT-6 is a broadleaf deciduous tropical tree.

We found significant differences between default and calibrated model parameter values (Table 2). For ELMv1-BeTR-ECA, the calibrated substrate affinity parameters for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ are of the same magnitudes as their default counterparts. The maximum uptake rates for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and the maximum nitrogen fixation rates are of the same order as their default values. For both PFT-4 (broadleaf evergreen tropical trees) and PFT-6 (broadleaf deciduous tropical trees), a higher amount of fixed nitrogen is directly acquired by plants in ELMv1-BeTR-ECA. When these new parameters were accordingly applied to tropical grid cells, we found smaller differences between ELMv1-ECA and ELMv1-BeTR-ECA for GPP and ER in the tropics (Fig. 3a and c), while differences in other variables (NBP, vegetation carbon, soil carbon, and total ecosystem carbon) remain quite significant (Fig. 3). Therefore, we assert that different numerical implementations result in different model predictions, and such different model predictions are likely not to be rectified by parameter calibration.

To further quantify the modeling uncertainty induced by different numerical implementations, we evaluated the four global simulations using the ILAMB benchmark software (Collier et al., 2018; version 2.5). Consistent with the results shown in Fig. 3, ELMv1-ECA and ELMv1-ECA-V have almost identical performance (Table 3). ELMv1-BeTR-ECA0 performed similarly or worse than ELMv1-ECA, but significantly better for the metric of Global Net Ecosystem Carbon Balance (or Net Biome Production, NBP; see also Fig. 4). Out of eight carbon metrics, ELMv1-BeTR-ECA performed slightly worse than ELMv1-ECA in net ecosystem exchange and significantly worse in vegetation biomass. For the remaining four hydrological cycle metrics, five radiation and energy metrics, and four climate forcing metrics, their performances are almost the same. ELMv1-BeTR-ECA performed slightly better for the metric of Global Net Ecosystem Carbon Balance (also see Fig. 4). For functional relationship benchmarks, ELMv1-BeTR-ECA0, ELMv1-BeTR-ECA, and ELMv1-ECA performed quite similarly (Fig. 5). Notably, ELMv1-BeTR-ECA was able to simulate higher evapotranspiration, resulting in better performance for the functional relationship between GPP and evapotranspiration (red line in Fig. 5a). ELMv1-BeTR-ECA also performed slightly better in some cases, e.g., GPP vs. precipitation (Fig. 5b), evapotranspiration vs. precipitation (Fig. 5d), and leaf area index vs. precipitation (Fig. 5e). 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.

Figure 4 Benchmark of simulated NBP with two global syntheses (a) compared with data from Hoffman et al. (2014) and (b) compared with data from Le Quere et al. (2016). In the linear regression, x is the benchmark and y is the model output.

Figure 5 Benchmark of simulated vs. empirical data-derived variable relationships. The empirical benchmarks are derived from FluxCom output for the period from January 1980 through December 2013 (Jung et al., 2017; Tramontana et al., 2016) and used in ILAMB.

Table 3 Comparison of ILAMB scores for 23 simulated variables by three different model versions. The three entries in bold for ELMv1-BeTR-ECA (and ELMv1-BeTR-ECA0) differ by greater than 5 % from ELMv1-ECA. All metrics are normalized to the range from 0 to 1, where greater values indicate better performance.

Figure 6 Comparison of simulated plant nutrient stresses (on a scale from zero to one) at year 2010. Panels (a), (c), and (e) are for nitrogen stress, and panels (b), (d), and (e) are for phosphorus stress. Greater values indicate stronger limitation for the corresponding nutrient.

Nevertheless, there are significant differences in the simulated nutrient dynamics (as reflected in the fast variables that have a significant seasonal cycle), which over a longer time period can affect the cumulative carbon state variables (e.g., right column in Fig. 3). For instance, for the simulation year 2010, the hydrological nitrogen loss estimated by ELMv1-ECA is 56.8 TgN yr⁻¹, while that estimated by ELMv1-BeTR-ECA is ∼10 times smaller (4.2; the value estimated by ELMv1-BeTR-ECA0 is 4.5 TgN yr⁻¹). The global N₂O emission estimated by ELMv1-ECA is 13.2, whereas ELMv1-BeTR-ECA estimated 14.3 (and the value by ELMv1-BeTR-ECA0 is 12.5 TgN yr⁻¹).

Further, ELMv1-ECA and ELMv1-BeTR-ECA predicted different spatial distributions of plant nutrient limitations. For example, ELMv1-ECA predicted higher nitrogen limitation to plant productivity in northern high-latitude regions and some regions in northwestern Australia (Fig. 6e) for the simulation year 2010, even though the general patterns of plant nitrogen limitation are very similar (Fig. 6a, b). However, compared to ELMv1-BeTR-ECA, ELMv1-ECA predicted much more widespread plant phosphorus limitation (Fig. 6f). 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). This type of discrepancy indicates that simulations under a future climate with these two seemingly similar models may very well diverge.

In summary, adopting different numerical coupling strategies resulted in different model behavior, even with many similar C cycle metrics. In particular, we found important differences in N and P dynamics that could impact 21st century simulations when atmospheric CO₂ increases are expected to cause progressive nutrient limitations.

4 Conclusion

In this study, we present BeTR-v2, a reactive transport module for flexible implementation of soil biogeochemistry in ecosystem models. We benchmarked the numerical solver with two analytical solutions and demonstrated that BeTR-v2 can support soil biogeochemical modeling in single-layer, 1-D vertically resolved single column and global simulation modes. We implemented the soil biogeochemistry formulation from ELMv1-ECA into ELMv1-BeTR-ECA and explored the role of numerical solver robustness and calibration on carbon and nutrient cycling. 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), ELMv1-BeTR-ECA has to use different model parameters to produce similar model behavior (as compared to ELMv1-ECA) for carbon and water cycling. However, the calibration was not able to correct the important differences in predicted N and P cycling, indicating that 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.