In many geoscience applications, the data extracted from environmental variables are very limited. Multiple-point geostatistical (MPS) approaches simulate these variables and associated uncertainties at unknown locations by using an exemplar model for the field, called the training image (TI). Existing MPS approaches aim at simulating the field in a way consistent with both available conditional data and TI properties. The inevitably limited size of the training database usually leads to an underestimated variability between different realizations as compared to the variability of the real phenomenon. Furthermore, in over-conditioned regions, patch-based methods often tend to paste the same patch in all realizations. In this paper, we suggest an optimization-based approach for MPS simulation that simulates a bunch of realizations simultaneously. In addition to maintaining consistency with both conditional data and TI properties, the proposed method aims at maximizing the variability between different realizations. Our experiments show that the proposed strategy enhances the variability of the realizations to better conform with real variabilities. The idea of targeting variance maximization can potentially be applied to other MPS simulation methods by simulating a bunch of realizations simultaneously with a constraint to avoid similar patterns at the same location in different realizations.