Soil biogeochemical models (SBMs) simulate element transfer processes between organic soil pools. These models can be used to specify falsifiable quantitative assertions about soil system dynamics and their responses to global surface temperature warming. To determine whether SBMs are useful for representing and forecasting data-generating processes in soils, it is important to conduct data assimilation and fitting of SBMs conditioned on soil pool and flux measurements to validate model predictive accuracy. SBM data assimilation has previously been carried out in approaches ranging from visual qualitative tuning of model output against data to more statistically rigorous Bayesian inferences that estimate posterior parameter distributions with Markov chain Monte Carlo (MCMC) methods. MCMC inference is better able to account for data and parameter uncertainty, but the computational inefficiency of MCMC methods limits their ability to scale assimilations to larger data sets. With formulation of efficient and statistically rigorous SBM inference frameworks remaining an open problem, we demonstrate the novel application of a variational inference framework that uses a method called normalizing flows to approximate SBMs that have been discretized into state space models. We fit the approximated SBMs to synthetic data sourced from known data-generating processes to identify discrepancies between the inference results and true parameter values and ensure functionality of our method. Our approach trades estimation accuracy for algorithmic efficiency gains that make SBM data assimilation more tractable and achievable under computational time and resource limitations.