The earthquake hazard and seismic risk in Iceland are the highest in the Southwest due to transform faulting in the South Iceland Seismic Zone (SISZ) and Reykjanes Peninsula Oblique Rift (RPOR). The reliable probabilistic seismic hazard assessment (PSHA) is therefore critical in this region which in turn requires two of the key elements for a reliable PSHA. Namely, the most appropriate and physically consistent ground motion models (GMMs) and the specification of seismic sources. In this study, we address this by employing a suite of new hybrid Bayesian empirical GMMs and a new physics-based finite-fault system model for the SISZ-RPOR system. By ranking the GMMs using the deviance information criterion against the Icelandic strong-motion dataset we propose backbone GMMs along with their Δ-factors that comprehensively capture the characteristics of Icelandic data. We then simulate suites of synthetic finite-fault earthquake catalogues that are consistent with the physics-based fault system. We then carry out a Monte-Carlo PSHA using the backbone approach for two representative near-fault and far-field sites in Southwest Iceland, but this approach provides a more realistic and robust treatment of the extent and effects of GMM's epistemic uncertainty on the PSHA. We show that the backbone PSHA is consistent with point estimates using the classical PSHA approach, but importantly they reveal the "body" i.e., the range, of realistically expected hazard values that dwarf any differences between the results from the two approaches. Thus, the backbone models of this study along with their Δ-factors are the ideal candidates for use in future PSHA for Iceland.