<p>Simulation models of nonlinear dynamic systems may not accurately predict the complex dynamic behavior of the system due to model simplification, assumption, and/or uncertain model parameters. While approaches have been developed to address model-form error in the context of physics-guided machine learning, the simultaneous model discrepancy quantification and model parameter calibration remains a long-standing challenge for nonlinear dynamic systems. Inspired by the seminal work of Kennedy and O’Hagan (KOH), this paper proposes a modularized, versatile model uncertainty quantification framework specifically designed for nonlinear dynamic systems. It does not require the explicit form of the differential equations that describe the system dynamics and can be applied to systems modeled by complex computer simulations or commercial software. The framework consists of four main modules. In the first module, a Spectral-normalized Neural Gaussian Process (SNGP) model is developed to represent the simulation model of the nonlinear dynamic system as a nonlinear autoregressive exogenous (NARX) model while accounting for surrogate prediction uncertainty. In the second module, a hybrid surrogate emulator is created by training a second SNGP model to learn the model bias of the NARX model in the presence of model parameter uncertainty based on experimental data. These two SNGP models are integrated to form a hybrid surrogate emulator capable of predicting system responses with quantified uncertainty for any given realization of uncertain model parameters. In the third module, the resulting hybrid surrogate emulator is employed as a forward model to calibrate the uncertain model parameters using a novel likelihood-free and computationally efficient Bayesian inference method based on normalizing flow. Finally, the fourth module combines the updated model parameters and the hybrid surrogate emulator to predict system responses for untested input conditions while accounting for uncertainty in both model parameters and the surrogate models. The proposed methodology is demonstrated and validated using two engineering benchmark problems, namely a single-degree-of-freedom nonlinear oscillator and a six-story nonlinear shear-building model. The results show that the proposed method can effectively enhance the prediction accuracy of the simulation model by quantifying both the model-form error and model parameter uncertainty.</p>