The inconsistency between the fixed-order (FO) and contour-improved (CI) representation of the QCD corrections to the inclusive hadronic tau decay width limits the precision to which the strong coupling can be determined from this process. It has recently been shown that subtracting the infrared renormalon divergence related to the gluon
condensate resolves the discrepancy. Here we suggest to employ the gradient flow to define gauge-invariant regularized operators and to use the corresponding condensates in the operator product expansion. The associated rearrangement of the perturbative series results in automatic renormalon subtraction without the need to determine explicitly the Stokes constants that normalize the divergent asymptotic series. Applying this method to the gluon condensate, we
find that the CI series is modified and now agrees with the (unmodified) FO series. This conclusively demonstrates the preference for the fixed-order approach, as has been advocated long ago.