We calculate the dimension-six $\Delta Q=0$ four-quark matrix elements describing heavy-meson lifetime ratios using the gradient flow with its short flow-time expansion as a renormalization procedure. On six RBC/UKQCD 2+1-flavor domain-wall fermion ensembles, we determine flowed bag parameters for physical charm and strange quarks and match to the $\overline{\mathrm{\text{MS}}}$ scheme with perturbative short flow-time expansion coefficients through next-to-next-to-leading order (NNLO). A multi-scale matching procedure using renormalization-group running improves the extrapolation to zero flow time. For the operators relevant to $\tau (D_s)/\tau (D^0)$ at the SU(3)$_{\rm F}$ symmetric point, we obtain $B_1^{\overline{\text{MS}}}(3\,{\rm GeV})=1.0524(97)$,$B_2^{\overline{\text{MS}}}(3\,{\rm GeV})=0.9621(70)$, $\epsilon_1^{\overline{\text{MS}}}(3\,{\rm GeV})=-0.2275(76)$, and $\epsilon_2^{\overline{\text{MS}}}(3\,{\rm GeV})=-0.0005(8)$ using a specific choice of evanescent operators. This is the first lattice-QCD determination of $\Delta Q=0$ four-quark operators with a full error budget. It opens the path towards higher-precision predictions of heavy-meson lifetimes and similar quantities exhibiting operator mixing under renormalization.