The flow over an open cavity is an example of supercritical Hopf bifurcation leading to periodic limit-cycle oscillations. One of its distinctive features is the existence of strong higher harmonics, which results in the time-averaged mean flow being strongly linearly unstable. For this class of flows, a simplified formalism capable of unravelling how exactly the instability grows and saturates is lacking. This study builds on previous work by Mantic-Lugo et al. (Phys. Rev. Lett., vol. 113, 2014, 084501) to fill in the gap using a parametrized approximation of an instantaneous, phase-averaged mean flow, coupled in a quasi-static manner to multiple linear harmonic disturbances interacting nonlinearly with one another and feeding back on the mean flow via their Reynolds stresses. This provides a self-consistent modelling of the mean flow-fluctuation interaction, in the sense that all perturbation structures are those whose Reynolds stresses force the mean flow in such a way that the mean flow generates exactly the aforementioned perturbations. The first harmonic is sought as the superposition of two components, a linear component generated by the instability and aligned along the leading eigenmode of the mean flow, and a nonlinear orthogonal component generated by the higher harmonics, which progressively distorts the linear growth rate and eigenfrequency of the eigenmode. Saturation occurs when the growth rate of the first harmonic is zero, at which point the stabilizing effect of the second harmonic balances exactly the linear instability of the eigenmode. The model does not require any input from numerical or experimental data, and accurately predicts the transient development and the saturation of the instability, as established from comparison to time and phase averages of direct numerical simulation data.