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 $tto0$ 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.