Numerical experiments performed with all possible combinations of approximations for the equations of a onedimensional planetary boundary layer model employing a one-equation turbulence closure scheme and a staggered, homogeneous, ยฎnite difference grid show that spikes and jittering in the vertical proยฎles of predicted variables are caused by an inappropriate numerical description of the turbulent kinetic energy equation. Jittering arises when using both moderate and large time steps, while spikes occur for large time steps only. There exist two types of jittering: transient and permanent. The former has a lifetime of about 1 week and occurs with large time steps. Permanent jittering, on the other hand, is observed when using moderate time steps and has a lifetime of more than a 1000 days. Introducing an iterative procedure for the eddy viscosity eliminates spikes and permanent jittering but is unsuccessful in removing transient jittering. It is further found that the transition to instability of the solution can be either a sudden or a gradual one. In the case of a sudden transition no numerical peculiarities are observed, whereas for a gradual transition to instability, jittering is present during the entire transition period. In this sense, jittering may be regarded as a herald of instability. Finally, a combination which employs implicit approximations consistently for all the terms in the model equations, regardless of whether using an implicit or semi-implicit approximation for the Coriolis term, proves to be devoid of any numerical artefacts without the need for introducing the iterative procedure for the eddy viscosity.