Modification and Validation of the Soil–Snow Module in the INM RAS Climate Model

: This paper describes the modification of a simple land snow cover module of the INM RAS climate model. The possible liquid water and refreezing of meltwater in the snow layer are taken into account by the proposed parameterization. This is particularly important for modelling the transition season, as this phenomenon is mainly observed during the formation and melting of the snow cover when the surface temperature fluctuates around 0 ◦ C. The snow density evolution simulation is also added. This parameterization is implemented in the INM-CM snow module and verified on observation data using the ESM-SnowMIP-like protocol. As a result, the INM-CM mean climate snow melt periods are refined, particularly in middle and high latitudes. The snow-covered area according to the model is also improved. In the future, a modified version of the land snow module can be used, coupled with a snow albedo model that takes into account snow metamorphism. This module can also be applied to sea ice snow.


Introduction
Simulations using climate models are primary tools for studying climate change.This approach supplements natural observations with numerical experiments, allowing for significant advancements in the study of Earth's climate.Snow cover has a powerful impact on the climate system, affecting many processes on the Earth.
Seasonal snow cover greatly influences energy balance at the surface (in particular, albedo) [1], the hydrological cycle [2], and other phenomena, for example, winter-time circulation anomalies over mid-high latitudes [3,4].
Snow cover is a complex system controlled by various factors.For example, during transitional seasons with temperature fluctuations around 0 • C, when the snow melts, not all of the meltwater drains into the soil, but some of its quantity is retained in the snow layer and can be refrozen.These phenomena change the morphology of the snow cover, which can eventually lead to later snowmelt and, thus, affect the hydrological cycle and surface radiation balance.On the other hand, over time, snow becomes denser, absorbs moisture, and can also be contaminated, leading to changes in its optical properties and surface radiation imbalance.
Snow cover models are actively developed both as stand-alone products and as parts of larger projects, such as land surface models (LSMs) or global climate models (GCMs).Verification of snow models takes place in several stages.For example, there are independent model intercomparison projects for land surface (LS3MIP, land surface, snow, and soil moisture model intercomparison project) [5] and snow (ESM-SnowMIP, the Earth system model-snow model intercomparison project) [6] linked to the global climate model comparison project (CMIP6, coupled model intercomparison project-Phase 6) [7,8].
ESM-SnowMIP [6] is an international project that aims to simulate and evaluate modern snow modules integrated into the Earth system models.Within the framework of this project, model data is compared with both local (station) and global observations.One of the objectives is to identify the physical processes that have the greatest impact on the simulation of snow cover features.Another goal is to better assess the Earth system feedback associated with snow in multi-model experiments.
There is a family of global climate models developed at the Marchuk Institute of Numerical Mathematics RAS (INM-CM).This group of models is constantly evolving, and today it includes a set of coupled models with different spatial resolutions, as well as various special modules.Model version INM-CM-48 [9] participates in CMIP6 [7] and shows good results compared to other models [10].This version includes a soil module with a previous generation of snow cover parameterization.Its description and verification results are given in [11,12].Future climate change modelling with the INM RAS climate model is discussed in [13], and changes in the Russian Federation are presented in [14].There is a new version of the INM RAS climate model, INM-CM-60 [15].This version differs from the previous one, INM-CM-48, with changes in the cloud and condensation scheme, aerosol evolution calculations and their indirect effects, newly implemented land snow and atmospheric boundary-layer parameterizations, and other schemes.
When a new version of any climate model is prepared, it is often not easy to state the reason for climate change from a previous one, because a lot of updates are introduced to it.Here, the response to only one improvement of land snow cover parameterization for the INM RAS climate model is estimated.
This paper proposes a very simple parameterization of snow water refreezing and snow compaction, capable of improving results in single-column simulations proposed in ESM-SnowMIP [6] and as part of the INM RAS Earth system model.The described parameterization is applied to land surface snow only.The new INM-CM snow module is verified against meteorological observations in stand-alone tests and reanalyses, as well as satellite data implemented as part of the INM RAS climate model.It should be noted that the original INM-CM 1D soil module was tested earlier, but the focus was not on snow properties but on surface flows, soil temperature, and freezing depth, which are only indirectly dependent on snow properties [11].However, here, a direct investigation of snow properties in a single-column mode is presented for both snow parameterization versions.The main novelty of this study is that the implementation of the modified snow parameterization leads to an improvement in the snow simulation in the global climate model.

