Density irregularities have been observed in subauroral polarization streams (SAPS). One hypothesis of the cause of this ionospheric turbulence, based on the background morphology, is the gradient drift instability (GDI). This work models the GDI using a 2D electrostatic fluid model to determine if it is a viable cause of turbulence generation in SAPS. The model solves a perturbed set of continuity, energy, and current closure equations using a pseudo-spectral method. A statistical study of different velocity profiles, based on SuperDARN radar and GPS total electron content data, is used to prescribe parameters in the numerical model. The parameter space of different SAPS profiles is explored to study the effect on GDI development. As the velocity shear is initialized closer to the unstable density gradient, the GDI becomes increasingly damped. For these cases, the density and electric potential turbulence cascades obtained from the numerical model follow power laws of about -5/3 or -2, which is in agreement with observational data. If the velocity shear significantly overlaps the unstable density gradient, the GDI becomes stabilized. Decreasing the velocity gradient scale length can cause instabilities that grow inside SAPS which have turbulence cascades with power laws of -6 for the density and -8 for the electric potential. In all parameter regimes explored, the instabilities are unable to propagate through the velocity shear. Turbulence is generated for a variety of SAPS relevant conditions; therefore, the GDI has been shown to be a viable candidate for generating ionospheric irregularities in SAPS.