The gradient flow exponentially suppresses ultraviolet field fluctuations and removes ultraviolet divergences (up to a multiplicative fermionic wavefunction renormalization). It can be used to describe real-space Wilsonian renormalization group transformations and determine the corresponding beta function. We propose a new nonperturbative renormalization scheme for local composite fermionic operators that uses the gradient flow and is amenable to lattice QCD calculations. We present preliminary nonperturbative results for the running of quark bilinear operators in this scheme and outline the calculation of perturbative matching to the MS-bar scheme.