Modeling and assessing spatial uncertainty of mineral deposits is critical for capital investments in mining projects. However, traditional approaches for modeling geological domains and geostatistical estimation provide smoothed representations of the deposit and ignore spatial variability, hence mislead downstream decisions. Spatial variability and related uncertainty in modelling mineral deposits can be quantitatively described by stochastic spatial simulations. This is demonstrated through an application to the Puma deposit, a major nickel lateritic asset in Brazil. To integrate the variability of the regolith profiles of the deposit, their thicknesses are calculated after an unwrinkling process is applied and the deposit is then jointly simulated using min/max autocorrelation factors (MAF). The realizations serve as geological boundaries within which Ni, Co, Fe, SiO2, MgO and Dry-tonnage factor (DTF) are subsequently jointly simulated directly at block support scales. The final result is a series of equiprobable representations of the deposit that incorporate both grade and tonnage uncertainty. These simulations can be used to assess the uncertainty about key aspects of the project, such as the strict control of the ore's quality that feeds the ferronickel processing plant. The framework explored herein takes advantage of the MAF and direct block simulation approaches which facilitate the joint simulation of large multivariate deposits for industrial environments.