Physics-based simulations of earthquake ground motion are useful to complement recorded ground motions. However, the computational expense of performing numerical simulations hinders their applicability to tasks that require real-time solutions or ensembles of solutions for different earthquake sources. To enable rapid physics-based solutions, we present a reduced-order modeling approach based on interpolated proper orthogonal decomposition (POD) to predict peak ground velocities (PGVs). As a demonstrator, we consider PGVs from regional 3D wave propagation simulations at the location of the 2008 MW 5.4 Chino Hills earthquake using double-couple sources with varying depth and focal mechanisms. These simulations resolve frequencies ≤1.0 Hz and include topography, viscoelastic attenuation, and S-wave speeds ≥500 m/s. We evaluate the accuracy of the interpolated POD reduced-order model (ROM) as a function of the approximation method. Comparing the radial basis function (RBF), multilayer perceptron neural network, random forest, and k-nearest neighbor, we find that the RBF interpolation gives the lowest error (≈0.1 cm/s) when tested against an independent data set. We also find that evaluating the ROM is 107–108 times faster than the wave propagation simulations. We use the ROM to generate PGV maps for 1 million different focal mechanisms, in which we identify potentially damaging ground motions and quantify correlations between focal mechanism, depth, and accuracy of the predicted PGV. Our results demonstrate that the ROM can rapidly and accurately approximate the PGV from wave propagation simulations with variable source properties, topography, and complex subsurface structure.