We study the possibility of using multilevel algorithms for the computation of correlation functions of gradient flow observables. For each point in the correlation function, an approximate flow is defined which depends only on links in a subset of the lattice. Together with a local action, this allows for independent updates and consequently a convergence of the Monte Carlo process faster than the inverse square root of the number of measurements. We demonstrate the feasibility of this idea in the correlation functions of the topological charge and the energy density.