Improving models for turbulent fluctuations and engineering quantities of interest (e.g. drag, equivalent sand-grain roughness) over general rough surfaces is important because fully resolved simulations of high-Reynolds-number flows are computationally impractical and classical drag correlations may exhibit high error when applied to surfaces outside of the training cases. Data-driven and physics-based modeling methods can provide improved results over classical drag correlations, with the caveat that out-of-sample predictions and model interpretability remain challenging for machine learning-based models. These methods do not predict the spatially-varying modulation of the turbulence by the roughness, which contributes strongly to the momentum, energy, and scalar transport effects observed in rough wall flows. Models of roughness sublayer turbulence have been developed for simplified single-scale roughness geometries using laminar models and resolvent analysis. Engineering-relevant roughness is often multiscale, and providing predictions for these surfaces is an ongoing challenge. In this work, a physics-based reduced-order model for mean and fluctuating quantities in roughness sublayer turbulence over engineering-relevant multiscale roughness is developed using resolvent analysis with a volume penalization term (RAVP) which is uniquely determined for any given roughness geometry. The efficacy of the model for predicting wake field and mean quantities is evaluated over a canonical sand-grain roughness geometry. By decomposing the volume penalizing term into mean and spatially varying components, RAVP can be computed on a wavenumber-by-wavenumber basis, enabling computationally inexpensive and quantitatively useful predictions for turbulent fluctuations and their statistics. An iterative framework for mean flow modeling and drag prediction is presented which incorporates the volume-penalizing term and dispersive stresses from RAVP as well as an eddy viscosity adjusted for surface roughness.