Review article: Parameterizations of snow-related physical processes in land surface models

Snow on land surface plays a vital role in the interaction between land and atmosphere in the state-of-the-art land surface models (LSMs) and the real world. Since the snow cover affects the snow albedo and the ground and soil heat fluxes, it is crucial to detect snow cover changes accurately. It is challenging to acquire observation data for snow cover, snow albedo, and snow depth; thus, an excellent alternative is to use the simulation data produced by the LSMs that calculate the snow5 related physical processes. The LSMs show significant differences in the complexities of the snow parameterizations in terms of variables and processes considered. Thus, the synthetic intercomparisons of the snow physics in the LSMs will help the improvement of each LSM. This study revealed and discussed the differences in the parameterizations among LSMs related to snow cover fraction, snow albedo, and snow density. We selected the most popular and well-documented LSMs embedded in the Earth System Model or operational forecasting systems. We examined single layer schemes, including the Unified Noah 10 Land Surface Model (Noah LSM), the Hydrology Tiled ECMWF Scheme of Surface Exchanges over Land (HTESSEL), the Biosphere-Atmosphere Transfer Scheme (BATS), the Canadian Land Surface Scheme (CLASS), and multilayer schemes of intermediate complexity including the Community Noah Land Surface Model with Multi-Parameterization Options (NoahMP), the Community Land Model version 5 (CLM 5), the Joint UK Land Environment Simulator (JULES), and the Interaction Soil-Biosphere-Atmosphere (ISBA). First, we identified that BATS, Noah-MP, JULES, and ISBA reflect the snow depth and 15 roughness length to parameterize snow cover fraction, and CLM 5 accounts for the standard deviation of the elevation value for the snow cover decay function. Second, CLM 5 and BATS are relatively complex, so that they explicitly take into account the solar zenith angle, black carbon, mineral dust, organic carbon, and ice grain size for the determinations of snow albedo. Besides, JULES and ISBA are also complicated model which concerns ice grain size, solar zenith angle, new snow depth, fresh snowfall rate, and surface temperature for the albedo scheme. Third, HTESSEL, CLM 5, and ISBA considered the effects of 20 both wind and temperature in the determinations of the new snow density. Especially, ISBA and JULES considered internal snow characteristics such as snow viscosity, snow temperature, and vertical stress for parameterizing new snow density. The future outlook discussed geomorphic and vegetation-related variables for the further improvement of the LSMs. Previous studies clearly show that spatio-temporal variation of snow is due to the influence of altitude, slope, and vegetation condition. Therefore, we recommended applying geomorphic and vegetation factors such as elevation, slope, time-varying roughness 25 1 https://doi.org/10.5194/tc-2021-319 Preprint. Discussion started: 30 November 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Physical processes related to snow play an essential role in interacting with the heat and moisture flux between the atmosphere and the land surface and between the land surface and soil. The snow's spatial and temporal heterogeneity causes difficulties 30 in observation, understanding, and modeling (Webster et al., 2018). Changes in snow cover and snow albedo are challenging to measure directly from the ground observation system. Moreover, snow depth is not easy to measure spatially. Therefore, Snow cover plays a vital role in surface albedo, which leads to the exchange of energy and water flux between land and atmosphere. If the snow cover fractions are estimated lower, the land surface albedo will be underestimated (Roesch et al., 2001). In the spring, when solar radiation is at its most intense, snow cover significantly affects the global radiation balance, and 40 snow cover change significantly impacts nature and human systems. Snow is a crucial water source in rivers and groundwater, and it is a more critical resource especially in semi-arid regions (Armstrong and Brun, 2008). Since snow cover is valuable for use as a leisure and travel resource in mid-latitude mountainous regions, it is crucial to monitor snow cover changes accurately. Besides, snow albedo is an important physical property that causes energy interactions between the land surface and the atmosphere, along with the thermal barrier properties of snow and the possibility of phase changes due to latent heat 45 (Armstrong and Brun, 2008). In particular, albedo-the reflectivity of the ground during incident solar radiation-is vital for the land surface's energy balance (Warren and Wiscombe, 1981). Wiscombe and Warren (1980) and Warren and Wiscombe (1980) created early models of snow albedo, and Zhong et al. (2017) improved the snow albedo simulation. Since snow depth is a concept that directly expresses the amount of snowfall, it is imperative to predict or observe accurately to prevent accidents caused by snow. Snow causes delays for both ground and air traffic and heavy snow disasters for crops and livestock. In 50 addition, it causes accidents, injuries, and deaths due to snowfall or snow avalanches (Armstrong and Brun, 2008). Therefore, LSM's parameterization method to determine snow depth using other atmospheric variables as input variables can be essential information to prevent snow-related disasters in areas with no observed snow depth.
As a result of analyzing previous studies, there were studies related to model development or improvement for one or two LSM models. For example, there were snow cover studies (Brun et al., 1992;Yang et al., 1997;Roesch et al., 2001;Brown 55 et al., 2006;Roesch and Roeckner, 2006;Brun et al., 2008;Dutra et al., 2009;Ma et al., 2019), and snow albedo studies (Roesch and Roeckner, 2006;Mölders et al., 2007;Dutra et al., 2009;Wang and Zeng, 2010;Malik et al., 2014;Bartlett and Verseghy, 2015;Zhong et al., 2017;Wang et al., 2020), and studies related to snow depth (Gottlib, 1980;Longley, 1960;Dutra Parameterization Options (Noah-MP;Niu et al., 2011); and 3) the Biosphere-Atmosphere Transfer Scheme (BATS;Dickinson et al., 1993). The Noah LSM and Noah-MP are the current and future generations of Korean Integrated Model (KIM; Hong 95 et al., 2018) operational forecasting systems of the Republic of Korea, respectively. We need to select representative models to derive ideas for future model improvement and concretely compare and analyze for the parameterization of the snow physics process. we analyzed BATS model because albedo parameterization of the Common Land Model (CoLM; Li et al., 2017), JSBACH (Reick et al., 2021), and JSBACH3-PF (Ekici et al., 2014) as well as Noah-MP is adapted from BATS method (Dickinson et al., 1986). 100 Therefore, this study attempts to understand the variables and parameters considered in each land surface model's current situation by looking at the parameters of three significant snow-related variables-snow cover, snow albedo, and snow density-through eight LSMs: Noah LSM, HTESSEL, BATS, CLASS, Noah-MP, CLM 5, JULES, and ISBA. In addition, we intend to provide insights into variables that may be considered in different models through a comparison of parameterization methods. Sections 2, 3, and 4 compare and analyze eight LSMs in relation to the parameterization of snow cover, snow albedo, 105 and snow density, respectively. We focused on the new snow density in the snow density part because snow compaction due to destructive metamorphism or overburden, melting, and drifting snow is too complicated. Section 5 provides insights for LSM improvement by discussing the comparison results and significant variables in other previous studies and presents summaries in Section 6.

