Abstract: Soil biogeochemical models (SBMs) are an important tool used by Earth scientists to quantify the impact of rising global surface temperatures. SBMs represent the soil carbon and microbial dynamics across time as differential equations, and inference on model parameters is conducted to project changes in parameter values under warming climate conditions. Traditionally, the field has relied on MCMC algorithms for posterior inference, often implemented via probabilistic programming languages like Stan. However, computational cost makes it difficult to scale MCMC methods to more complex SBM models and large-scale datasets. In this paper, we develop variational inference methods for time-discretized SBMs as an alternative to MCMC. We propose an efficient family of variational approximations based on Gauss-Markov distributions that leverages the temporal structure of sequential models, scaling linearly in both time and space with respect to the sequence length. We show in experiments with simulated data and real CO$_2$ response ratios that our approach converges faster, and recovers posterior that more accurately captures uncertainty than previous variational methods. Our black-box inference approach is designed to integrate with probabilistic programming languages to enable future scientific applications.