Physics-based computational models increasingly inform high-stakes predictions across a wide range of scientific and engineering applications. However, multiple competing models at different abstraction levels and theoretical frameworks necessitate model selection for reliable prediction. Bayesian model selection (BMS) tackles this challenge by quantifying model evidence given data via the marginal likelihood. This quantity inherently penalizes complexity, thereby favoring the simplest model that explains the data [1]. Estimating the marginal likelihood involves high-dimensional integration over parameter space, which becomes computationally intensive for physics-based models with complex posterior geometries. Consequently, physics-based engineering applications often either rely on suboptimal approximations or employ specialized alternatives, which impose rigid constraints on the workflow and thereby limit their integration into existing Bayesian workflows. We present a unified Bayesian model selection workflow for physics-based models that integrates advanced marginal likelihood estimation with Gaussian process (GP) emulation in a modular, end-to-end software pipeline. Through systematic evaluation of state-of-the-art marginal likelihood estimators, we identified the Learned Harmonic Mean Estimator with Normalizing Flows (LHME) [2] as a robust and scalable choice. LHME estimates the marginal likelihood using posterior samples and is agnostic to the underlying MCMC algorithm, enabling seamless integration into existing Bayesian workflows. However, posterior sampling becomes prohibitively expensive for physics-based models. We overcome this limitation by leveraging GP emulation to replace expensive models with fast-to-evaluate surrogates, enabling efficient and flexible BMS. We demonstrate the model selection workflow through three specific engineering applications: cryobot thermal modeling, blood hemolysis modeling, and avalanche runout simulation. Together, these applications illustrate the applicability of our approach across domains, real-world datasets, and computationally expensive simulators. Integrating powerful and versatile marginal likelihood estimation with surrogate modelling in a unified framework makes rigorous Bayesian model selection practical for complex engineering problem settings that were previously inaccessible.