Abstract. Increasing evidence is becoming available about not only the surface differential rotation of rapidly rotating cool stars but, in a small number of cases, also about temporal variations, which possibly are analogous to the solar torsional oscillations. Given the present difficulties in resolving the precise nature of such variations, due to both the short length and poor resolution of the available data, theoretical input is vital to help assess the modes of behaviour that might be expected, and will facilitate interpretation of the observations. Here we take a first step in this direction by studying the variations in the convection zones of such stars, using a two dimensional axisymmetric mean field dynamo model operating in a spherical shell in which the only nonlinearity is the action of the azimuthal component of the Lorentz force of the dynamo generated magnetic field on the stellar angular velocity. We consider three families of models with different depths of dynamo-active regions. For moderately supercritical dynamo numbers we find torsional oscillations that penetrate all the way down to the bottom of the convection zones, similar to the case of the Sun. For larger dynamo numbers we find fragmentation in some cases and sometimes there are other dynamical modes of behaviour, including quasi-periodicity and chaos. We find that the largest deviations in the angular velocity distribution caused by the Lorentz force are of the order of few percent, implying that the original assumed "background" rotation field is not strongly distorted.