The use of physics-based earthquake simulators is increasingly common in earthquake forecasting for seismic hazard. Their popularity is related to their ability to overcome completeness limitations of real seismicity catalogues, while reproducing complex earthquake rupture behavior and interaction patterns through the modelling of the physical processes involved in earthquake nucleation and rupture propagation. One common challenge when designing earthquake simulations is the selection of the input parameters that will produce the most feasible models in terms resemblance to natural earthquake processes and relationships, e.g., the rate-and-state frictional parameters – a, b – and the initial normal stress. The frequent lack of empirical data on such parameters, often bases their selection on non-systematic testing and qualitative model performance analysis, thus potentially reducing the objectivity of the modelling. We present a new quantitative approach to evaluate and rank the performance of multiple earthquake simulation models based on a workflow that scores each synthetic catalog according to their combined fit to objective seismological benchmarks. These benchmarks rely on widely used empirical earthquake data: 1) scaling relationships, 2) shape of the magnitude-frequency distribution and 3) rates of surface ruptures from paleoseismology. The approach permits an objective and effective approximation to model performance evaluation, allowing to easily identify which models (and input parameter combinations) simulate better natural earthquake relations and behavior. The algorithm-based approach also facilitates the exhaustive analysis of many input parameter combinations and allows the identification of systematic correlations between parameters and model performance. We validate the approach with earthquake simulations on a theoretical planar fault and with published simulations at the Eastern Betics Shear Zone (EBSZ) in southeastern Spain. In both cases, the method successfully ranks the models with better resemblance to natural catalogues, while avoiding self-correlation of the benchmark scores, i.e., the best model is not the best in all benchmarks but the better balanced across them. In the case of the EBSZ, our ranking analysis replicates the qualitative and manual analyses previously published, which reinforces the usefulness of the approach. We also identify very clear correlations between the model performance and the rate-and-state a and b parameters. In particular, we observe that larger differences between a and b tend to better model performance. Conversely, the initial normal stress does not correlate with the performance. Overall, we estimate that the approach can ease researchers on earthquake simulation design and building, and on better understanding the impact of selected input parameters into the physics-based models. Moreover, the model ranking results can be employed to drive further analysis such as weighting of earthquake forecast models in seismic hazard logic trees.