110
Model studies show that snow cover changes significantly influence the land surface energy balance and ground temperature (Thomas and Rowntree, 1992;Viterbo and Betts, 1999). Since the snow-covered ground albedo is much higher than the snowfree ground albedo, accurate prediction of the presence or absence of snow is crucial for understanding near-surface energy balance. Snow cover can change grassland albedo by approximately 3-4 and mountain albedo by 2-3 ( Thomas and Rowntree, 1992;Betts and Ball, 1997;Jin et al., 2002;Barlage et al., 2010). Armstrong and Brun (2008) presented snow depth, snow 115 density, and SWE as fundamental properties when describing snow cover. Table 1 expressed the variables considered in the eight LSMs, including the main variables such as snow depth, snow density, roughness length and SWE. All LSMs except BATS, JULES, ISBA used the SWE to parameterize snow cover. The snow cover parameterization of Noah LSM and CLM 5 considered the threshold of SWE. The Noah LSM's snow cover formula reflects the land surface type through the SWE threshold parameter. CLM 5 also considers the land surface type when using the SWE 120 threshold, and it especially takes into account sub-grid geomorphic variation, the main difference among LSMs. HTESSEL is the simplest of the eight LSMs and considers only the SWE and snow density. The CLASS model used the snow depth and the snow depth threshold instead of the SWE and the SWE threshold, and the snow depth is calculated from SWE and snow density. The BATS, JULES, ISBA and Noah-MP model are different from the others in that it utilizes snow depth and roughness length to parameterize snow cover. Besides, they used different values of roughness length according to the land 125 surface type. Table 2 presents the different formulations for snow cover fraction in the eight LSMs. The SWE, the SWE threshold value (W max ) when the snow cover is 1, and the distribution shape parameter (P s ) determine the snow cover parameterization of the Noah LSM (Eq. (a) in Table 2) (Koren et al., 1999). If the SWE exceeds W max , the snow cover is assumed to be 1. Anderson (1973) considered the range of P s as a value between 2 and 4 in the empirical equation for snow cover fraction, and the Noah 130 LSM applied a fixed value of 2.6. The primary empirical values-the W max and the P s -affect the land surface albedo because the snow cover fraction variation controls the snow albedo. The Noah LSM provides W max values according to the land surface types based on the modified International Geosphere-Biosphere Program (IGBP) MODIS Noah classification. The Noah LSM designate W max values as 0.08 m, 0.04 m, 0.025 m, 0.02 m, and 0.01m for forests, short vegetation, tundra, low vegetation areas such as snow, glacial areas, and water bodies, respectively. Wang and Zeng (2010) showed that the W max value for 135 grassland and forest tend to have similar simulation results to the observed values when adjusted to 0.01 m and 0.2 m instead of 0.04 m and 0.08 m for each vegetation type. Livneh et al. (2010) improved the W max value to 0.02 m for non-forest and 0.04 m for the forest, reflecting the study area's vegetation properties. Therefore, we need to find an appropriate W max value that can guide the snow cover fraction following its observation value depending on the land surface type.
HTESSEL parameterized the snow cover fraction using a formulation of SWE and low and high densities of snow, which 140 showed a different route for the snow cover parameterization. The existing HTESSEL snow cover was calculated only as an SWE function (Eq. (b)). However, the new scheme introduced snow density in addition to the SWE (Eq. (c)). The existing parameterization method was expressed as extreme values of both snow densities, while the modified parameterization method was improved to reflect the changing snow density. The new scheme is straightforward, but it reflects the variability of snow cover between the beginning of a cold season (low snow density) and the end of an ablation season (high snow density) (Dutra 145 et al., 2009). At the beginning of the cold season, it is easy to cover the entire grid even if there is a small SWE. On the other hand, since snow melting that creates non-snow patches at the snow depletion period, more SWE is required for snow cover to be 1 in one grid. HTESSEL introduced snow density to the simulation of seasonal snow cover variation, thus devising a snow cover parameterization scheme that could reflect the physical process of snow.
In BATS, the snow cover fraction is measured in terms of the snow depth, so it changes over time by subtracting the 150 sublimation rate and snow melting from the snowfall rate. BATS parameterized the snow cover fraction, divided into areas with and without vegetation (Dickinson et al., 1993). For land surfaces without vegetation, it was parameterized with 0.1 m and snow depth (Eq. (d)), while for land surface with vegetation, it was parameterized as a function of the snow depth and the roughness length multiplied by 10 (Eq. (e)). Yang et al. (1997) described the total fraction of the grid square covered by snow where σ f is green vegetation fraction, f s,g and f s,v are fraction of ground covered by snow and fraction of vegetation covered by snow, respectively. Eq. (d) from Table 2 clarified as where D sn is snow depth and z 0g is roughness length for bare soil and the value is designated as 0.01 m (Yang et al., 1997).

