Miscible recovery techniques are among the most efficient and widely used enhanced oil recovery methods both in conventional and in heavy oil reservoirs (miscible gas flooding and solvent injection, respectively). The efficiency of a miscible flood is significantly dependent on the degree of mixing that occurs between the fluids. Mixing of the injected gas/solvent with oil leads to solvent dilution and thus reduces the effective strength of miscible displacement. The primary mechanisms of miscible mixing in pore space are molecular diffusion and convection (mechanical spreading) due to bulk flow velocity. The combined effect of these two mechanisms on the degree of mixing can be characterized by the dispersion coefficient. Laboratory measurement of this parameter is difficult due to the fact that it is strongly dependent on flow and porous media properties. An alternative approach is to numerically calculate dispersion coefficient using pore scale digital core analysis. In this study, pore level miscible displacements in heterogeneous porous media are computationally modeled through simultaneous solving of NavierStokes, Continuity, and Convection-Diffusion equations on virtual unconsolidated granular porous media. These virtual media are constructed by a pattern generator using the concept of random packing of grains. The heterogeneity level of these media is characterized by a coefficient of variations defined as the ratio of standard deviation in grain diameter to the mean grain diameter. Using the results of numerical simulations, longitudinal (along the flow direction) dispersion coefficient is calculated at different values of flow velocity, viscosity ratio, and pore scale medium heterogeneity. The input diffusion coefficients are either constant or functions of concentration. The results of the simulations show that when viscosity ratio is unity, at very low Peclet numbers the dominant transport mechanism is molecular diffusion. At high Peclet numbers, however, advection is dominant and the magnitude of longitudinal dispersion coefficient scales with Np^1.2. Higher heterogeneity of media will results in higher values of longitudinal dispersion coefficient. Also, it is found that using an appropriate average diffusion coefficient over the range of concertation interval can adequately account for concentration dependency of diffusion coefficient. The effect of viscosity contrast on longitudinal dispersion coefficient is shown to be significant. Viscous fingering due to an unfavorable viscosity ratio leads to higher values for longitudinal dispersion coefficient.