<p>The determination of constitutive modeling is an important theoretical basis for solving geotechnical engineering problems. On the other hand, soil heterogeneity affects the prediction of constitutive models. In this study, the simulation of the SANISAND constitutive model, for the single-element triaxial test was coded in MATLAB. In the next step, sensitivity analysis was performed to evaluate the most effective model parameter, and stochastic analysis of the constitutive model was done based on these parameters using the Monte Carlo Simulation (MCS). In another part of the study, an evaluation of the stochastic static liquefaction potential of sandy soils and a comparative study between several constitutive models were presented. The results of the stochastic simulation of the model show that the soil heterogeneity has the greatest effect on the pore water pressure in undrained conditions, among the deviatoric stress and pore water pressure. While, among the deviatoric stress, volume strain, and void ratio, soil heterogeneity creates more uncertainty in volumetric strain in drained conditions. The results showed that as the relative density decreases, the difference between the simulation and experimental results increases. Stochastic static liquefaction analysis shows that, with an increase in confining pressure, the uncertainty of the deviatoric stress increases.</p>