We demonstrate the use of the gradient flow combined with the short flow-time expansion (GF+SFTX) as a renormalization procedure for four-quark operator matrix elements and associated bag parameters relevant to neutral heavy-meson mixing ($\Delta Q=2$) and heavy-meson lifetimes ($\Delta Q=0$). Using six RBC/UKQCD 2+1-flavor domain-wall fermion ensembles, we calculate for a charm-strange system with physical quark masses flowed bag parameters and match them to the $\overline{\mathrm{\text{MS}}}$ scheme using perturbative SFTX coefficients up to next-to-next-to-leading order in QCD. We employ a multi-scale matching strategy and a renormalization-group improved flow-time evolution which allows for a reliable estimate of systematic uncertainties. For a fictitious neutral $D_s$ meson, we obtain the $\Delta Q=2$ $\overline{\mathrm{\text{MS}}}$ bag parameter ${\cal B}^{\overline{\text{MS}}}_1(3\,{\rm GeV})=0.7673(123)$, consistent with existing short-distance $D^0$ mixing determinations. For the $\Delta Q=0$ lifetime-ratio operator basis, we find the $\overline{\mathrm{\text{MS}}}$ results $B^{\overline{\text{MS}}}_1(3\,{\rm GeV})=1.0524(97)$, $B^{\overline{\text{MS}}}_2(3\,{\rm GeV})=0.9621(71)$, $\epsilon^{\overline{\text{MS}}}_1(3\,{\rm GeV})=-0.2275(76)$, and $\epsilon^{\overline{\text{MS}}}_2(3\,{\rm GeV})=-0.0005(8)$. We provide conversion formulae to re-express these results for an arbitrary choice of evanescent operators. These results demonstrate that GF+SFTX can deliver precise determinations of dimension-six four-quark operators and establish a framework for future lattice computations including more complex operator bases, where the challenge of power-divergent mixing is shifted to the continuum and handled in the SFTX.