  We study thermodynamic properties of Nf=2+1 QCD on the lattice adopting
O(a)-improved Wilson quark action and Iwasaki gauge action. To cope with the
problems due to explicit violation of the Poincare and chiral symmetries, we
apply the Small Flow-time eXpansion (SFtX) method based on the gradient flow,
which is a general method to correctly calculate any renormalized observables
on the lattice. In this method, the matching coefficients in front of operators
in the small flow-time expansion are calculated by perturbation theory. In a
previous study using one-loop matching coefficients, we found that the SFtX
method works well for the equation of state, chiral condensates and
susceptibilities. In this paper, we study the effect of two-loop matching
coefficients by Harlander et al. We also test the influence of the
renormalization scale in the SFtX method. We find that, by adopting the mu_0
renormalization scale of Harlander et al. instead of the conventional
mu_d=1/sqrt{8t} scale, the linear behavior at large t is improved so that we
can perform the t -> 0 extrapolation of the SFtX method more confidently. In
the calculation of the two-loop matching coefficients by Harlander et al., the
equation of motion for quark fields was used. For the entropy density in which
the equation of motion has no effects, we find that the results using the
two-loop coefficients agree well with those using one-loop coefficients. On the
other hand, for the trace anomaly which is affected by the equation of motion,
we find discrepancies between the one- and two-loop results at high
temperatures. By comparing the results of one-loop coefficients with and
without using the equation of motion, the main origin of the discrepancies is
suggested to be attributed to O((aT)^2)=O(1/N_t^2) discretization errors in the
equation of motion at N_t =< 10.
