This paper presents the design and numerical simulation of vertical-cavity surface-emitting laser (VCSEL) incorporating a high-contrast grating (HCG) by using a three-dimensional (3-D) finite difference time domain (FDTD) method. The polarization dependence of HCG parameters as well as HCG reflector properties are numerically studied. We investigate the resonant spectra behavior and modal properties with respect to different aperture size. Furthermore, the quality factor (Q factor) related to modal loss is calculated, and is used to explain the mode selection of the HCG-VCSELs. Based on the simulation and analysis, it is found that the designed HCG-VCSELs exhibit larger mode discrimination compared to the conventional VCSELs. Our work could provide guidelines for designing and optimizing polarization-stable, single-mode VCSELs for use in chip-scale Cs atomic clocks.