In this work, we investigate the effects of subsurface grain morphology on the microplasticity on the free surface of an interstitial-free steel alloy and, a nickel-based alloy under tensile deformation. We propose a method to generate volume elements (VE) for crystal plasticity (CP-FFT), using a stochastic Monte Carlo Markov Chain (MCMC) to vary the subsurface microstructure. The approach utilises weighted Voronoi tessellations, in which the subsurface grains are incrementally perturbed whilst ensuring the surface remains invariant. This method allows for sampling subsurface microstructures with large variance for a given surface observation. We demonstrate that the surface microplasticity is significantly influenced by the subsurface microstructure, as previously reported, however we further reveal that small perturbations in the grain morphology gives rise to noticeable variation in the evolution of deformation on the surface. This approach demonstrates the need to consider a distribution of predictions from simulations of several equivalent volume elements. Using this statistical approach, it becomes possible to identify of regions of deformation that are least affected by subsurface uncertainty, which can be used to inform model calibration. By direct comparison with experimental results obtained using high-resolution digital image correlation (HRDIC), we can discriminate the uncertainty from subsurface effects from the inherent uncertainty present in the constitutive model, thereby allowing for the discovery of deformation mechanisms that are presently not considered in the material model.