Using reactive transport codes to provide mechanistic biogeochemistry representations in land surface models: A proof of concept with CLM-PFLOTRAN

We explore coupling to a configurable subsurface reactive transport code as a flexible and extensible approach to biogeochemistry in land surface models. A reaction network with the CLM-CN decomposition, nitrification, denitrification, and plant uptake similar to CLM4.5 is used as an example. We implement the reactions in the open-source PFLOTRAN code, coupled with the Community Land Model (CLM), and test at Arctic, temperate, and tropical sites. To make the reaction network designed for use in explicit time stepping in CLM compatible with the implicit time stepping used in PFLOTRAN, the Monod substrate rate-limiting function with a residual concentration is used to represent the limitation of nitrogen availability on plant uptake and immobilization. Switching from explicit to implicit methods increases numerical rigor but introduces computational challenges. Our objective is to achieve accurate, efficient, and robust numerical solutions to demonstrate the feasibility of CLM-PFLOTRAN soil biogeochemistry. Our results suggest that care needs to be taken to use scaling, clipping, or log transformation to avoid negative concentrations during the Newton iterations. With a tight relative update tolerance to avoid false convergence, an accurate solution can be achieved with about 50 {%} more computing time than CLM in point mode site simulations using either the scaling or clipping methods. The log transformation method takes 60-100 {%} more computing time than CLM. The computing time increases slightly for clipping and scaling; it increases substantially for log transformation for half saturation decrease from 10^-3 to 10^-9 \unit{mol m^-3}, which normally results in decreasing nitrogen concentrations. The frequent occurrence of very low concentrations (e.g. below nanomolar) can increase the computing time for clipping or scaling by about 20 {%}; computing time can be doubled for log transformation. Caution needs to be taken in choosing the appropriate scaling factor because a small value caused by a negative update to a small concentration may diminish the update and result in false convergence even with very tight relative update tolerance. As some biogeochemical processes (e.g., methane and nitrous oxide production and consumption) involve very low half saturation and threshold concentrations, this work provides insights for addressing nonphysical negativity issues and facilitates the representation of a mechanistic biogeochemical description in earth system models to reduce climate prediction uncertainty.