INM-CM Modification
The INM RAS climate model consists of three main modules: atmosphere, ocean, and aerosol dynamics.In this study, the INM-CM-48 [9] version is taken as the base.The spatial resolution of its global atmosphere circulation model is 2 • × 1.5 • in longitude and latitude and it has 21 vertical σ levels.The interactive aerosol module [16], which simulates the evolution of concentrations of 10 substances, is included in the atmosphere model.The ocean global circulation model has a horizontal resolution of 1 • × 0.5 • in longitude and latitude and 40 vertical σ levels.It includes dynamic and thermodynamic modules for the sea ice with elastic-viscous-plastic rheology [17].

Initial Version of INM-CM Soil-Snow Module
The soil-snow module of the INM-CM family is a part of the atmosphere dynamics model.This one-dimensional module [11,12], describes processes of heat and moisture transfer in soil, vegetation, and snow cover, as well as heat and moisture exchange with the atmosphere.The snow module described here is used as the top layer of the soil-snow module.
The main characteristics of the snow cover are the SWE (snow water equivalent) layer thickness (S, [m]), its density (ρ sn , [kg/m 3 ]), and its depth (H sn , [m]).In the original version of this module, the SWE calculation is based on the balance equation (Equation ( 1)), which takes into account snowfall (P-atmospheric precipitation at a negative temperature of underlying surface), snowmelt (M-snowmelt intensity), as well as snow sublimation (E -latent heat flux at the surface spent on snow sublimation, L-specific heat of snow sublimation, ρ w -water density).
Snow starts to melt when the underlying surface temperature becomes equal to 0 • C, and the heat balance at the surface is positive.In this scenario, all melted snow, along with any rain that falls on the snow cover, enters the soil immediately.The snowmelt intensity, M, is determined by excess heat flux at the surface, ∆E > 0, and the specific heat of snow melting, λ (Equation ( 2)).
In this version, the snow density is assumed to be a constant, i.e., ρ sn = 185.4kg/m 3 .Note that the snow albedo (α sn ) in the INM-CM is determined by a simple parameterization, only depending on the surface temperature, as follows:

Description of Modified Snow Module in INM-CM
The modification of the INM-CM soil-snow module described here can be divided into several parts: snow melt accumulation processes and changes in snow cover density over time.Firstly, the water equivalent of the snow cover is calculated based on the surface heat balance equation.Then, using the updated amount and composition data of snow cover, the density is updated.Figure 1 shows the conceptual flowchart of the described parameterization.

