This article proposes a physical-statistical modeling approach for spatio-temporal data arising from a class of stochastic convection-diffusion processes. Such processes are widely found in scientific and engineering applications where fundamental physics imposes critical constraints on how data can be modeled and how models should be interpreted. The idea of spectrum decomposition is employed to approximate a physical spatio-temporal process by the linear combination of spatial basis functions and a multivariate random process of spectral coefficients. Unlike existing approaches assuming spatially and temporally invariant convection-diffusion, this article considers a more general scenario with spatially varying convection-diffusion and nonzero-mean source-sink. As a result, the temporal dynamics of spectral coefficients is coupled with each other, which can be interpreted as the nonlinear energy redistribution across multiple scales from the perspective of physics. Because of the spatially varying convection-diffusion, the space-time covariance is nonstationary in space. The theoretical results are integrated into a hierarchical dynamical spatio-temporal model. The connection is established between the proposed model and the existing models based on integro-difference equations. Computational efficiency and scalability are also investigated to make the proposed approach practical. The advantages of the proposed methodology are demonstrated by numerical examples, a case study, and comprehensive comparison studies. Computer code is available on GitHub. Supplementary materials for this article are available online.