Wide-band gap SiC is a promising semiconductor material for microelectronic applications due to its superior electronic properties, high thermal conductivity, chemical and radiation stability, and extremely high break-down voltage. Over the past several years, tremendous advances have been made in SiC crystal growth technology. Nevertheless, SiC will not reach its anticipated potential until a variety of problems are solved, one of the problem is step bunching during step flow growth of SiC, because it could lead to uneven distribution of impurity and less smooth surfaces. In this paper, step bunching morphologies on vicinal 4H-SiC (0001) surfaces with the miscut toward $\left[ {1\bar 100} \right]$ or $\left[ {11\bar 20} \right]$ directions are studied with a three-dimensional kinetic Monte Carlo model, and then compared with the analytic model based on the theory of Burton-Cabera-Frank. In the kinetic Monte Carlo model, based on the crystal lattice of 4H-SiC, a lattice mesh is established to fix the positions of atoms and bond partners. The events considered in the model are adsorption and diffusion of adatoms on the terraces, attachment, detachment and interlayer transport of adatoms at the step edges. The effects of Ehrlich-Schwoebel barriers at downward step edges and inverse Schwoebel barrier at upwards step edges are also considered. In addition, to obtain more elaborate information about the behavior of atoms in the crystal surface, silicon and carbon atoms are treated as the minimal diffusing species. Finally, the periodic boundary conditions are applied to the lateral direction while the " helicoidal boundary conditions” are used in the direction of crystal growth. The simulation results show that four bilayer-height steps are formed on the vicinal 4H-SiC (0001) surfaces with the miscut toward $\left[ {1\bar 100} \right]$ direction, while along the $\left[ {11\bar 20} \right]$ direction, only bunches with two-bilayer-height are formed. Moreover, zigzag shaped edges are observed for 4H-SiC (0001) vicinal surfaces with the miscut toward $\left[ {11\bar 20} \right]$ direction. The formation of these step bunching morphologies on vicinal surfaces with different miscut directions are related to the extra energy and step barrier. The different extra energy for each bilayer plane results in step bunches with two-bilayer-height on the vicinal 4H-SiC (0001) surface. And the step barriers finally lead to the formation of step bunches with four-bilayer-height. Finally, the formation mechanism of the stepped morphology is also analyzed by a one-dimensional Burton-Cabera-Frank analytic model. In the model, the parameters are corresponding to those used in the kinetic Monte Carlo model, and then solved numerically. The evolution characteristic of step bunching calculated by the Burton-Cabera-Frank model is consistent with the results obtained by the kinetic Monte Carlo simulation.