Snow accumulation-melting processes.
The new version of the soil-snow module takes into account the porous structure of the snow layer.It is assumed that water gradually seeps through the snow layer, rather than instantly draining into the soil during snow melt.As the meltwater percolates, it can refreeze and release heat to the surrounding snow.Rainwater is presumed to infiltrate the soil without refreezing.
It is assumed that the amount of liquid water in the snow layer (S wat , [m]) cannot exceed a critical value, and any excess water immediately drains into the soil.The critical value of liquid water in snow (S max wat ) depends on its amount (e.g., SWE) and its porosity (ε sn ).These quantities are estimated using the following empirical hydrological formulae [18,19]: The consideration of a possible liquid fraction in the snow layer requires the correction of the balance equation for SWE.Firstly, at each time step in the soil-snow module, the surface heat balance is calculated without taking into account possible phase transitions.Then, in the case of phase transitions, possible energy excess in the heat balance, ∆E > 0, is spent on snow melting (as in the original version, the intensity of snow melting is determined by Equation ( 2)), whereas the energy deficit, ∆E < 0, is spent on the refreezing of meltwater contained in snow pores (S r f rz , [m]).The refreezing intensity in this case is determined by the following equation: Taking into account the described phenomena, the SWE balance equation is written as follows: Here, S atm corresponds to snow in its usual form, falling from the atmosphere as snowflakes, M = max 0, M − S max wat ∆t is the flow of excess meltwater, which is not stored in snow pores and drains into the soil, Note that in the modified version, the condition for the onset of melting remains unchanged: the underlying surface temperature is 0 • C, and the surface heat balance is positive.
Snow density evolution.In the modified version of the INM-CM soil-snow module, the snow cover is represented as a certain substance consisting of four fractions.It can simultaneously include snow in its usual form (both old and freshly fallen) as well as meltwater and refrozen snow.It is important to differentiate fresh and old snow: over time, it becomes denser due to snow grain metamorphism into a more rounded shape, compaction under the weight of over-layers, and refreezing of meltwater [20].Thus, the density of freshly fallen snow is 50-100 kg/m 3 , and it can reach 500-700 kg/m 3 for old snow (at the age of 1-2 months).
Snow compaction over time is described with the parameterization, following [21], where there is slow densification under loads resisted by the compactive viscosity: Here, ρ sn is the snow density, calculated in kg/m 3 , M sn is the snow mass per unit area (in kg/m 2 ), g = 10 m/s 2 is the acceleration of gravity.The parameter η 0 is the viscosity coefficient at ρ sn = 0 and T sn = T C , T sn -snow layer temperature (in K), T C = 273.2K.The constants C 1 and C 2 equal 0.08 K −1 and 21 × 10 −3 m 3 /kg, respectively.The Equations ( 8) and ( 9) can be rewritten in discrete form as follows: Here, S is the SWE in m, t • sn is the snow layer temperature in Celsius, and the value C 0 is 0.12 × 10 −3 1/(m•s) (this corresponds to η 0 = 8.3 × 10 7 Pa•s) similar to the SWAP (soil-water-atmosphere-plant) land surface model [18].
Snow that has not fallen during the current model time step is considered to be old.The average snow cover density over the layer, taking into account all four described fractions, is calculated as an SWE-weighted average: Here, ρ old is the old snow density calculated using the evolutionary model ( 10), ρ new = 100 kg/m 3 is freshly fallen snow density, ρ w = 1000 kg/m 3 is the water density, ρ ice = 917 kg/m 3 denotes ice density, δ old , δ new , δ wat , δ r f rz is the ratio (in total water equivalent) of old and fresh snow, meltwater, and refrozen, respectively.The value for fresh snow density is chosen according to the limit of the critical value of liquid water in snow, as follows:

Verification of INM-CM Soil-Snow Module
To verify the INM-CM soil-snow module, numerical experiments similar to the ESM-SnowMIP project [6] are performed.For this purpose, both local simulations for 10 reference sites (Table 1) and global simulations are used.The calculations are conducted with the original version of the INM-CM soil-snow module as well as with the modified version described in this paper.Note that INM-CM is not a contributor to ESM-SnowMIP.

Local Simulations
For local simulations, observation data from 10 meteorological stations (Table 1) are used.They are located in different regions of the world: alpine, arctic, boreal, and marine regions.For these sites, all necessary data are available [29], for use as input to the model and for evaluating the quality of the simulation.Atmospheric forcings (temperature, pressure, humidity, precipitation, radiation fluxes, wind speed) are used as input data.Snow cover and soil parameters are required as outputs for comparison with observations.In these experiments, the model is run in a single-column mode using only the INM-CM soil-snow module.

