Gradient flow is a smoothing procedure that suppresses ultraviolet fluctuations of gauge fields. It is often used for high-precision scale setting and renormalization of operators in lattice QCD calculations. The gradient flow equation is defined on the SU(3) manifold and therefore requires geometric, or structure-preserving, integration methods to obtain its numerical solutions. I discuss the properties of the three-stage third-order Runge-Kutta integrator introduced by Luescher (that became almost the default choice in lattice QCD applications) and its relation to structure-preserving integrators available in the literature. I demonstrate how classical low-storage Runge-Kutta methods can be turned into structure-preserving integration methods and how schemes of order higher than three can be built. Based on the properties of the low-storage schemes I discuss how the methods can be tuned for optimal performance in lattice QCD or any other applications.