Snow is a complex material. It resists elastic normal and shear deformations, while some deformations are plastic. Snow can deform and break. It can be significantly compressed and gets harder under compression. Existing snow solvers produce impressive results. E.g., hybrid Lagrangian/Eulerian techniques have been used to capture all material properties of snow. The auxiliary grid, however, makes it challenging to handle small volumes. In particular, snow fall and accumulation on surfaces have not been demonstrated with these solvers yet. Existing particle-based snow solvers, on the other hand, can naturally handle small snow volumes. However, existing solutions consider simplified material properties. In particular, shear deformation and the hardening effect are typically omitted. We present a novel Lagrangian snow approach based on Smoothed Particle Hydrodynamics (SPH). Snow is modeled as an elastoplastic continuous material that captures all above-mentioned effects. The compression of snow is handled by a novel compressible pressure solver, where the typically employed state equation is replaced by an implicit formulation. Acceleration due to shear stress is computed using a second implicit formulation. The linear solvers of the two implicit formulations for accelerations due to shear and normal stress are realized with matrix-free implementations. Using implicit formulations and solving them with matrix-free solvers allows to couple the snow to other phases and is beneficial to the stability and the time step size, i.e., performance of the approach. Solid boundaries are represented with particles and a novel implicit formulation is used to handle friction at solid boundaries. We show that our approach can simulate accumulation, deformation, breaking, compression and hardening of snow. Furthermore, we demonstrate two-way coupling with rigid bodies, interaction with incompressible and highly viscous fluids and phase change from fluid to snow.