Global Simulations
For the global experiment, simulations were performed using the global climate model INM-CM-48 with original and modified soil-snow modules.An ensemble calculation with 6 runs was carried out for 20 years, starting from a prepared control point on 1 January 1997.The start date and duration of the simulations were chosen as such because observational data on global snow cover have been more reliable since the mid-1990s, as well as due to the available periods of CAMS [30] (2003-2017) and NSIDC-Blended5 [31] (1981-2010) reanalysis data for the comparisons.
For global verification, the spatial distribution of SWE, as well as global metrics, such as the snow-covered area, are used.The simulation results are compared with the results of the ESM-SnowMIP experiments [6], reanalysis data, and observations.
In order to compare the snow-covered area, it is necessary to estimate the snow fraction ratio (SFR) of the model grid cell.For this purpose, the following empirical parameterization [32] is used in the output post-processing stage: Here, H sn denotes snow depth (in meters), ρ sn and ρ new [kg/m 3 ] denote current and fresh snow densities, and C = 2.5 z 0 , z 0 = 10 −2 m-ground roughness length.Here, the values 2.5 and 1.6 are empirical constants according to [33]).The optimal limit value for the SWE of 5 mm is used, corresponding to minimum grid cell coverage [6,34].If the SWE is less than this value, the snow cover fraction is assumed to be 0.

Local Simulations
The results of calculations using the soil-snow module of the INM RAS climate model with input data from the ESM-SnowMIP sites show good agreement with the observed data for snow cover and surface and soil temperatures at these stations.Both versions of the snow model (the original INM-CM-48 and the modified one) correctly reproduce the thickness (Figures 2 and 3) and mass (Figures 2 and 4) of the snow cover and the surface temperature (Figure 5).The modified version shows better agreement with observations for the majority of sites, particularly during the spring months when refreezing and changes in snow structure usually occur.According to this version, the snow cover melts an average of 1 month later, which is in agreement with the data from the sites.However, there is an overestimation of the snow amount at some sites.The reasons for this may be the use of a single-layer snow scheme and an overly simple snow albedo parameterization, which leads to slower melting.Nevertheless, as a result of modifications, peak values of SWE and H sn are corrected (most noticeably for the stations Sapporo, Swamp Angel, and Weissfluhjoch), which indicates a more accurate simulation of snow accumulation by the new version of the model.Note that the clarification of modeled snow cover mass is also associated with more accurate modelling of its density, but not only with changes in its depth.
To compare the accuracy of snow cover modelling by the INM-CM and similar models participating in the ESM-SnowMIP project [6,35], basic quality metrics for soil time series data (NRMSE-root mean square error normalized by the standard deviation of observation data, bias, and correlation coefficient) [36] are calculated for the daily height (H sn ), snow water equivalent (SWE), and surface temperature (T sur f ). Figure 6 shows values for the original and modified versions of the model with respect to observations for available stations, as well as averaged over all such stations.These mean magnitudes are also shown in Table 2.Note that only the days with snow depth greater than 0.1 m and surface temperature below 0 • C are included.According to similar plots for snow model errors in [35], the modified version of the INM-CM has an accuracy comparable to other known models.During testing of the INM-CM soil-snow module, it was noticed that for the mountainous area (stations Col de Porte, Reynolds Mountain East, Senator Beck, Swamp Angel, and Weissfluhjoch), the maximum snow cover strongly differs from the observations, although there is a good agreement on average over the year.To correct this effect, additional model tuning was carried out.Initially, in the INM RAS model, the sensible heat flux between the Earth's surface and the air F heat was calculated using the following equation [11]: Here, ρ a [kg/m 3 ] is the air density at 2 m, C T is a dimensionless coefficient of heat transfer between the surface and atmosphere, V a [m/s] is the wind speed module at 2 m, and θ a and θ s are the potential temperatures of the air at 2 m and of the surface, respectively.Under typical snowmelt conditions (where the surface temperature is below 0 • C, the snow temperature is at 0 • C, and there is stable stratification), C T is approximately 10 −3 or less.
Although the classical Monin-Obukhov theory assumes stable stratification, there is strong turbulence in the mountains, leading to an increased exchange between the surface and atmosphere.Thus, the overestimation of maximum snow cover, as well as the delayed start of melting for stations in the mountains, may be associated with an underestimated sensible heat flux arising from too low of an exchange coefficient.Therefore, for heights over 500 m above sea level, the dimensionless coefficient of heat exchange is bounded from below: C T ≥ 5 × 10 −3 .As a result, the maximum snow amount is corrected (Figure 7), and the quality value averaged over the snow season is also refined (Figure 8).The INM-CM soil-snow module also outputs soil temperature at 23 horizons from 1 cm to 10 m in depth.For most of the stations considered, temperature observation data are available for at least one level within the soil.So, the soil temperature is compared with observations at a depth of 10 cm.There is a strong agreement (Figure 9) with all available site observations for both model versions.Figure 10 shows that the updated snow model does not lead to a sensible decline in the quality of the upper soil layer temperature simulation.However, the positive bias in the soil temperature at 10 cm is changed by a negative one.Snow refreezing in the new version leads to later snowmelt followed by late soil heating in the spring, resulting in an increase in soil water as well as surface latent heat flux during the spring and early summer.Moreover, the correlation coefficient for the temperature of the soil's top layer decreases from 0.96 to 0.86.This may be due to the soil settings specific to the original version of the snow model.

