## 3 Model Development
### Updated Snow Parametrization  
A multi-layer snow scheme (ML) is introduced in ECLand, replacing the previous single-layer (SL) snow scheme. ML is an intermediate complexity snow scheme following the definition by (2017), representing the vertical structure and prognostic evolution of snow temperature, density, liquid water content and albedo. A schematic of ML is shown in Figure 5. A detailed description of ML is reported in [PERSON] et al. (2018) and in the following only the main characteristcs of ML are reported. The heat budget of the snowpack is solved in two-steps. At first, a preliminary snow temperature is computed solving the heat budget in absence of melting, considering the vertical diffusion of heat and the absorption of solar radiation by each layer. For the latter, a parametrization adapted from [PERSON] (2017) is used, in which the penetration of solar radiation decreases exponentially with depth and the extinction coefficient is a function of snow density. Then, melting and freezing processes are diagnosed by computing the heat content of each snow layer, also considering the liquid water content of each snow layer, and the snow temperature is modified accordingly to the latent heat released/absorbed. The movement of liquid water between layers is taken into account using a bucket-type formulation ([PERSON] et al., 2017). The heat flux at the snow-soil interface is computed as the sum of the residual heat flux at the end of the heat budget computation and the basal heat flux (\(G_{B}\) in Figure 5). The number of active layers (\(N\)) varies depending on snow depth \(D_{sn}\), from a minimum number of 1 layer to a maximum of \(N_{max}=5\). Multiple layers are used when the grid box is completely snow-covered, that is for \(D_{sn}\geq 0.1\) m. As opposed to [PERSON] et al. (2018), a different discretization algorithm is applied over flat terrain and complex terrain regions, the latter defined as regions with the standard deviation of the sub-grid-scale orography greater than 50 m. A variable vertical discretization is used to maintain a relatively high vertical resolution for the snow layers responding to fast time scales also for a deep snowpack. Over flat terrain, the thickness of the topmost layer in contact with the atmosphere is fixed to 0.05 m; the second and third topmost layers, as well as the one in contact with the soil underneath, can increase up to a certain maximum \(D_{i}^{max}\) as described in (2017); for \(D_{sn}>0.45\) m, the layer \(N_{max}-1\) is used as an accumulation layer. An example of the vertical discretization of a 1.0 m-thick snowpack is shown in Figure 5b. Over complex terrains, for \(D_{sn}<0.25\) m the same vertical discretization as described for a flat terrain is used. For a snowpack with \(D_{sn}\geq 0.25\) m, the minimum (\(D_{min,i}\)) and maximum (\(D_{max,i}\)) thicknesses of each layer (i.e., the minimum and maximum allowed thickness for an active snow layer) varies with snow depth as follows:  
\[D_{min,i} =\min(0.25,0.10+\alpha_{o}\left(D_{sn}-0.25\right))\quad\text{if}\ i=1 \ldots,N_{max}\] \[D_{max,i} = \begin{cases}\min(0.25,0.10+\alpha_{o}\left(D_{sn}-0.25\right)) \quad\text{if}\ i=1\\ \min(0.30,0.15+\alpha_{o}\left(D_{sn}-0.25\right))\quad\text{if}\ i=2\ldots,N_{ max},\end{cases} \]  
with \(\alpha_{o}\) set to 0.10. This implies that for a thick snowpack over complex terrain the snow layers are thicker than a snowpack with same depth over a flat region. A full description of the vertical discretization algorithm is detailed in [PERSON] et al. [41]. A number of snow physical processes are also changed in ML compared to SL. The snow heat conductivity is parametrized following [PERSON] et al. [43], with an additional term considering water vapor diffusion effects [44]. Changes in snow density due to wind blowing are an important process in particular at high latitude [45; 46]. In ML a parametrization of wind-driven compaction is introduced, using the formulation proposed by [PERSON] et al. [46].