The land surface processes simulated within the E3SM model are central to understanding the spatio-temporal variations of climate change over the next century. However, given the modeling complexity, the E3SM Land Model (ELM) utilizes several simplifying physical assumptions and parameterizations that lead to both significant spread in the predictions and deviations from observed data. This study introduces a non-intrusive, statistical approach to both calibrating and quantifying the modeling uncertainties in the ELM. Specifically, the Karhunen-Loève Decomposition (KLD) procedure is used to identify a linear low-dimensional latent space that captures the high-dimensional dynamics of the spatio-temporal data, followed by a computationally inexpensive surrogate approximation in the latent space. A statistical model error term is then introduced and embedded directly into the surrogate model. A Bayesian framework is then used with low-cost samples from the surrogate and observational data from select FLUXNET sites across the US to calibrate both the model error parametrization and ELM input parameters values that best predict observed data. This approach is tested on several quantities of interest, temporal scales, and in sparse datasets to understand its robustness and performance. The results from this study not only yield insight into the relative magnitudes for various forms of uncertainty that impact ELM predictions, but also showcase a novel approach to performing structural error estimation in computationally expensive Earth system models.