We report on a sub-percent scale determination using the omega baryon mass and gradient-flow methods. The calculations are performed on 22 ensembles of $N_f=2+1+1$ highly improved, rooted staggered sea-quark configurations generated by the MILC and CalLat Collaborations. The valence quark action used is Möbius Domain-Wall fermions solved on these configurations after a gradient-flow smearing is applied with a flowtime of $t_\rm gf=1$ in lattice units. The ensembles span four lattice spacings in the range $0.06 \lesssim a \lesssim 0.15$ fm, six pion masses in the range $130 \lesssim m_\pi \lesssim 400$ MeV and multiple lattice volumes. On each ensemble, the gradient-flow scales $t_0/a^2$ and $w_0/a$ and the omega baryon mass $a m_\Omega$ are computed. The dimensionless product of these quantities is then extrapolated to the continuum and infinite volume limits and interpolated to the physical light, strange and charm quark mass point in the isospin limit, resulting in the determination of $\sqrtt_0=0.1422(14)$ fm and $w_0 = 0.1709(11)$ fm with all sources of statistical and systematic uncertainty accounted for. The dominant uncertainty in this result is the stochastic uncertainty, providing a clear path for a few-per-mille uncertainty, as recently obtained by the Budapest-Marseille-Wuppertal Collaboration.