The gradient-flow operator product expansion for QCD current correlators including operators up to mass dimension four is calculated through NNLO. This paves an alternative way for efficient lattice evaluations of hadronic vacuum polarization functions. In addition, flow-time evolution equations for flowed composite operators are derived. Their explicit form for the non-trivial dimension-four operators of QCD is given through order $alpha_s^3$.