We describe the $10$-dimensional space of $Sp(2)$-invariant $G_2$-structures on the homogeneous $7$-sphere $S^7=Sp(2)/Sp(1)$ as $mathbb{R}^+times Gl^+(3,mathbb{R})$. In those terms, we formulate a general Ansatz for $G_2$-structures, which realises representatives in each of the $7$ possible isometric classes of homogeneous $G_2$-structures. Moreover, the well-known nearly parallel round and squashed metrics occur naturally as opposite poles in an $S^3$-family, the equator of which is a new $S^2$-family of coclosed $G_2$-structures satisfying the harmonicity condition $div T=0$. We show general existence of harmonic representatives of $G_2$-structures in each isometric class through explicit solutions of the associated flow and describe the qualitative behaviour of the flow. We study the stability of the Dirichlet gradient flow near these critical points, showing explicit examples of degenerate and nondegenerate local maxima and minima, at various regimes of the general Ansatz. Finally, for metrics outside of the Ansatz, we identify families of harmonic $G_2$-structures, prove long-time existence of the flow and study the stability properties of some well-chosen examples.