160
The formation of equation (2) is same as Eq. (e) in Table 2. Yang et al. (1997) upgraded the snow cover fraction scheme using different roughness lengths according to the land surface type and snow depth (Eq. (f)). This new scheme reflects the aspect from Jordan et al. (2008) that snow cover and albedo are related to the vegetation and roughness length. Yang et al. (1997) showed that the simulated snow cover fraction and surface albedo has a close agreement with the observed, especially in short vegetation. Figure 1 demonstrates the comparison of BATS 165 formulation (Eqs. (d), (e), (f)) for snow cover fraction in the case of short grass land cover. Eq. (f) may be appropriate for grasses and agricultural lands (Yang et al., 1997) because it illustrates that the higher the snow depth, the faster cover the short grass. BATS focused on the snow-ground process and did not distinguish between the soil temperature and the lowest layer of snow. Geothermal heat can reach the upper surface of the snow by conduction. Nevertheless, it was assumed to ignore snow melting at the bottom of the snowpack (Dickinson et al., 1993).

170
CLASS parameterized the snow cover fraction using the snow depth (D sn ) and the threshold value of snow depth (D sn,lim ) (Eq. (g)) (Verseghy, 2012). The threshold of the snow depth (D sn,lim ) applied to 0.1 m. If the snow depth (D sn ) is >0.1 m, the snow cover is considered 1. When the snow depth is <0.1 m, parameterization for snow cover fraction uses Eq. (g) in  (2003). In CLASS, the snow depth threshold was assumed to be the same value regardless of the land surface type, so there could be a limitation in that it cannot follow the variability of the snow cover fraction according to the land cover type.
Noah-MP parameterizes the snow cover fraction using roughness length, snow depth as well as snow density. Yang et al. (1997) proposed the Hyperbolic Tangent(tanh) formulation (Eq. (f)) because they need to reflect two stages of snow depth and snow albedo relationship: sharply increasing albedo stage and slowly increasing albedo stage per changing snow depth.

180
This formulation has a limitation of too fast covering snow, so Eq. (f) does not consider the seasonal variation of snow cover patterns. Niu and Yang (2007) perceived this limitation of Eq. (f) which is fit for shorter vegetation, and revised formulation Eq. (j) using a melting factor for snow cover fraction (F melt ). F melt (Eq. (k)) calculates using fresh snow density, 100 kg m −3 , and current snow density (Eq. (l)) with determining melting curves parameter M which can be adjustable.
CLM 5 parameterizes the accumulation and depletion of snow cover fraction separately (Lawrence et al., 2018). The param-185 eterization of the snow cover fraction during the snow accumulation phase (Eq. (m)) reflects the new snow amount (q snow ∆t), a constant (k accum ) whose default value is 0.1, and the snow cover fraction (f n s ) from the previous time step. The snow cover fraction parameterization during the snow ablation phase (Eq. (n)) shows similarity with Noah LSM in that SWE and SWE thresholds influence the depletion curve. Swenson and Lawrence (2012) developed the empirically derived expression Eqs.  Table 2 (Swenson and Lawrence, 2012;Lawrence et al., 2018). roughly demonstrates the results of Jordan et al. (2008) in that the topographic variation of the ground causes changes in the snow cover depth.
Snow cover fraction parameterization in JULES (Eq. (p)) is identical with BATS snow cover fraction for vegetation (Eq. (e)). JULES adopted only snow cover fraction using snow density and roughness length for bare soil (Best et al., 2011) and ignored the vegetation fraction for calculating the snow cover fraction. ISBA separates snow cover fraction parameterization for bare soil (Eq. (q)) and vegetation (Eq. (r)) similarly to BATS (Table 2). ISBA formulation of snow cover fraction for bare 200 soil is similar to CLASS using snow depth and snow depth thresholds. Besides, the ISBA formulation structure of snow cover fraction for vegetation is the same as BATS; the coefficient is different between ISBA and BATS as 2 and 10, respectively. The vegetation roughness length in ISBA designated as specific value depending on each vegetation type (Decharme et al., 2019).
In particular, ISBA used equation (1) with GVF for the average between the snow cover fraction in the vegetated (f s,v ) and non-vegetated (f s,g ) as a part of a grid cell concept similar to BATS (Decharme et al., 2016). Among the eight LSMs, only 205 BATS and ISBA consider GVF directly for snow cover parameterization, and BATS, Noah-MP, JULES, and ISBA account for roughness length.

