Probabilistic modeling of multidimensional spatiotemporal data is critical to many real-world applications. However, real-world spatiotemporal data often exhibits complex dependencies that are nonstationary, i.e., correlation structure varies with location/time, and nonseparable, i.e., dependencies exist between space and time. 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 paper, we propose a new statistical framework -- Bayesian Complementary Kernelized Learning (BCKL) -- to achieve scalable probabilistic modeling for multidimensional spatiotemporal data. To effectively describe complex dependencies, BCKL integrates kernelized low-rank factorization with short-range spatiotemporal Gaussian processes (GP), in which the two components complement each other. Specifically, 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. We develop an efficient Markov chain Monte Carlo (MCMC) algorithm for model inference and evaluate the proposed BCKL framework on both synthetic and real-world spatiotemporal datasets. Our results confirm the superior performance of BCKL in providing accurate posterior mean and high-quality uncertainty estimates.