Developing effective and computationally efficient statistical models to accommodate nonstationary/nonseparable processes containing both long-range and short-scale variations becomes a challenging task, especially for large-scale datasets with various corruption/missing structures. In this talk, we will introduce Bayesian Complementary Kernelized Learning (BCKL) to achieve scalable probabilistic modeling for multidimensional spatiotemporal data. BCKL integrates kernelized low-rank factorization with short-range spatiotemporal Gaussian processes (GP), in which the two components complement each other. We use a multi-linear low-rank factorization component to capture the global/long-range correlations in the data and introduce an additive short-scale GP based on compactly supported kernel functions to characterize the remaining local variabilities. Our results confirm the superior performance of BCKL in providing accurate posterior mean and high-quality uncertainty estimates.