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 Mw6.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 at long periods allows us to quantify the role of 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 towards non-ergodic, physics-based seismic hazard assessment.