  The use of spectral proper orthogonal decomposition (SPOD) to construct
low-order models for broadband turbulent flows is explored. The choice of SPOD
modes as basis vectors is motivated by their optimality and space-time
coherence properties for statistically stationary flows. This work follows the
modeling paradigm that complex nonlinear fluid dynamics can be approximated as
stochastically forced linear systems. The proposed stochastic two-level
SPOD-Galerkin model governs a compound state consisting of the modal expansion
coefficients and forcing coefficients. In the first level, the modal expansion
coefficients are advanced by the forced linearized Navier-Stokes operator under
the linear time-invariant assumption. The second level governs the forcing
coefficients, which compensate for the offset between the linear approximation
and the true state. At this level, least squares regression is used to achieve
closure by modeling nonlinear interactions between modes. The statistics of the
remaining residue are used to construct a dewhitening filter that facilitates
the use of white noise to drive the model. If the data residue is used as the
sole input, the model accurately recovers the original flow trajectory for all
times. If the residue is modeled as stochastic input, then the model generates
surrogate data that accurately reproduces the second-order statistics and
dynamics of the original data. The stochastic model uncertainty,
predictability, and stability are quantified analytically and through Monte
Carlo simulations. The model is demonstrated on large eddy simulation data of a
turbulent jet at Mach number $M=0.9$ and Reynolds number of $Re_D\approx 10^6$.
