  Recently, Harlander et al.\ [Eur.\ Phys.\ J.\ C {\bf 78}, 944 (2018)] have
computed the two-loop order (i.e., NNLO) coefficients in the gradient-flow
representation of the energy--momentum tensor (EMT) in vector-like gauge
theories. In this paper, we study the effect of the two-loop order corrections
(and the three-loop order correction for the trace part of the EMT, which is
available through the trace anomaly) on the lattice computation of
thermodynamic quantities in quenched QCD. The use of the two-loop order
coefficients generally reduces the $t$~dependence of the expectation values of
the EMT in the gradient-flow representation, where $t$~is the flow time. With
the use of the two-loop order coefficients, therefore, the $t\to0$
extrapolation becomes less sensitive to the fit function, the fit range, and
the choice of the renormalization scale; the systematic error associated with
these factors is considerably reduced.
