Advances in physics-based earthquake simulations, utilizing high-performance computing, have been exploited to better understand the generation and characteristics of the high-frequency seismic wavefield. However, direct comparison to ground motion observations of a specific earthquake is challenging. We here propose a new approach to simulate data-fused broadband ground motion synthetics using 3D dynamic rupture modeling of the 2016 Mw 6.2 Amatrice, Italy earthquake. We augment a smooth, best-fitting model from Bayesian dynamic rupture source inversion of strong-motion data (<1 Hz) with fractal fault roughness, frictional heterogeneities, viscoelastic attenuation, and topography. The required consistency to match long periods allows us to quantify the role of small-scale dynamic source heterogeneities, such as the 3D roughness drag, from observational broadband seismic waveforms. We demonstrate that 3D data-constrained fully dynamic rupture synthetics show good agreement with various observed ground-motion metrics up to ∼5 Hz and are an important avenue toward non-ergodic, physics-based seismic hazard assessment.