We perform a scale setting procedure of a mixed action setup consisting of valence Wilson twisted mass fermions at maximal twist on CLS ensembles with $N_f=2+1$ flavours of $O(a)$-improved Wilson sea quarks. We determine the gradient flow scale $t_0$ using pion and kaon masses and decay constants in the isospin symmetric limit of QCD as external physical input. We provide an update for our results with ensembles at the physical pion mass and lattice spacings down to $a\approx0.039$ fm. We employ model variation techniques to probe the systematic uncertainties in the extraction of the ground state signal of lattice observables, as well as for the continuum-chiral extrapolations used to determine $t_0$ at the physical point. We determine $t_0$ from a combination of the unitary Wilson lattice regularization and the mixed action regularization, which leads to an improved control of the systematic and statistical uncertainties.