The effective mechanical behavior of snow can be deduced from microstructural homogenization through numerical simulations. Although such numerical upscaling of elasticity and strength of snow microstructures is standard (using FEM), numerical schemes to study generic features of the transition from small to large strain situations that involve yielding and failure are scarce. This prevents the development of accurate homogenized constitutive models valid for the post-failure and large deformation regimes. It has been shown that treating this transition is feasible using DEM under the assumption of particulate microstructures. However, this requires snow microstructures to be segmented into a granular collection of (usually spherical) cohesive elements. Here, we suggest generating random porous microstructures by level-cutting Gaussian random fields and using the material point method to numerically simulate them under mechanical loading. This allows investigating both small and large deformation characteristics of irregular porous media, such as snow, where a segmentation into grains and bonds can be ambiguous. We demonstrate our approach by examining elasticity and failure as a function of a wide range of solid volume fractions, from 20% (low-density snow) to 80% (high-density firn), as the most important control on the mechanical behavior. Observing that onset of failure can be well described through the second order work, we show that the failure strength follows a power law similar to that of the elastic moduli. Moreover, we propose that the failure envelope can be approximated by a porosity-dependent quadratic curve in the space of the two first stress invariants. Furthermore, we observe that plastic deformation appears to be governed by an associative plastic flow rule. Finally, these results combined with a viscoplastic Perzyna model and a sintering (hardening) model should allow us to develop a universal homogenized snow constitutive model.