A coupled treatment of transport processes, phase changes and mechanical settling is the core of any detailed snowpack model. A key concept underlying the majority of these models is the notion of layers as deforming material elements that carry the information on their physical state. Thereby an explicit numerical solution of the ice mass continuity equation can be circumvented, although with the downside of virtual no flexibility in implementing different coupling schemes for densification, phase changes and transport. As a remedy we consistently recast the numerical core of a snowpack model into an extendable Eulerian-Lagrangian framework for solving the coupled non-linear processes. In the proposed scheme, we explicitly solve the most general form of the ice mass balance using the method of characteristics, a Lagrangian method. The underlying coordinate transformation is employed to state a finite-difference formulation for the superimposed (vapor and heat) transport equations which are treated in their Eulerian form on a moving, spatially non-uniform grid that includes the snow surface as a free upper boundary. This formulation allows us to unify the different existing viewpoints of densification in snow or firn models in a flexible way and yields a stable coupling of the advection-dominated mechanical settling with the remaining equations. The flexibility of the scheme is demonstrated within several numerical experiments using a modular solver strategy. We focus on emerging heterogeneities in (two-layer) snowpacks, the coupling of (solid-vapor) phase changes with settling at layer interfaces and the impact of switching to a non-linear mechanical constitutive law. Lastly, we discuss the potential of the scheme for extensions like a dynamical equation for the surface mass balance or the coupling to liquid water flow.