High frequency (HF) radio waves propagating through regions of ionospheric plasma irregularities exhibit random fluctuations in both amplitude and phase referred to as scintillations, which heavily influences the detection performance of HF sky‐wave radar and hybrid sky‐surface wave radar systems. In this paper, a numerical multiple phase screen model for HF sky‐wave propagation is developed, which can simulate the HF wave field scintillations caused by the irregularities in inhomogeneous media. In this model, the random phase changes due to irregularities and the diffraction effect due to phase mixing are considered separately and computed sequentially, where the phase change is numerically calculated using a geometric optical phase delay to allow for the consideration of horizontal electron density gradients. Utilizing the proposed model, we first simulate the observations of scintillation from the direct wave of the HF hybrid sky‐surface wave radar. The statistics of simulation and measured data, including power spectrum and scintillation indices on amplitude and phase, are considered to verify the effectiveness of the proposed model. Then, we analyze the effects of traveling ionospheric disturbances (TIDs) on the intensity of scintillation. The simulation indicates that the horizontal electron density gradients due to TIDs cause the fluctuation of the amplitude scintillation index S4.