Accurate prediction of turbulent flow is crucial for urban built environments. Traditional reduced-order models (ROMs), such as proper orthogonal decomposition (POD), often struggle to capture the nonlinear temporal dynamics, while recent machine-learning-based spatiotemporal predictors suffer from high training costs and limited physical interpretability. To address these challenges, we propose a POD-based Fourier Neural Operator (POD–FNO), which evolves dynamics in an interpretable POD modal space while retaining the nonlinear expressiveness of neural operators. This design enables accurate and stable temporal prediction with substantially improved computational efficiency, making it particularly suitable for data-intensive urban flow scenarios. The model is evaluated in the multi-bluff-body wake region and compared with representative models, including FNO, Long Short-Term Memory (LSTM), and transformer-based architectures, as well as autoencoder- and POD-based ROM formulations. The results demonstrate that POD–FNO effectively combines the computational efficiency of POD-based modeling with the superior accuracy and long-horizon stability of FNO, consistently outperforming other temporal predictors. In addition, we analyze computational complexity, showing how spectral convolution and mode truncation drive efficiency. Beyond predictive performance and efficiency, we further investigate the physics-aware interpretability of POD–FNO by analyzing the mode importance and coupling effects, offering insights into the learned nonlinear temporal dynamics under complex flow. The framework is further applied to multiple-bluff-body layout configurations, demonstrating robust performance and scalable efficiency.