  We study the time-evolution of axisymmetric oscillations of superfluid
magnetars with a poloidal magnetic field and an elastic crust, working in
Newtonian gravity. Extending earlier models, we study the effects of
composition gradients and entrainment on the magneto-elastic wave spectrum and
on the potential identification of the observed quasi periodic oscillations
(QPOs). We use two-fluid polytropic equations of state to construct our stellar
models, which mimic realistic composition gradient configurations. The basic
features of the axial axisymmetric spectrum of normal fluid stars are
reproduced by our results and in addition we find several magneto-elastic waves
with a mixed character. In the core, these oscillations mimic the shear mode
pattern of the crust as a result of the strong dynamical coupling between these
two regions. Incorporating the most recent entrainment configurations in our
models, we find that they have a double effect on the spectrum: the magnetic
oscillations of the core have a frequency enhancement, while the mixed
magneto-elastic waves originating in the crust are moved towards the
frequencies of the single-fluid case. The distribution of lower-frequency
magneto-elastic oscillations for our models is qualitatively similar to the
observed magnetar QPOs with $\nu < 155 $Hz. \emph{In particular, some of these
QPOs could represent mixed magneto-elastic oscillations with frequencies not
greatly different from the crustal modes of an unmagnetised star.} We find that
many QPOs could even be accounted for using a model with a relatively weak
polar field of $B_{p} \simeq 3\cdot 10^{14}$G, because of the superfluid
enhancement of magnetic oscillations. Finally, we discuss the possible
identification of 625 and 1837~Hz QPOs either with non-axisymmetric modes or
with high-frequency axisymmetric QPOs excited by crustal mode overtones.