Snow albedo parameterization in LSMs
Snow albedo is the rate of reflection of incident radiation on snowpack, which plays an essential role in the energy balance of the land surface over snow during the melting season (Warren and Wiscombe, 1980;Wiscombe and Warren, 1980). One of the 210 land surface model's main objectives is the accurate distribution of incident radiation into radiation, sensible heat, and latent heat fluxes. Net radiation is a concept that includes both upward & downward shortwave radiation and upward & downward longwave radiation in the surface energy balance. The surface albedo-the ratio of reflected solar radiation to the incident solar radiation-determines the energy transfer between the air and land surface. Typical albedo values are as high as 80-90 % on a non-melting snow surface (Armstrong and Brun, 2008) and are expressed as various values according to the state of snow 215 surface (fresh, dry snow; old, dry snow; wet snow; melting ice/snow) or the land surface type (snow-covered forest, snow-free vegetation/soil, water) ( Table 3).
One of the most rapid changes in land surface albedo is the presence or absence of snow cover. Noah LSM, CLM 5, and JULES use the snow cover fraction as a significant variable in the albedo parameterization (Table 4). Thus, they calculate snow albedo by summing the snow-free albedo and snow-covered albedo according to the snow cover ratio. In Jordan et al. 220 (2008), the interception of radiation by vegetation strongly changes the albedo in the snow-covered area. The vegetation type, vegetation density, snow accumulation on the canopy, and incident radiation can affect snow albedo alteration (Jordan et al., 2008). For instance, Noah LSM, HTESSEL, BATS, and CLM 5 consider the vegetation effects on snow albedo parameterization indirectly by using the maximum snow albedo for each vegetation type or applying the green vegetation fraction value. The albedo for each grid point can be calculated by linearly weighting the snow-free albedo and snow and vegetation albedo in 225 the snow-covered area. In particular, BATS, CLM 5, and JULES were relatively complex land surface models because they identify direct & diffuse albedo, visible & near-infrared albedo separately, taking into account all the effects of the solar zenith angle, and ice grain size. Additionally, BATS and CLM 5 considers dust and carbon. Table 5 presents snow albedo parameterization in eight different LSMs. In the Noah LSM, the albedo is parameterized by a weighted average of the snow-free albedo and the snow-covered albedo (Eq. (a) in Table 5). The snow albedo is closely 230 related to the snow cover fraction . The snow-covered patch can be distributed differently due to the slope or altitude, the spatial direction of a snowdrift, and the snow melting period. Moreover, the snow albedo on the surface can increase noticeably as the vegetation rate decreases. When there is short or large vegetation, the snow-covered area may exist in the form of a patch. The surface albedo is parameterized to account for this vegetation effect, as shown in Eq. (b), using the green vegetation fraction (σ f ) value. at the time of snowfall ablation season was set to 0.82-0.91, the optimal value was almost the same range at 0.83-0.91 (Sun et al., 2019). In addition, Malik et al. (2014) suggested improved values only for the time of snow accumulation, and A was in the range 0.9-0.94, and B was in the range 0.5-0.66 for different areas.
The original scheme of the snow albedo for HTESSEL follows the formulation of Baker et al. (1990), Verseghy (1991), Douville et al. (1995) and separates into melting and non-melting cases. Its formulation structure is the same as the past version 260 of ISBA parameterization (Mölders et al., 2007). If the snow does not melt, the albedo uses a linearly decaying coefficient (γ a = 0.008). When the snow melting exceeds 0 or when the snow surface temperature is > 2 K lower than the freezing point temperature, HTESSEL parameterized the albedo as Eq. (e) in Table 5 with the exponential decay factor (γ f = 0.24) and the minimum albedo (α min,H = 0.5). In other words, the higher the temperature, the more quickly the snow albedo decreases.
When the amount of new snowfall exceeds 1 kg m −2 hr −1 , the maximum snow albedo is set to 0.85. However, Pedersen and 265 Winther (2005) and Mölders et al. (2007) pointed out that the original scheme has a disadvantage of too low threshold value when used for numerical forecasting in snow albedo parameterization. The parameterization was improved as shown in Eq.
(f), so the snow albedo was continuously reset. HTESSEL assumed that 10 kg m −2 of new snowfall is required to reset the maximum snow albedo to 0.85 consistently.
The original scheme uses a fixed value of 0.15 for albedo under a canopy regarding the vegetation's influence. However, this 270 is insufficient to describe the variability of snow albedos according to the vegetation characteristics. As a result of observations in a forest area with snow intercepted by vegetation, the albedo appeared to be 0.3 immediately after heavy snowfall and changed to 0.2 after a few days according to the forest type (Betts and Ball, 1997;Viterbo and Betts, 1999). The snow intercepted by vegetation disappears due to wind blowing in cold temperatures or snow melting in warm temperatures. Therefore, HTESSEL set the albedo for each vegetation condition differently according to the 16 types of IGBP vegetation. HTESSEL 275 used the white-sky albedo climate values (2000)(2001)(2002)(2003)(2004) in the 0.3-5 µm band of the northern hemisphere with snow on the ground suggested by Moody et al. (2007).
BATS calculates the snow albedo by classifying the diffuse and direct radiation and parameterizes each albedo by separating it into visible and near-infrared regions. Concerning the diffuse radiation, when the solar zenith angle is less than 60°, the new snow albedo for visible radiation incident is 0.95, and the new snow albedo for near-infrared solar radiation is 0.65. The albedo 280 decay parameterization method applies over time (Dickinson et al., 1993) (Eqs. (g) and (h) in Table 5). As a factor associated with snow decay, BATS calculates the decreasing term F age through equations (3)-(11) (Wang et al., 2020): F age give the fractional reduction of snow albedo due to growth of snow grain size and accumulation of dirt and soot by where tauss t is a nondimensional age of snow, defined as with swe t representing the snow water equivalents of the current time step. The nondimensional age of snow sge t depends on a model prognostic variable δ a as follows where δ s is expressed as Here, swe t − swe t−1 is the change of snow water equivalent (in mm or kg m −2 ) in one time step ∆t, and swe c indicates the critical value of new snow water equivalent to fully cover the old snow, assumed to be 1 mm; and δ a is parameterized as In (7), Age 1 represents the effect of growing snow particle size by vapor diffusion which is formulated as where arg = 5000 1 with T f the freezing temperature (273.16 K) and T g the ground temperature (K); Age 2 reflects the effect of recrystallization of melting moisture or close to it and is represented as and Age 3 represents the effects of dust and soot as follows: BATS designated Age 3 as 0.01 for Antarctica and 0.3 for other areas (Dickinson et al., 1993).
Concerning the direct radiation, BATS parameterizes the new snow albedo using factors such as the diffuse snow albedo and 305 the solar zenith angle (Eqs. (i) and (j) in Table 5). F zenith is a factor to account for solar zenith angle impact, calculated as where b can be adjustable to best available data, and it is designated b = 2.0 in BATS (Dickinson et al., 1993). Equation (12) has the property for all b that it vanishes at cos(zen) = 0.5, so F zenith values equal to 0. The value of b is unity at cos(zen) = 0 which is the same as the solar zenith angle in the horizontal plane, so F zenith is equal to 1. In Wang et al. (2020), the method 310 of parameterizing F zenith was further simplified as follows: where F zenith is a value between 0 and 1; for example, when the solar zenith angle exceeds 60 degrees, the snow albedo for visibility area increases.
albedo is supposed to be 0.84 in the original CLASS. CLASS calculated the new snow albedo as the weighted average of 0.84 and the ground albedo before snowfall. Fresh snowfall albedo in CLASS (Eq. (k)) has the same structure as HTESSEL's new snowfall albedo parameterization (Eq. (f)). However, there is a difference in that CLASS assigned a value of 1 mm for the SWE threshold of 100% snow cover. In addition, CLASS assumed that the snow albedo in the current time step is over 0.55 (Eq. (l)).
Recently, Melton et al. (2020)  Melton and Arora, 2016). CLASSIC model parameterized snow albedo decay based on the data of Aguado (1985), Robinson and Kukla (1984), and Dirmhirn and Eaton (1975). Melton et al. (2019) showed that it is the same as empirical exponential decay functions applied to CLASS-CTEM by Verseghy (2017). The background old snow value was 0.55 in the original 325 CLASS (Eq. (l)). However, the background old snow value (α sn,total,old ) applied as 0.5 if the snow melting ratio was >0 (melting snow condition) unless it applied as 0.7 in the case of dry snow in the updated CLASS in CLASSIC (Eq. (m)). The CLASSIC classified three categories: total albedo, visible albedo, and near-infrared albedo that have typical values for fresh snow, old dry snow, and melting snow as shown in Table 6 (Aguado, 1985;Robinson and Kukla, 1984;Dirmhirn and Eaton, 1975 Noah-MP implemented two options for parameterizing snow albedo from BATS and CLASS schemes. The BATS scheme of snow albedo accounts for fresh snow albedo, variation in snow age, solar zenith angle, snow grain size growth, and dirt or soot on snow (Niu et al., 2011). The CLASS scheme of snow albedo considers fresh snow albedo and snow age (Niu et al., 2011). Noah-MP uses an old version of the CLASS scheme as Eqs. (k) Table 5). Frozen lake albedos have the same fixed value as glacier albedo (Bonan, 1996), direct and diffusive albedos of an unfrozen lake are parameterized following the cosine of the solar zenith angle (Eq.
Soil albedo depends on the soil color and varies with the moisture content of the first soil layer (Bonan, 1996). Soil albedo is calculated differently according to the change in soil moisture content (∆ = 0.11 − 0.40θ 1 > 0) based on the saturated albedo and dry albedo determined according to the soil color (Eq. (w)). The soil color of the CLM 5 is pre-specified for the CLM grid using the method of Lawrence and Chase (2007). The soil color is adjusted to one of 20 averaging the fitted monthly soil colors over all snow-free months and compared with the MODIS monthly local solar noon all-sky surface albedo described in Strahler et al. (1999) and Schaaf et al. (2002).
CLM 5 simulates the snow albedo and solar absorption within the snow layer with the Snow, Ice, and Aerosol Radiative Model (SNICAR), which was implemented in CLM version 3 (Flanner, 2005;Flanner et al., 2007). SNICAR includes the twostream radiative transfer solution of Toon et al. (1989) and uses the theory of Wiscombe and Warren (1980). The albedo and the 355 vertical absorption profile depend on the solar zenith angle, the albedo of the surface under the snow, the mass concentration of atmospheric-deposited aerosols (black carbon, mineral dust, and organic carbon), and the ice effective grain size (Lawrence et al., 2018).
JULES has two types of scheme: the diagnostic albedo scheme and the prognostic albedo scheme. The diagnostic albedo calculates snow albedo as Eq. (x1) in Table 5 using surface temperature (T g ), snow aging parameter (k 1 ), snow albedo threshold calculated using the cosine of the zenith angle (µ) and grain size (r). The current state of grain size is calculated every time step with snowfall rate S f by

370
where d 0 is fresh snow depth and r 0 is 50µm for fresh snow (2.5 kg m −2 fresh snow mass required to refresh the albedo) and G r is the empirical growth rate of snow grain area expressed as In Equation (16), A s is set to 0.23 * 10 6 µ m 2 S −1 . The prognostic albedo scheme in JULES consider direct and diffuse albedo from visible and near-infrared area same as BATS and CLM5. This scheme incorporates ice grain size affected by surface 375 temperature and snowfall rate as the significant factor to calculate snow albedo.
ISBA calculates snow albedo from 3 distinct bands (Eqs. (z1), (z3), and (z5) in Table 5), the first band with the ultraviolet and visible range (0.3-0.8 µm), the second and the third band with near-infrared ranges (0.8-1.5 and 1.5-2.8 µm) following Brun et al. (1992). For the first snow layer, snow albedo account for snow optical diameter, the near surface atmospheric pressure, and the age of the first snow layer in days (Eq. (z1)). Only first snow layer considers atmospheric pressure and time of snow 380 albedo decay ratio by impurities. The albedo scheme for the second and the third snow layer take into account snow grain size.

Snow density parameterization in LSMs
Snowfall amounts are measured with the snow depth and SWE when the snow melts. Typical snow rates are 1 cm h −1 or 0.8 390 mm h −1 . The LSMs normally calculate snow depth by dividing the SWE by the snow density. The density of newly fallen snow is approximately 20-300 kg m −3 , and the dry snow density is about 60-120 kg m −3 with the condition of low-to-moderate winds, while it increases with the condition of wet snow, sleet, and wind-packed snow. When precipitation mixes snow and rain, the new snow density exceeds 300 kg m −3 . When snow falls by strong winds, the snow particles become fragmented and small-grained. As the wind speed of 1 m s −1 increases, the density gradually increases to 20 kg m −3 (Jordan et al., 1999).

395
Therefore, snow characteristics (dry/wet) and wind conditions in addition to temperature affect the new snow density.
The LSMs considered different variables for calculating the snow density and snow depth (Table 7). Noah LSM, HTESSEL, CLM5, ISBA, and Noah-MP used air temperature to calculate the new snow density. Among these LSMs, HTESSEL, CLM 5, and ISBA had in common that they considered both the temperature and wind effects when calculating the new snow density. Jordan et al. (2008) also recognized the effect of wind speed, and the variation of wind speed caused by air pressure affects 400 the snow depth on the slope, leading to snow erosion on the slope's upwind side and snow deposition on its lee side. BATS includes complex snow physics such as SWE, snow grain size, and the effects of dust and soot for calculating snow density changes. CLASS parameterizes the snow density updated over time with the maximum snow density calculated by the snow depth and the empirical constant according to the snow temperature. The time tendency of snow density in ISBA is affected by the snow density, vertical stress, snow viscosity, and wind-driven compaction rate in each layer. The snow density over a 405 timestep in JULES is updated by compactive viscosity, ice and liquid water contents for the mass of snow above the layer, snow density, and the temperature of the k-th snow layer.
The LSMs parameterize the new snow density differently to estimate the snow depth as in Table 8. The new snow depth in Noah LSM is calculated by dividing the new SWE by the new snow density (Eq. (a)). If the sum of the existing and the new snow depth is < 1.0E −3 , the LSM took the higher value of the snow density between the new and the existing snow density.  Table 8). According to the parameterization method, Figure 3 illustrates the relationship between temperature and new 415 snow density. The minimum value of density of new snowfall is 0.05 g cm 3 , and the density of new snowfall tends to increase as the temperature increases.
In HTESSEL, the snow depth value is inversely proportional to the snow density and snow cover fraction and is proportional to SWE (Dutra et al., 2009) (Eq. (c)). When calculating the snow depth, using the snow density and snow mass (similar to the SWE) is similar to those of Noah LSM, and the consideration of the snow cover fraction is similar to that in CLM 5. Eq. (d) 420 shows the equation for obtaining the new snow density in HTESSEL. This equation is converted from a constant value to an equation using the equation of CROCUS (Brun et al., 1989(Brun et al., , 1992 and is a parameterization method that expresses the new snow density as a function of the temperature and wind speed (Dutra et al., 2009). In Eq. (d), a sn is 109 kg m −3 , b sn is 6 kg m −3 , and c sn is 26 kg m −7/2 s 1/2 . At this time, the snow density is limited to the range 50-450 kg m −3 .
The BATS model calculates the snow depth as shown in Eq. (e) in Table 8  CLASS also derives snow depth dividing SWE by snow density (Eq. (g) in Table 8). In the case of the maximum snow 430 density, CLASS 2.7 used 300 kg m −3 as a fixed value, whereas CLASS 3.6 parameterized as a function of the snow depth (Brown et al., 2006;Verseghy, 2012). If snow temperature below 0 • C, snow density decreases exponentially with time from the maximum mean density which is the calculated background snow density (Eq. (h)). The equation is derived from field measurements of Longley (1960) and Gold (1958) and express similarly to the albedo. Bartlett et al. (2006) found that the (i)). This empirical constant is consistent with the results of Wakahama (1968) and Colbeck (1973), which both showed a tendency to increase the compression rate of snow density with wet snow (Brown et al., 2006).
Noah-MP calculates snow density or snow depth considering destructive metamorphism by temperature change, compaction 440 by the weight of the overlying snow layers, and snow melting metamorphism according to Anderson et al. (1976). Noah-MP's snow layer determines by the total snow depth (Yang and Niu, 2003). If the snow depth is <0.045 m, snow layer does not exists.
If the snow depth is ≥ 0.045 m, the first snow layer has a layer thickness as the total snow depth. When snow depth is ≥ 0.05 m, two layers produces, and when snow depth is ≥ 0.15 m, a third layer is created. The third snow layer is thick due to the compaction following Sun et al. (1999) (Niu et al., 2011). Snow depth is dividing snow on the ground by bulk density snowfall which is set to 120 when fresh snow density is smaller than 120 and parameterized as the following fresh snow density (Eqs. (j) and (k)). Fresh snow density in Noah-MP is calculated according to the Hedstrom and Pomeroy (1998) by Eq. (l) in Table   8. Fresh snow density mainly account for the surface air temperature and freezing temperature.
CLM 5 calculates the snow depth next time, adding new snow depth to the snow depth at the current time (Lawrence et al., 2018). CLM 5 parameterizes the new snow depth as the rate of change (∆z sno , m), which is a function of the bulk density of in Table 8. Snow density scheme in JULES account for the reference snow temperature (k s = 4000 K), the snow melting point 460 (T m = 273.15 K), the temperature of the k-th snow layer (T k , K), reference snow density (ρ 0 = 50 kg m −3 ), snow density of the k-th snow layer (ρ k , kg m −3 ), and compactive viscosity (η 0 = 10 7 Pa s), and the mass of snow above the middle layer where I k (kg m −2 ) and W k (kg m −2 ) are an ice content and a liquid water content, respectively (Best et al., 2011). The layer thickness (d k , m) which have the relationship with density and mass expressed by ρ k d k = I k + W k . This formulation is identical meaning with Noah LSM, BATS, CLASS, and Noah-MP.

465
The newly fallen snow density in ISBA is identical with fresh snow density scheme in HTESSEL (Eq. (d) in Table 8).
ISBA expresses the time tendency of snow density as Eq. (r) like JULES, but the specific terms are different. The snow density changes by snow compaction due to vertical stress, the snow viscosity change (Brun et al., 1989), and wind-induced densification of snow layers (Brun et al., 1997). Especially, compaction rate of wind-driven snow densification parameterized as τ W (i) which describes precisely in Appendix B (Wind-induced densification of near-surface snow layers) from Decharme and vertical stress is also considered in JULES. It accounts for layer snow depth changes per time step (∆z), layer snow density (ρ sn ) and the acceleration of gravity (g) in every snow layer, but only the first snow layer used half mass of the uppermost layer as vertical stress effect. The snow viscosity in Eq. (u) is calculated by reference snow viscosity (η 0 = 7622370 Pa s), reference snow density (ρ 0 = 250 kg m−3), snow density in i-th layer (ρ sn (i)), snow temperature in i-th layer (T sn (i)), and the constants 475 (a η = 0.1 K −1 , b η = 0.023 m 3 kg −1 , and ∆T η = 5 K). The decrease of viscosity with the presence of liquid water (f W (i))is calculated by Eq. (v), where W l max (i) (kg m −2 ) and W l (i) (kg m −2 ) refer to the maximum liquid-water-holding capacity and liquid water content in i-th layer.
ables of snow models in the supplementary table from Menard et al. (2021). Figure 4 shows that only a few LSMs used specific variables to discern the snow cover fraction, so snow model developers can easily compare differences among LSMs. The eight land surface models examined in this study calculate the snow cover fraction through "the threshold of the SWE and the SWE," for Noah LSM, and HTESSEL, "the critical value of the snow depth and snow depth," for CLASS. In addition, BATS, ISBA, JULES and Noah MP considers "the snow depth and roughness length" for deriving snow cover fraction (Figure 4). Noah MP 485 additionally concerned snow density as the melting factor, and BATS and ISBA used the vegetation fraction separating the bare soil ground and vegetation covered ground for applying different roughness length according to the vegetation type and snow depth. Therefore, it is crucial to properly set the parameter values for "the threshold of snow depth or SWE" and the roughness length depending on the land surface type.
It is significant to derive an optimal parameter value that can reflect the local characteristics based on the same land type 490 from the observation data. For instance, accurate snow cover estimation is crucial because snow cover fraction affects snow sublimation, soil heat flux, surface emissivity, and albedo in Noah LSM. That is why parameter improvements were attempted in Livneh et al. (2010) and Wang and Zeng (2010). The parameter values according to the land cover type provided in Noah LSM were set differently according to the height of the vegetation. The current version of the Noah LSM reflects the SWE for snow cover parameterization but does not consider the snow depth, snow density, roughness length, and green vegetation frac-495 tion. Therefore, considering snow depth (CLASS, RUC, ORCHIDEE-I, SPONSOR in Figure 4) or snow density (HTESSEL, EC-EARTH, Noah-MP, CoLM in Figure 4) for snow cover fraction parameterization could be better because these variables have timely varying characteristics due to the temperature. In addition, current state-of-the-art LSMs considers fixed roughness length value according to the vegetation cover (JULES, BATS, ISBA, and Noah MP). Therefore, a seasonal variation of roughness length needs to calculate snow cover fraction for reflecting vegetation effect.  The land surface model's newly added snow cover parameterization method reflects the sub-grid topography variation in CLM 5. Ma et al. (2019) used this method which considers the sub-grid unit's terrain characteristics to calculate the snow 510 cover fraction. However, this method has a limitation due to only taking a variable called "elevation" into account as the standard deviation for the feature of the geomorphic setting. Roesch and Roeckner (2006) also measured snow cover fraction in the ECHAM climate model using the standard deviation of the sub-grid elevation and SWE. Douville et al. (1995) proposed a formula with the irregular snow cover distribution for the mountainous areas, including snow depth, roughness length, and the standard deviation of the subgrid orography (m). 515 Clark et al. (2011) presented a list of predictors related to the cause of spatial variation in snow depth according to spatial scale, as shown in Table 9. The list included both the elevation and slope among the significant predictors for snow temporal and spatial distribution. Since the spatial variation of snow depth is highly relevant with snow cover, we can apply essential predictors to the snow cover fraction parameterization. Anderson et al. (2014) also provided insights into the snow physics process that regulates the relationship between snow distribution and topographic characteristics and argued that snow cover 520 correlates with altitude, aspect, canopy density, and snow accumulation/ablation. Besides, Mott et al. (2018) suggested that the characteristics of snow cover change according to the mountain range scale, the ridge, and the slope scale may differ.
This idea provides worth variables considering when applying scales to interesting areas in LSMs. First, snow accumulation depends on the climate, elevation, and vegetation at the mountain-range scale (López-Moreno et al., 2008;Clark et al., 2011;Anderson et al., 2014). Especially, elevation is the most significant variable for precipitation pattern at the mountain range scale 525 (Mott et al., 2018). Second, snow deposition patterns at the ridge-scale reveal more considerable spatial variability of snow or enhanced snow deposition over leeward slopes by preferential deposition of snowfall (Mott and Lehning, 2010;Gerber et al., 2017). Third, at the slope scale, the highest model resolutions show the snow redistribution processes such as saltation and turbulent suspension (Mott et al., 2018).
Noah LSM, CLM 5, and JULES calculates the land surface albedo by summing the ground albedo without snow and the 530 snow albedo using the snow cover fraction. Other LSM models using the Noah LSM's snow albedo decay scheme were VIC and DHSVM. Livneh et al. (2010) introduced the parameter value of albedo decay parameterization. Storck (2000), Sun et al. (2019), and Malik et al. (2014) also enhanced this parameter to fit the characteristics of the study area. Furthermore, maximum snow albedo was approximately 0.6-0.95 in other LSMs and journal articles (Sproles et al., 2013;Verseghy, 1991;Wigmosta et al., 1994;Andreadis et al., 2009;Yang et al., 1997;Bartelt and Lehning, 2002). In Noah LSM, snow albedo is affected by 535 snow cover; it affects the latent heat of the freezing rain, the heat flux of the snow surface, and the phase change heat flux for snow melting. Therefore, the maximum snow albedo is the main parameter of snow albedo parameterization, so setting the maximum snow albedo's optimal value is essential.
In addition to the HTESSEL, BATS, and CLASS models, ISBA includes a prognostic albedo scheme in which the previous snow albedo value affects the current snow albedo value. At first, CLASS did not consider the snow albedo changes according 540 to the spectral range, similar to HTESSEL (Eqs. (k) and (l) in Table 5). Noah-MP also used this scheme (Wang et al., 2020).
After that, CLASS applied the difference between visible and near-infrared snow albedo according to spectral range like CLM 5, BATS, JULES, and ISBA; it used an empirical formulation that has the difference between dry snow (Eqs. (n) and (o) in Table 5) and melting snow (Eqs. (p) and (q) in Table 5). Melton et al. (2019) used the spectral method instead of the empirical exponential decay function in order to apply the results of Wiscombe and Warren (1980) that the growth of snow particles, 545 soot/dirt precipitation affects snow albedo. BATS's snow albedo parameterization method was actively accepted in CLASS, and there was a slight difference in applying to freshly fallen snow with a near-infrared value of 0.73, unlike BATS's 0.65. In addition, the structure of the equation reflecting the effect of the solar zenith angle was slightly different. The CLASS model in CLASSIC shows the most significant difference from other land surface models in that they reflected canopy albedo (Bartlett and Verseghy, 2015;Melton et al., 2020). In particular, CLASS has the lowest albedo bias among the JULES-I, JULES-GL7, 550 JULES-UKESM, Surfex-ISBA, CLM 5, CLASS, and HTESSEL from Menard et al. (2021) (Figure 5). Thus, in future studies, the improvement of snow albedo simulation that reflects the effect of vegetation is crucial.
Since many LSMs calculate the snow depth as the SWE ratio to the snow density, it is vital to parameterize the new snow density affected by fresh snowfall. The new snow density's minimum value differed slightly depending on the type of LSM or the study area, so it was necessary to adjust the parameter value. In the snow density parameterization, although CLM 555 5, HTESSEL, and ISBA took wind and temperature into consideration when calculating the new snow density (Noah LSM, Noah-MP only considered temperature), various other variables may affect the new snow density. Föhn (1976Föhn ( , 1985 raised the question of the representativeness of the measurement of the new snow depth in precipitation gauges in the mountains because they ignored the topographical characteristics such as the mountain peak, the wind-slope, and the leeward slope or depressions. Since records from rainfall gauges surrounded by valley plains help infer the depth 560 of new snow depths over large areas (Föhn, 1976(Föhn, , 1985. Therefore, it is imperative to reflect the influence of topographic conditions' variation (elevation and slope) for mountainous areas' snow depth. Clark et al. (2011) proposed to derive snow depth by considering topographic variations such as wind slope vs. leeward slope, depression, crests, and the variation of critical elevation. Spatial variation of snow depth and temporal consistency of snow cover diversified due to snow freezing temperature, slope, aspect, and radiation. Hojatimalekshah et al. (2020) examined the relationship between tree canopy and 565 snow depth and Kantzas et al. (2014) investigated the change of the snow regime according to vegetation dynamics through a ground model; these studies show that consideration of the vegetation effect is necessary for the snow parameterization method. Therefore, snow density parameterization regarding snow depth needs to include vegetation density, Leaf Area Index (LAI), Stem Area Index (SAI), vegetation type, and snow behavior in sheltered areas. In breif, geomorphic and vegetation-related variables can affect the spatial-temporal variation of the snow depth.

Conclusions
Examining the eight LSM's snow physics parameterizations allowed us to find each parameterization's vulnerabilities. Depending on the purpose of each LSM, we presented a comparison table so that readers can identify the variables already reflected in various LSMs which have enough information on snow-related schemes. We mainly focused on parameterization methods related to snow cover, snow albedo, and snow density. This study confirmed that snow cover, snow albedo, and snow 575 density parameterization differ for each land surface model in terms of complexity and significant variables or parameters. The comparative analysis for each element is summarized as follows.
In the snow cover case, the snow decay scheme's parameterization method in CLM 5 considers the subgrid's topographical variation. Besides JULES, BATS, ISBA, and Noah MP used snow depth and roughness length for calculating the snow cover fraction. BATS and ISBA consider vegetation fraction, while Noah MP accounts for the melting factor using snow density. In terms of snow density, while Noah LSM and Noah-MP consider only temperature, HTESSEL, CLM 5, and ISBA reflect temperature and wind to parameterize the density of new snowfall. ISBA and JULES include the internal snow processes such as vertical stress, snow temperature, and snow viscosity. In particular, ISBA accounts for the wind-driven snow compaction precisely. BATS calculates snow density using the growth of snow particle size and the accumulation of dust and soot. CLASS 590 calculates the maximum snow density as a function of snow depth, which increases exponentially over time with wet snow.
Noah LSM, BATS, CLASS, Noah-MP, HTESSEL, and CLM 5 considers SWE or snowfall rate and snow density for parameterizing snow depth, and HTESSEL and CLM 5 use the snow cover fraction additionally to the snow depth parameterization.
b The horizontal scale associated with variability in freezing levels is often much larger than the vertical scale.