The difference between a mechanistic model of a spatio‐temporal process and associated observations can reasonably be represented by a random process with possible dependence structure in space and time. Often, such error processes are assumed to be stationary in time, so that a unique spatial error covariance structure is applicable for all time. However, model performance can be related to large‐scale time‐varying features, causing the error process to be nonstationary, leading to a time‐varying spatial covariance. Because of dimensionality constraints, one cannot estimate such time‐varying covariances reliably, and therefore they must be modeled. Such models for covariances should be able to accommodate influences from exogenous covariates. However, the inclusion of covariates in a standard covariance model is challenging due to the fact that the covariance must be positive definite. We mitigate this issue by modeling the time‐varying spatial dependence through the Fourier coefficients of the log spectrum, as originally developed for time series representations. The convenient feature of these so‐called “cepstral” models (or “exponential models”) is that the associated parameters are unrestricted yet still are guaranteed to produce a positive definite covariance. Thus, it is fairly straightforward to model these parameters as time‐varying random processes that depend on other covariates within the Bayesian hierarchical modeling paradigm. We illustrate the effectiveness of this model through various simulated examples and by applying it to differences from a long‐lead sea surface temperature forecast model. In this case, we consider the Pacific Decadal Oscillation and Southern Oscillation Index as possible factors that influence the error covariance. Copyright © 2014 John Wiley & Sons, Ltd.