Forward modeling of ground penetrating radar (GPR) is an important part to the inversion/modeling of the observed data. The aim of this study is to establish specific numerical schemes for forward modeling of GPR data by finite difference frequency domain (FDFD) method which were originally developed for seismic or finite difference time domain (FDTD) method. A total number of six modified and improved FDFD techniques have been used to discretize the two-dimensional (2D) transverse electric (TE)-mode scalar wave equation in order to find the suitable method for this. These techniques include five-point classical to nine-point mixed unstaggered-grid configurations. The numerical schemes for three unsplit perfectly matched layer (PML) for nine-point mixed unstaggered-grid configurations are also presented. The applicability of these techniques is tested by using the underground models of relative permittivity and conductivity for the two cases of homogeneous and 2-cross models. GPR shot gather data for these two models are also produced for this study. The relative reflection errors of the numerical schemes are also estimated for the homogeneous model to comprehend the appropriate method for the modeling. The algorithm for complex-frequency shifted PML (CFSPML) gives the least error in case of the forward modeling of the GPR data.