Global Simulations
Simulations using the global climate model show results similar to the single column one.The changes in the soil-snow module also lead to later snowmelt (up to a month later in some regions) and more intensive snow cover formation.These effects are especially notable in middle and high latitudes and in the highlands.Now, in the Arctic, as well as in the mountains along the west coast of Canada and Alaska and the Himalayas, snow cover persists into June, whereas in the old version, almost all snow melted in May (Figure 11).
The contribution of the described processes in snow cover formation is also noticeable in late autumn-early winter in the Arctic regions of Eurasia and the Himalayas.These results are consistent with observational data in various regions, e.g., [37].Figures 12 and 13 show good agreement in the forecasts of snow-covered areas with the CAMS global reanalysis (EAC4) by the European Centre for Medium-Range Weather Forecasts [30] and the NSIDC-Blended5 5-product-based SWE database [31].The different averaging periods in Figures 12 and 13 are due to the time limits of the available model and reanalysis data.

Discussion
The results of the INM-CM soil-snow module are comparable with other snow-covered models.According to the ESM-SnowMIP model ranking by errors [35], the INM RAS climate model version with the new simple snow physics implementation could be placed at about 20th out of 27 models on the quality of SWE simulation.For surface temperature, this version could be ranked between 10th and 12th place out of 27.
In [36], the same soil and surface block [11] was incorporated into the weather forecast model SLAV, and data assimilation results for soil temperature were compared with observations from several groups of European stations.It was found that, in general, the correlation coefficient for the model and observed soil temperature ranged from 0.80 to 0.95, while the RMSE was about 0.5-1.5 degrees.This study reports similar correlation coefficient values for the model and observed soil temperatures, with similar or lower RMSE values.
The simulation of SWE and soil temperature in the intermediate complexity climate model IAP RAS is presented in [38,39].In the IAP RAS model, snow parameterization is similar to the original snow parameterization considered in this study.Detailed parameterization of soil heat and water conductivity is used.In [38], it is shown that the total amount of SWE in their model is simulated reasonably, but the IAP RAS climate model tends to underestimate SWE in Arctic regions and tends to overestimate SWE in middle latitudes, probably due to too simple parameterization of advection of moisture.In this study, the spatial distribution of SWE is simulated better due to the explicit treatment of moisture advection.Also, the IAP RAS climate model tends to overestimate the upper soil temperature in Arctic and Northern Siberia regions, and underestimate the soil temperature in central European Russia, likely because of the simplified treatment of atmospheric advection.The results presented in this study show a better simulation of upper soil temperature because of the explicit treatment of atmosphere dynamics.
Taking into account the different snow layer fractions allows us to describe the snowflake metamorphism, not only for dry but also for wet snow, with the modified version of the model (for example, similar to [40,41]).Snow particle size modelling is one of the approaches used to simulate the snow albedo affecting the radiation balance on the surface.The INM-CM could also take into account the aerosol contamination of snow [42], which is important for the albedo.However, this feature is not currently implemented in the current version of the model.All of these issues could become the subject of future studies that continue the research under discussion.
This study is devoted to snow lying on land, but its results can also be applied to describe sea ice snow with the necessary adjustments.This also accounts for snow blowing and redistribution due to the sea ice movement, puddling, etc.

