Using numerical simulations, we investigate the near-field temporal variability of axisymmetric helium plumes as a function of inlet-based Richardson (
${Ri}_0$
) and Reynolds (
${Re}_0$
) numbers. Previous studies have shown that
${Ri}_0$
plays a leading-order role in determining the frequency at which large-scale vortices are produced (commonly called the ‘puffing’ frequency). By contrast,
${Re}_0$
dictates the strength of localized gradients, which are important during the transition from laminar to turbulent flow. The simulations presented here span a range of
${Ri}_0$
and
${Re}_0$
, and use adaptive mesh refinement to achieve high spatial resolutions. We find that as
${Re}_0$
increases for a given
${Ri}_0$
, the puffing motion undergoes a transition at a critical
${Re}_0$
, marking the onset of chaotic dynamics. Moreover, the critical
${Re}_0$
decreases as
${Ri}_0$
increases. When the puffing instability is non-chaotic, time series of different variables are well-correlated, exhibiting only modest changes in the dynamics (including period doubling and flapping). Once the flow becomes chaotic, denser ambient fluid penetrates the core of the plume, similar to penetrating ‘spikes’ formed by Rayleigh–Taylor instabilities, leading to only moderately correlated flow variables. These changes result in a non-trivial dependence of the puffing frequency on
${Re}_0$
. Specifically, at sufficiently low
${Re}_0$
, the puffing frequency falls below the prediction from Wimer et al. (J. Fluid Mech., vol. 895, 2020). As
${Re}_0$
increases beyond the critical
${Re}_0$
, the puffing frequency increases and then drops back down to the predicted scaling. The dependence of the puffing frequency on
${Re}_0$
provides a possible explanation for previously observed changes in the scaling of the puffing frequency for high
${Ri}_0$
.