Flaring activity following gamma-ray bursts (GRBs), observed in both long and short GRBs, signals a long-term activity of the central engine. However, its production mechanism has remained elusive. Here we develop a quantitative model of the idea proposed by Perna et al. of a disc whose outer regions fragment due to the onset of gravitational instability. The self-gravitating clumps migrate through the disc and begin to evolve viscously when tidal and shearing torques break them apart. Our model consists of two ingredients: theoretical bolometric flare lightcurves whose shape (width, skewness) is largely insensitive to the model parameters, and a spectral correction to match the bandpass of the available observations, that is calibrated using the observed spectra of the flares. This simple model reproduces, with excellent agreement, the empirical statistical properties of the flares as measured by their width-toarrival time ratio and skewness (ratio between decay and rise time). We present model fits to the observed lightcurves of two well-monitored flares, of GRB 060418 and of GRB 060904B. To the best of our knowledge, this is the first quantitative model able to reproduce the flare lightcurves and explain their global statistical properties.