Conclusions
The modified land snow module of the INM RAS climate model is presented.It takes into account the influence of liquid water in the snow and the snow layer densification over time.Here, several well-known parameterizations are used, but their combination may be too new to be incorporated into one snow model.
The land surface module with old and updated snow models was tested in singlecolumn experiments on data from 10 widely scattered sites, covering almost all land areas where snow cover can occur.Snow simulation is improved overall in the modified model, as evidenced by all key metrics for soil modules (except for topsoil temperature) at local sites.The simulation quality of local snow cover data is consistent with the spread of snow models [6], although it is not the best.
The improved version of the snow module was incorporated into the climate model, and long-term global numerical experiments were also carried out.Performed numerical simulations show that the quality of snow cover reproduction is improved with the new version of the INM-CM.As a result, the average climate dates of snow cover melting and the reproduction of a snow-covered surface area, as well as its mass, are improved.Additionally, the updated description of snow layer physics allows for more accurate modelling of related phenomena, such as the radiation balance at the surface.
The main conclusion of this research is that the incorporation of the selected snow parameterization into the climate model leads to improvements in snow simulations in both local and global climate model runs.Furthermore, the presented snow layer parameterization shows results comparable to other models.

Figure 1 .
Figure 1.Conceptual flowchart of the modified snow model.

Figure 2 .
Figure 2. Monthly average depth (a) and mass (b) of snow cover during the entire simulation period for the Col de Porte site, according to the INM-CM snow module (blue-original version, red-modified) and observations (black).

Figure 3 .
Figure 3. Annual variation of snow depths averaged over the simulation period, according to the INM-CM snow module (blue-original version, red-modified) and site observations (black).

2 WFJ_1996_2016Figure 4 .Figure 5 .
Figure 4. Annual variation of snow mass averaged over the simulation period, according to the INM-CM snow module (blue-original version, red-modified) and site observations (black).

Figure 7 .
Figure 7. Monthly average snow depth during the simulation period for the Col de Porte site, according to the different INM-CM snow module versions (including the version with the correction of the heat transfer coefficient, C T ) and observation data.

Figure 8 .
Figure 8. Quality metrics (normalized root-mean-square errors, bias, and correlation coefficient) for snow depth and SWE for mountain sites (for the original, modified, and modified with the additional tuning of the INM-CM snow module versions).

Figure 9 .Figure 10 .
Figure 9. Monthly average soil temperature at a depth of 10 cm during the simulation period for the Reynolds Mountain East site, according to the INM-CM snow module (blue-original version, red-modified) and observations (black).

Figure 11 .Figure 12 .Figure 13 .
Figure 11.(a) Mean monthly SWE according to the INM RAS global climate model with the modified version of the soil-snow module, (b) difference with the original version (INM-CM48), and (c) with the NSIDC-Blended5 reanalysis data (in [mm]; model data averaged over the ensemble and all simulation period; month-May).

Table 2 .
Quality metrics averaged over stations.
Figure 6.Quality metrics (normalized root-mean-square errors, bias, and correlation coefficient) for snow depth, snow water equivalent thickness, and surface temperature.