The evaluation of extreme snowfalls is an important challenge for hazard management in mountainous regions. In this paper, the extreme snowfall data acquired from 40 meteorological stations in the French Alps since 1966 are analyzed using spatial extreme statistics. They are then modeled within the formal framework of max-stable processes (MSPs) which are the generalization of the univariate extreme value theory to the spatial multivariate case. The three main MSPs now available are compared using composite likelihood maximization, and the most flexible Brown-Resnick one is retained on the basis of the Takeuchi information criterion, taking into account anisotropy by space transformation. Furthermore, different models with smooth trends (linear and splines) for the spatial evolution of the generalized extreme value (GEV) parameters are tested to allow snowfall maps for different return periods to be produced. After altitudinal correction that separates spatial and orographic effects, the different spatial models are fitted to the data within the max-stable framework, allowing inference of the GEV margins and the extremal dependence simultaneously. Finally, a nested model selection procedure is employed to select the best linear and spline models. Results show that the best linear model produces reasonable quantile maps (assessed by cross-validation using other stations), but it is outperformed by the best spline model which better captures the complex evolution of GEV parameters with space. For a given return period and at fixed elevation of 2000 m, extreme 3 day snowfalls are higher in the NE and SE of the French Alps. Maxima of the location parameter of the GEV margins are located in the north and south, while maxima of the scale parameter are located in the SE which corresponds to the Mediterranean influence that tends to bring more variability. Besides, the dependence of extreme snowfalls is shown to be stronger on the local orientation of the Alps, an important result for meteorological variables confirming previous studies. Computations are performed for different accumulation durations which enable obtaining magnitude-frequency curves and show that the intensity of the extremal directional dependence effect is all the more important when the duration is short. Finally, we show how the fitted model can be used to evaluate joint exceedence probabilities and conditional return level maps, which can be useful for operational risk management.