We present for the first time Next-to-Leading (NLO) QCD renormalization group (RG) evolution matrices for nonleptonic $\Delta F=2$ transitions in the Standard Model effective field theory (SMEFT). To this end we transform first the known two-loop QCD anomalous dimension matrices (ADMs) of the BSM (Beyond the SM) operators in the so-called Buras Misiak Urban basis into the ones in the common weak effective theory (WET) basis (the so-called Jenkins Manohar Stoffer basis) for which tree-level and one-loop matching to the SMEFT are already known. This subsequently allows us to find the two-loop QCD ADMs for the SMEFT nonleptonic $\Delta F=2$ operators in the Warsaw basis. Having all these ingredients we investigate the impact of these NLO QCD effects on the QCD RG evolution of SMEFT Wilson coefficients for nonleptonic $\Delta F=2$ transitions from the new physics scale $\Lambda$ down to the electroweak scale ${\mu }_{\mathrm{ew}}$. The main benefit of these new contributions is that they allow one to remove renormalization scheme dependences present in the one-loop matchings both between the WET and SMEFT and also between SMEFT and a chosen UV completion. But the Next-to-Leading (NLO) QCD effects, calculated here in the Naive dimensional regularisation minimal subtraction scheme, turn out to be small, in the ballpark of a few percent but larger than one-loop Yukawa top effects when only the $\Delta F=2$ operators are considered. The more complicated class of nonleptonic $\Delta F=1$ decays will be presented soon in another publication.