Modeling Snow Surface Spectral Reﬂectance in a Land Surface Model Targeting Satellite Remote Sensing Observations

: Snow surface spectral reﬂectance is very important in the Earth’s climate system. Traditional land surface models with parameterized schemes can simulate broadband snow surface albedo but cannot accurately simulate snow surface spectral reﬂectance with continuous and ﬁne spectral wavebands, which constitute the major observations of current satellite sensors; consequently, there is an obvious gap between land surface model simulations and remote sensing observations. Here, we suggest a new integrated scheme that couples a radiative transfer model with a land surface model to simulate high spectral resolution snow surface reﬂectance information speciﬁcally targeting multisource satellite remote sensing observations. Our results indicate that the new integrated model can accurately simulate snow surface reﬂectance information over a large spatial scale and continuous time series. The integrated model extends the range of snow spectral reﬂectance simulation to the whole shortwave band and can predict snow spectral reﬂectance changes in the solar spectrum region based on meteorological element data. The kappa coe ﬃ cients (K) of both the narrowband snow albedo targeting Moderate Resolution Imaging Spectroradiometer (MODIS) data simulated by the new integrated model and the retrieved snow albedo based on MODIS reﬂectance data are 0.5, and both exhibit good spatial consistency. Our proposed narrowband snow albedo simulation scheme targeting satellite remote sensing observations is consistent with remote sensing satellite observations in time series and can predict narrowband snow albedo even during periods of missing remote sensing observations. This new integrated model is a signiﬁcant improvement over traditional land surface models for the direct spectral observations of satellite remote sensing. The proposed model could contribute to the e ﬀ ective combination of snow surface reﬂectance information from multisource remote sensing observations with land surface models. spatial variation characteristics of land surface snow spectral albedo. Based on the new integrated model, we extend the simulation of snow spectral albedo to the whole solar spectrum region and reﬁne the wavelength interval of the snow spectral albedo simulation to 1 µ m. The narrowband snow albedo targeting MOD09GA reﬂectance data simulated by the integrated model spatially agrees well with the narrowband snow albedo retrieved by utilizing the MOD09GA reﬂectance data and has a smaller deviation in the time series. The results of our research prove the feasibility of simulating the snow radiation observed by satellite remote sensing while utilizing an integrated snow hydrological model and snow radiative transfer model based on ground-driven data. Our results suggest that if real-time meteorological forcing data and snow optical data can be obtained, the new integrated model will be able to predict spatiotemporally continuous snow spectral albedo and further predict narrowband snow albedo by targeting di ﬀ erent optical remote sensing satellite observations. We developed a new method representing an important attempt to simulate snow reﬂectance information directly from land surface models and ground-driven data targeting satellite remote sensing observations, e ﬀ ectively eliminating the impacts of multiple intermediate transformations on the accuracy of directly simulated remote sensing observations. The results of this study will be beneﬁcial for the application of multisource remote sensing data to snow radiation monitoring in land surface integrated models. Moreover, the new land surface integrated model can compensate for the shortcomings of satellite remote sensing observations in time series.


Introduction
Snow surface albedo can significantly affect climate and hydrological cycles at different scales [1]. Climate change and snow albedo participate in a positive feedback cycle: when global warming causes a decrease in the areas of snow cover and sea ice, the surface albedo decreases, the solar radiation absorbed by the Earth-atmosphere system increases, and the temperature rises, thereby increasing the global warming effect and causing further melting of snow cover and sea ice. Previous studies have shown that this positive feedback between snow albedo and climate change has led to significant declines biogeochemistry, human dimensions, and ecosystem dynamics [26]. The SNICAR model, based on the snow radiative transfer mechanism, was first introduced into CLM4.0 [41], greatly improving the snow albedo simulation accuracy of the CLM. In the above studies, with the aim of simulating the overall energy budget of snow radiation [25,26,42,43], most of the researchers considered snow albedo in only the VIS and NIR bands under direct or scattering conditions when simulating the land surface snow energy balance process [44] but lacked the ability to simulate the fine spectral reflectance of snow. In addition, these studies did not combine modeling with satellite remote sensing observations because their purpose was mainly to understand snow physical processes.
Based on the above analysis, there is an obvious gap; that is, the snow reflectance information observed by satellite remote sensing is the discrete spectral reflectance, while the snow reflectance information simulated by land surface models is mostly broadband albedo. To calibrate land surface models and compare them with satellite observations, the albedo must be estimated through a series of transformations of the satellite-observed discrete spectral reflectance and then compared with the modeled albedo. In summary, the traditional integrated models used to obtain snow surface reflectance information do not consider the role played by satellite remote sensing observations, and utilizing satellite remote sensing observations to obtain snow surface reflectance information rarely considers meteorological and snow evolution factors. Most of the studies on snow radiative transfer isolate the remote sensing observation process and the land surface forward simulation process; unfortunately, there are few effective combinations of remote sensing observations and land surface forward simulations in the available research on snow radiative transfer. Thus, the key to calibrating remote sensing data and comparing the data with land surface modeling results is in being able to directly simulate snow reflectance by targeting satellite remote sensing observations based on land surface models. Therefore, considering the spectral band observation range of different satellite sensors, an effective means to couple remote sensing observations with forward simulations would be to establish a set of high spectral resolution snow albedo simulation schemes targeting satellite remote sensing observations based on land surface models. This new integrated scheme can not only realize the simulation of satellite remote sensing observations with ground-driven data but also calibrate land surface models with remote sensing observations.
The GBEHM is an ecohydrological model specifically developed for the environment of alpine mountainous regions. This model proposes a new method for simulating the snow energy balance process [29]. Compared with other snow hydrological models, the GBEHM can better simulate snow hydrological processes in cold regions and achieves better accuracy when simulating the input parameters (such as the snow cover area, snow depth and snow water equivalent) of the snow radiative transfer process [29]. In addition, when simulating the snow radiative transfer process, the SNICAR model not only considers the detailed dynamic evolution of the snow grain size but also simulates the multiple scattering of multilayer snow cover, thereby yielding a better snow albedo simulation accuracy [32]. Therefore, the integration of these two models can effectively avoid the limitations of either model alone. In summary, the GBEHM can simulate snow energy balance and snow mass balance processes and can obtain the basic input parameters of the snow radiative transfer process, but it lacks the ability to simulate the snow grain size evolution and snow radiative transfer processes in detail; conversely, the SNICAR model can simulate the dynamic evolution of the snow grain size and the transfer process of solar radiation energy in snow, but it needs to calibrate and drive the snow mass and energy balance processes and snow process parameters and cannot independently simulate the spatiotemporal evolution of snow radiative transfer. Hence, the integration of these two models can effectively solve the separation between snow mass-energy balance processes and snow radiative transfer processes in snow albedo simulations. For the abovementioned reasons, we developed a new integrated snow radiative transfer simulation model utilizing both a snow hydrological module and a snow radiative transfer module together with ground-driven data, remote sensing data and meteorological element data. In addition, based on the new integrated model, we developed a set of snow spectral albedo simulation schemes with high spectral resolution in the solar spectrum region Remote Sens. 2020, 12, 3101 4 of 31 and realized a narrowband snow albedo simulation targeting satellite remote sensing observations. Our research results are expected to provide a basis for research on cryosphere change and even global climate change.

Methods
In the methodology section, we introduce the coupling principle of the integrated model. A new model based on snow hydrological processes and snow radiative transfer processes in cold regions is integrated to simulate the snow spectral albedo with high spectral resolution targeting satellite remote sensing observations. In the new integrated model, we couple the snow radiative transfer process based on the SNICAR model [32] with the snow hydrological process based on the GBEHM [29] and carefully track the transfer paths of solar radiation energy in multilayer snow cover; as a consequence, the simulation of the snow mass-energy balance process is more reasonable, and the limitations of single-model simulations are avoided.

Integration of Snow Energy Balance and Snow Radiative Transfer Processes
In the new integrated model, we couple the snow energy balance process with the snow radiative transfer process to track the transfer path of solar radiation energy in the snow layer. Assuming that the energy input absorbed by the surface snow layer is E sur ↓ , for the surface snow layer, the energy balance equation can be written as: where R snow,net ↑ is the net radiation flux of snow, E h ↓ is the sensible heat flux, E e ↓ is the latent heat flux, E p ↓ is the rainfall heating, and E g ↓ is the ground heat flow. The simulation of the equilibrium radiation of snow in the integrated model is based on the equilibrium radiation principle. Based on the equilibrium radiation equation, we can express the equilibrium radiation of snow as follows: R snow,net ↑ = R snow,direct + R snow,diffuse 1 − A snow,albedo − E snow,effective (2) Thus, the coupling between the snow energy balance process and the snow radiative transfer process in the integrated model can be expressed as: where R snow,direct is the direct solar radiation received from the snow surface, R snow,diffuse is the diffuse solar radiation received from the snow surface, A snow,albedo is the snow albedo, and E snow,effective is the effective radiation of snow. In the original GBEHM [29], the calculations of R snow,direct , R snow,diffuse and A snow,albedo on the right side of the equation all adopt the parameterized scheme of the BATS model without carefully considering the complex radiative transfer process inside the snow layer. In the new integrated model, on the basis of retaining the GBEHM to simulate the snow energy balance process, we introduce the SNICAR model with a more detailed description of the snow radiative transfer process to calculate R snow,direct , R snow,diffuse and A snow,albedo . The SNICAR model was initially established by Flanner [32] based on the WW model [45] and utilizes the two-stream radiative transfer solution of Toon [46].

Integration of the Snow Mass Balance and Snow Grain Size Evolution Processes
The snow grain size is the most important parameter when simulating snow albedo, and the snow grain size evolution is caused mainly by the compaction, densification, melting and sublimation of snow cover. In the new integrated model, the snow mass balance equation and the snow grain size evolution equation are combined through a specific functional relationship. The snow grain size and its growth rate are nontunable parameters, while the variations in the snow grain size and snow mass Remote Sens. 2020, 12, 3101 5 of 31 are partial derivatives with respect to time. The specific integration process is described below. In the GBEHM, the snow mass balance equation is: where θ i is the fraction volume of ice, θ l is the fraction volume of liquid, ρ i is the ice density, ρ l is the liquid density, U l is the liquid water flux, and z is the distance of a node from the ground surface.
In the new integrated model, we assume that the snow grain size evolution consists of four main processes of snow grain size change: those caused by dry snow, by the liquid water content, by new snowfall and by the recrystallization of liquid water. Assuming that the amount of snow grain size change due to snow mass change is ∂d ∂t , the snow grain size evolution process can be expressed as: where d is the snow grain size, t denotes time, ∆d dry,snow is the snow grain size change caused by dry snow, ∆d liquid,snow is the snow grain size change caused by the liquid water content, ∆d new f all,snow is the snow grain size change caused by new snowfall, and ∆d f reeze,snow is the snow grain size change caused by the recrystallization of liquid water. The dry snow evolution, liquid water change, new snowfall change and liquid water recrystallization are all simulated by the GBEHM in the snow grain size evolution process. The snow variables (soil and snow temperature, number of snow layers, liquid water content, ice content, snow fraction, snow depth, and snow density) that drive the evolution of the snow grain size are calculated by the snow hydrological module of the GBEHM. The change in snow grain size with temperature, precipitation and other microscopic parameters actually manifests as changes in the snow particle shape and volume, and a change in snow particle volume will lead to a change in snow mass. Therefore, snow grain growth is often closely related to the snow mass balance process. If the functional relationship between the snow mass change and snow grain size evolution is expressed by f (x), the integrated equation for the snow mass balance process and snow grain size evolution process can be expressed as follows:

Narrowband Snow Albedo Simulation Method Targeting MODIS Data
In this study, the spectral distribution characteristics of the incident solar flux ( Figure 1) suitable for typical mid-latitude winter atmospheric conditions are introduced based on the calculation of hyperspectral radiative transfer by Flanner et al. [32], and the snow spectral albedo, narrowband snow albedo and shortwave broadband snow albedo are simulated based on the weight coefficients of each band. To simulate the snow spectral albedo at each wavelength in the entire shortwave band in more detail, we divide the solar spectrum region (0.3-5.0 µm) into 470 spectral bands at a wavelength interval of 0.01 µm and simulate the snow spectral albedo based on the incident irradiance in each spectrum band. The weight coefficients of each band suitable for the calculation of snow spectral albedo in the integrated model are shown in Figure 1. According to the spectral ranges of the seven bands of MOD09GA reflectance data, a forward calculation method for narrowband snow albedo is proposed in this study for the MOD09GA reflectance data band. The core of the method is to calculate the narrowband snow albedo in a fixed band by utilizing both the wavelength-by-wavelength snow spectral albedo simulated by the integrated model and the incident solar flux spectral distribution function. The principle of the simulation method is as follows: For any wavelength segment between 1 λ and 2 λ , the equation for calculating the narrowband snow albedo observed by satellite remote sensing can be expressed as follows: where , snow albedo A is the narrowband snow albedo, i r is the snow spectral albedo for band i , and I ↓ is the incident solar radiation flux. By introducing the incident solar flux spectral distribution coefficients, when simulating the narrowband snow albedo targeting MODIS satellite remote sensing observations, we can simplify Equation (7) where N R is the narrowband snow albedo for band N simulated by the integrated model, i f represents the weights of the first i band applied to the spectral bands, which are specific to direct and diffuse cases, and m and n are the indexes corresponding to the starting and ending wavelengths of the narrowband snow albedo, respectively.

Parameter Settings and Transfer Process of the Integrated Model
In the model integration process, the snow cover is divided into a maximum of 5 snow layers according to the depth. For different snow layers, the multiple scattering of multilayer snow cover is considered. In the new integrated model, the VIS broadband snow albedo is defined in the wavelength range of 0.3-0.7 μm, the NIR broadband snow albedo is defined in the wavelength range of 0.7-3.0 μm, and the shortwave broadband snow albedo is defined in the wavelength range of 0.3-5.0 μm. The time step of the integrated model calculation is 1 h. The coupling framework of the integrated model is illustrated in Figure 2.
The parameter transfer process in the integrated model is as follows: 1. The beginning time loop initializes the snow status parameters, the atmospheric forcing data, the input parameters of the Mie scattering model, and the input parameters of the GBEHM. According to the spectral ranges of the seven bands of MOD09GA reflectance data, a forward calculation method for narrowband snow albedo is proposed in this study for the MOD09GA reflectance data band. The core of the method is to calculate the narrowband snow albedo in a fixed band by utilizing both the wavelength-by-wavelength snow spectral albedo simulated by the integrated model and the incident solar flux spectral distribution function. The principle of the simulation method is as follows: For any wavelength segment between λ 1 and λ 2 , the equation for calculating the narrowband snow albedo observed by satellite remote sensing can be expressed as follows: where A snow,albedo is the narrowband snow albedo, r i is the snow spectral albedo for band i, and I ↓ is the incident solar radiation flux. By introducing the incident solar flux spectral distribution coefficients, when simulating the narrowband snow albedo targeting MODIS satellite remote sensing observations, we can simplify Equation (7) to: [1,470], m, n ∈ N * ) (8) where R N is the narrowband snow albedo for band N simulated by the integrated model, f i represents the weights of the first i band applied to the spectral bands, which are specific to direct and diffuse cases, and m and n are the indexes corresponding to the starting and ending wavelengths of the narrowband snow albedo, respectively.

Parameter Settings and Transfer Process of the Integrated Model
In the model integration process, the snow cover is divided into a maximum of 5 snow layers according to the depth. For different snow layers, the multiple scattering of multilayer snow cover is considered. In the new integrated model, the VIS broadband snow albedo is defined in the wavelength range of 0.3-0.7 µm, the NIR broadband snow albedo is defined in the wavelength range of 0.7-3.0 µm, and the shortwave broadband snow albedo is defined in the wavelength range of 0.3-5.0 µm. The time step of the integrated model calculation is 1 h. The coupling framework of the integrated model is illustrated in Figure 2.

Conversion of MOD09GA Snow Reflectance into Broadband Snow Albedo
To verify the new integrated model and satellite remote sensing observations by using snow albedo data from ground observations, we use the asymptotic radiative transfer (ART) model to convert the MOD09GA snow reflectance into blue-sky snow albedo. The ART model, proposed in 2004 by Kokhanovsky [19], does not consider the effects of pollutants on snow albedo, but subsequent The parameter transfer process in the integrated model is as follows: 1.
The beginning time loop initializes the snow status parameters, the atmospheric forcing data, the input parameters of the Mie scattering model, and the input parameters of the GBEHM. Based on the snow mass balance and snow grain size evolution processes, the changes in snow grain size and other snow parameters are simulated dynamically. 4.
The snow grain size data simulated by the integrated model are combined with the snow energy balance and snow radiative transfer processes to simulate the transfer process of solar radiation energy in snow and track the energy changes caused by radiative scattering, absorption and reflection. 5.
Based on the simulation of the snow radiative transfer process by the integrated model, the snow spectral albedo in the solar spectrum region is estimated with a wavelength interval of 0.1 µm according to the incident solar flux spectral distribution. 6.
Aiming at the spectral waveband range of remote sensing satellite sensors, the incident solar flux spectral distribution and wavelength-by-wavelength snow spectral albedo are combined in the solar spectrum region to simulate the narrowband snow albedo targeting satellite remote sensing observations. 7.
The spatial and temporal variabilities in snow spectral albedo and narrowband snow albedo targeting MODIS observations are predicted based on the integrated model. The present loop is ended, and the next loop is entered.

Conversion of MOD09GA Snow Reflectance into Broadband Snow Albedo
To verify the new integrated model and satellite remote sensing observations by using snow albedo data from ground observations, we use the asymptotic radiative transfer (ART) model to convert the MOD09GA snow reflectance into blue-sky snow albedo. The ART model, proposed in 2004 by Kokhanovsky [19], does not consider the effects of pollutants on snow albedo, but subsequent researchers improved the ART model and proposed a new equation considering pollutant impacts [13]. The updated equation is as follows: where R is the snow reflectance, R 0 is the snow reflectance at no absorption, f is the angular function, g is the average cosine of the scattering angle, d is the snow grain size, λ is the wavelength, c is the volumetric concentration of grains, α is the bulk ice absorption coefficient, B is the absorption enhancement factor, κ is the absorption coefficient of polluted snow, and m is the absorption angstrom parameter. In this study, based on previous studies, we set the concentration of black carbon (BC) at 41.77 ± 6.36 ng g −1 [47].

Model Accuracy Validation
In this paper, to evaluate the simulation accuracy of the integrated model for the snow radiative transfer process, we compare the new model with the parameterized snow albedo scheme commonly used in the CoLM and BATS models. This parameterized scheme was originally proposed by Dickinson et al. [28] and considers two main variables, namely, the solar elevation angle and snow age, when simulating snow albedo.
The mean bias error (MBE), mean absolute error (MAE), root mean square error (RMSE), Pearson's correlation coefficient (R), coefficient of determination (R 2 ), and Nash-Sutcliffe efficiency Remote Sens. 2020, 12, 3101 9 of 31 coefficient (NSE) are used to assess the predictive power of the integrated model in simulating snow albedo [48]. These terms are defined as follows: where X i is the measured data, Y i is the modeled data, and X and Y are the average values of the measured data and modeled data, respectively.σ X and σ Y are the standard deviations of the measured data and modeled data, respectively. The kappa coefficient (K) is often used to check consistency and can also be used to measure the classification accuracy. To verify the spatial consistency between the narrowband snow albedo simulated by the integrated model targeting MODIS observations and the retrieved snow albedo (the MOD09GA data), the kappa coefficient is selected to evaluate their spatial consistency: where m is the total number of columns in the error matrix of remote sensing observations and the integrated model, x ii is the classification number of the confusion matrix composed of remote sensing observations and the integrated model simulation, and x ir and x ic are the total classification numbers. N is the total classification number used for the kappa coefficient accuracy evaluation.

Research Region
The Upstream Heihe River (UHR) basin is located on the Tibetan Plateau ( Figure 3) [10] at elevations in the range of 1637-5108 m [49]. The spatial extent of the study area is bounded by 98.57 • -101.17 • longitude and 37.72 • -39.12 • latitude. The distribution characteristics of snow cover in the UHR basin are characterized by instantaneous snow below 2700 m, patchy snow at 2700-3400 m and permanent snow above 3400 m [50]. Snowfall is more frequent in spring; spring and winter exhibit higher average snow albedo, whereas summer and autumn are characterized by lower average snow albedo [49]. The number of annual snow-covered days in most areas of the research region exceeds 60 days [51]. Thus, this area is an ideal place to study snow albedo. In addition, we have established two automatic weather stations in the study area, which can conduct extended snow observations and provide a large amount of observational data for snow radiative transfer process simulations over long time series and at a large spatial scale for verifying the authenticity of the proposed model. higher average snow albedo, whereas summer and autumn are characterized by lower average snow albedo [49]. The number of annual snow-covered days in most areas of the research region exceeds 60 days [51]. Thus, this area is an ideal place to study snow albedo. In addition, we have established two automatic weather stations in the study area, which can conduct extended snow observations and provide a large amount of observational data for snow radiative transfer process simulations over long time series and at a large spatial scale for verifying the authenticity of the proposed model.  The observation data from Yakou and Jingyangling stations in areas characterized by long snow days and large tracts of snow cover are selected in this research as the verification data for the integrated model [52,53]. The elevations of Yakou and Jingyangling stations exceed 3700 m, annual snowfall is abundant, and the underlying surface is meadow or grassland, although the surface type has little effect on the snow radiative transfer simulation; thus, these stations can provide stable observation data in the long term ( Figure 3).

Data
Multiple forms of data are used to drive and validate the integrated model, including remote sensing data, ground observation data, Weather Research and Forecasting (WRF) meteorological data, snow optical characteristics data and other data. The specific data types, attributes and formats are shown in Table 1. In previous studies, we found that the accuracy of MOD10A1 snow albedo data and MCD43 snow albedo data is poor over the Tibetan Plateau [49]. Therefore, it is difficult to fully verify the new integrated model simulation results when comparing MOD10A1 snow albedo data and MCD43 snow albedo data with the simulation results of the new integrated model. In contrast, we previously suggested that the snow albedo retrieved from MOD09GA reflectance data is highly accurate [49]. Hence, the MOD09GA reflectance product contains ideal snow albedo data because of its high temporal resolution. In addition, the research region is located in the northeastern part of the Tibetan Plateau, where the underlying surface is basically meadow and exhibits good land surface uniformity. Therefore, in the research region, the MOD09GA snow reflectance is relatively uniform on the pixel scale and is less affected by the presence of mixed pixels. For these reasons, we initially prefer MOD09GA reflectance data for simulating the snow radiative transfer process targeting satellite remote sensing observations. When selecting MOD09GA snow reflectance data, we refer to quality assurance (QA) information and extract only high-quality MOD09GA snow reflectance data without cloud cover. However, the MODIS QA data have a poor ability to recognize snow cover. Therefore, we combine MOD10A1 fractional snow cover (FSC) data with MOD10A1 normalized difference snow index (NDSI) data to obtain the most ideal pure snow pixels as possible. The MOD10A1 dataset is from the National Snow and Ice Data Center (NSIDC), and the MOD09GA reflectance dataset is from the National Aeronautics and Space Administration (NASA). The MOD09GA reflectance products provide 500 m land surface spectral reflectance values and 1 km angular grid values. The time-resolved MOD10A1 dataset and MOD09GA reflectance dataset are daily, and the spatial resolution is 500 m.

Ground Observation Data
The ground observation data mainly include upwelling shortwave radiation flux (USRF) and downward shortwave radiation flux (DSRF). In this study, the ground observation data are from Yakou and Jingyangling snow stations. The site verification data mainly include snow radiation data and snow characteristics, and the observation data employed herein are from days with snow cover at the station. The observation interval of the site data is 30 min. Among these observations, the snow reflectance data can be calculated from the albedo meter [49].

Meteorological Data
The meteorological data are produced by the WRF model [54]. The dataset is provided by the CARSDC. This dataset includes hourly temperature, wind speed, precipitation, relative humidity, atmospheric pressure, and longwave and shortwave radiation data [55]. To match the spatial resolution of the new integrated model, we resample these data to 1 km.
We verify the DSRF data in the WRF-simulated atmospheric forcing data using DSRF observations from Yakou station ( Figure 4). The MBE, RMSE and R are 23 W/m 2 , 98 W/m 2 and 0.90, respectively, between the simulated and observed DSRF. Furthermore, to prove the effectiveness of the radiation data driving the integrated model, we quote the verification results from the authors of the original data at three automatic weather stations in the UHR basin. The Pan et al. [56]  results show that at Guantan station, the MBE, RMSE and R are 31.5 W/m 2 , 89.7 W/m 2 and 0.92, respectively, between the simulated and observed DSRF; at Huazhaizi station, the MBE, RMSE and R are 35.5 W/m 2 , 160.8 W/m 2 and 0.88, respectively; and at Maliantan station, the MBE, RMSE and R are 47.8 W/m 2 , 174.0 W/m 2 and 0.86, respectively. The above verification results all confirm that the WRF-simulated DSRF data can reflect the DSRF changes in the UHR basin and can be used to drive the data of the climate models.

Other Data
The other data mainly include soil data, Digital Elevation Model (DEM) data and land use data. The other data are from the CARSDC. These data are a Chinese subset of global land cover data supported by the International Geosphere-Biosphere Program Data and Information System (IGBP-DIS) based on Advanced Very High Resolution Radiometer (AVHRR) remote sensing data; the data are referred to as IGBP-DIS [57]. The IGBP-DIS data use the United States Geological Survey (USGS) classification method, which adopts the IGBP classification system that divides the Earth into 17 categories with continental units. The soil dataset is provided by the CARSDC [58]. The soil data are in a grid format, and the projection is WGS84. The DEM data are from the Shuttle Radar Topography Mission (SRTM) with a 90 m resolution. The land use data, soil data and DEM data are all resampled to 1 km.

Snow Optical Characteristics Data
The snow optical characteristics data mainly include snow and aerosol Mie parameter data, snow grain size evolution lookup table data, and light-absorbing snow impurity lookup table data. The snow optical characteristics data are from the CESM. The Mie parameter data contain the singlescattering albedo at each wavelength, asymmetry factor, scattering efficiency, extinction efficiency, absorption efficiency, scattering phase function and other Mie parameter data. The snow grain size evolution lookup table data mainly include the lookup table data involving the snow grain size parameterized scheme. The light-absorbing snow impurity lookup table data mainly include BC, dust and other pollutant lookup table data.

Other Data
The other data mainly include soil data, Digital Elevation Model (DEM) data and land use data. The other data are from the CARSDC. These data are a Chinese subset of global land cover data supported by the International Geosphere-Biosphere Program Data and Information System (IGBP-DIS) based on Advanced Very High Resolution Radiometer (AVHRR) remote sensing data; the data are referred to as IGBP-DIS [57]. The IGBP-DIS data use the United States Geological Survey (USGS) classification method, which adopts the IGBP classification system that divides the Earth into 17 categories with continental units. The soil dataset is provided by the CARSDC [58]. The soil data are in a grid format, and the projection is WGS84. The DEM data are from the Shuttle Radar Topography Mission (SRTM) with a 90 m resolution. The land use data, soil data and DEM data are all resampled to 1 km.

Snow Optical Characteristics Data
The snow optical characteristics data mainly include snow and aerosol Mie parameter data, snow grain size evolution lookup table data, and light-absorbing snow impurity lookup table data. The snow optical characteristics data are from the CESM. The Mie parameter data contain the single-scattering albedo at each wavelength, asymmetry factor, scattering efficiency, extinction efficiency, absorption efficiency, scattering phase function and other Mie parameter data. The snow grain size evolution lookup table data mainly include the lookup table data involving the snow grain size parameterized scheme. The light-absorbing snow impurity lookup table data mainly include BC, dust and other pollutant lookup table data.

Accuracy Verification of the Integrated Model
Based on the simulation of the integrated model at a single point, we use the measured broadband snow albedo from Yakou and Jingyangling stations in this study to verify and evaluate the simulation results of the parameterized scheme and integrated model. Based on the verification results of the measured data, we compare the simulation accuracy of the parameterized scheme with that of the integrated model scheme.
At Yakou station, the MAE, RMSE, R, R 2 and NSE are 0.10, 0.11, 0.27, 0.07 and −1.42, respectively, between the parameterized blue-sky snow albedo and measured snow albedo, while those between the blue-sky snow albedo simulated by the new integrated model and the measured snow albedo are 0.03, 0.04, 0.85, 0.73 and 0.72, respectively (Table 2). At Jingyangling station, the MAE, RMSE, R, R 2 and NSE are 0.14, 0.16, 0.27, 0.07 and −1.05, respectively, between the parameterized blue-sky snow albedo and measured snow albedo and 0.04, 0.06, 0.80, 0.64 and 0.40, respectively, between the blue-sky snow albedo simulated by the new integrated model and the measured snow albedo (Table 2). At Yakou station, the estimated deviation of the parameterized blue-sky snow albedo is 10%, while that of the blue-sky snow albedo simulated by the new integrated model is only 3%. At Jingyangling station, the estimated deviation of the parameterized blue-sky snow albedo is 14%, while that of the blue-sky snow albedo simulated by the new integrated model is only 4%. In general, the new integrated model based on physical mechanisms greatly reduces the error in traditional radiation simulations. Moreover, the blue-sky snow albedo simulation in the new integrated model is more reasonable and closer to the measured results than the parameterized blue-sky albedo ( Figure 5). Consequently, the integrated model is more reasonable for simulating the dynamic variation trends of blue-sky snow albedo both spatially and temporally and provides the possibility of simulating fine snow spectral albedo.
results of the measured data, we compare the simulation accuracy of the parameterized scheme with that of the integrated model scheme.
At Yakou station, the MAE, RMSE, R, R 2 and NSE are 0.10, 0.11, 0.27, 0.07 and −1.42, respectively, between the parameterized blue-sky snow albedo and measured snow albedo, while those between the blue-sky snow albedo simulated by the new integrated model and the measured snow albedo are 0.03, 0.04, 0.85, 0.73 and 0.72, respectively (Table 2). At Jingyangling station, the MAE, RMSE, R, R 2 and NSE are 0.14, 0.16, 0.27, 0.07 and −1.05, respectively, between the parameterized blue-sky snow albedo and measured snow albedo and 0.04, 0.06, 0.80, 0.64 and 0.40, respectively, between the bluesky snow albedo simulated by the new integrated model and the measured snow albedo (Table 2). At Yakou station, the estimated deviation of the parameterized blue-sky snow albedo is 10%, while that of the blue-sky snow albedo simulated by the new integrated model is only 3%. At Jingyangling station, the estimated deviation of the parameterized blue-sky snow albedo is 14%, while that of the blue-sky snow albedo simulated by the new integrated model is only 4%. In general, the new integrated model based on physical mechanisms greatly reduces the error in traditional radiation simulations. Moreover, the blue-sky snow albedo simulation in the new integrated model is more reasonable and closer to the measured results than the parameterized blue-sky albedo ( Figure 5). Consequently, the integrated model is more reasonable for simulating the dynamic variation trends of blue-sky snow albedo both spatially and temporally and provides the possibility of simulating fine snow spectral albedo.  After validating the measured data from the snow stations, we find that the verification result of Yakou station is higher than that of Jingyangling station. The main reason for this result is that the snow observation stations are located on the Tibetan Plateau, and the scale of the plateau and thus the spatial representativeness of site observations often lead to differences in the verification results [49,59]. The estimation ability of the integrated model for snow surface reflectance information is higher than that of the parameterized scheme, which poorly simulates (seriously underestimates) the real snow albedo. Moreover, the integrated scheme, which combines multisource data and multiple models, has a greater accuracy than the traditional parameterized scheme in simulating blue-sky snow albedo; in other words, the integrated scheme effectively resolves the blue-sky snow albedo underestimation problem and achieves a superior blue-sky snow albedo simulation accuracy. Finally, the simulation accuracy of large-scale snow cover radiation is often poor, but with the new integrated method, the simulation accuracy of large-scale snow radiative transfer processes is greatly improved.

Simulation of Spatiotemporally Distributed Snow Spectral Albedo
The new integrated model based on the snow hydrological and snow radiative transfer processes in cold regions can simulate snow reflectance with a spectral resolution of 1 nm involving 470 bands in the solar spectrum region (in this section, as an example, the snow spectral albedo values of 20 single wavelength bands are between 0.6 and 2.5 µm); snow spectral reflectance simulations of this magnitude will undoubtedly enhance the monitoring of snow radiation characteristics. In addition, the new integrated model achieves completely continuous time series.
The snow spectral albedo simulated by the integrated model exceeds 0.6 in the VIS band ( Figure 6). After entering the NIR band, the snow spectral albedo decreases sharply, shows a trough at a wavelength of 1 µm and then continues to decline after a slight increase. The snow spectral albedo decreases to the lowest level at wavelengths between 1.5 and 1.6 µm. When the wavelength is greater than 1.6 µm, the snow spectral albedo increases slightly. Above 1.9-2.0 µm, the snow spectral albedo continues to fluctuate and decrease. Comparing the snow spectral variation curves ( Figure A1; details on the variation characteristics of the snow spectrum albedo are listed in Appendix A) simulated by the integrated model, the variation trends of the spatially distributed snow spectral albedo based on the integrated model simulation and the snow spectral reflectance curve predicted by the snow radiative transfer model are exactly the same in the shortwave band. Likewise, the variation trends of the snow spectral variation curve simulated by the integrated model and the snow spectral albedo curve measured indoors by the U.S. Army Cold Regions Research and Engineering Laboratory (CRREL) are also completely consistent [60] (Figure A2; details on the variations in the snow spectral albedo with wavelength simulated by CRREL are listed in Appendix B), and these trends respond to the snow spectral albedo changes in the shortwave band. In contrast to many previous studies on snow radiative transfer simulation, we simulate the spectral albedo for the whole solar spectrum region, which greatly expands the spectral range in the forward simulation. The real-time characteristics of land surface-driven data and site observation data are used to simulate the real-time changes in the snow spectral signature. The integrated model is important not only for simulating changes in snow spectral albedo but also for predicting the snow spectral albedo in solar spectrum regions based on meteorological forcing data. In the new integrated model, we add a snow spectral albedo simulation to the output parameters of the traditional snow radiative transfer process simulation and directly simulate the snow spectral albedo from the ground. As shown in Figure 7, when the meteorological forcing data are known, the new integrated model can predict the snow spectral albedo at any wavelength in any time period. The integrated model successfully predicts the spatial and temporal evolution trends of the snow spectral albedo at various wavelengths in the solar spectrum region. In addition, Figure 7 shows that the snow spectral albedo change trend in the NIR band is more significant than that in the VIS band, and these simulation results are supported by the fact that the snow spectral reflectance in the NIR band with a wavelength of 0.7-1.4 μm has a large downward trend, as shown in Figure A1 (Appendix A). In contrast to the traditional snow spectral curve simulation, the spatiotemporally distributed snow spectral albedo simulation greatly refines the spatial and temporal distributions and changes the characteristics of the snow spectrum. Refinement of the snow spectral albedo simulation process will greatly improve simulations of surface radiation balance and energy exchange processes and thus significantly affect future simulations and predictions of climate change [61]. The integrated model is important not only for simulating changes in snow spectral albedo but also for predicting the snow spectral albedo in solar spectrum regions based on meteorological forcing data. In the new integrated model, we add a snow spectral albedo simulation to the output parameters of the traditional snow radiative transfer process simulation and directly simulate the snow spectral albedo from the ground. As shown in Figure 7, when the meteorological forcing data are known, the new integrated model can predict the snow spectral albedo at any wavelength in any time period. The integrated model successfully predicts the spatial and temporal evolution trends of the snow spectral albedo at various wavelengths in the solar spectrum region. In addition, Figure 7 shows that the snow spectral albedo change trend in the NIR band is more significant than that in the VIS band, and these simulation results are supported by the fact that the snow spectral reflectance in the NIR band with a wavelength of 0.7-1.4 µm has a large downward trend, as shown in Figure A1 (Appendix A). In contrast to the traditional snow spectral curve simulation, the spatiotemporally distributed snow spectral albedo simulation greatly refines the spatial and temporal distributions and changes the characteristics of the snow spectrum. Refinement of the snow spectral albedo simulation process will greatly improve simulations of surface radiation balance and energy exchange processes and thus significantly affect future simulations and predictions of climate change [61].

Narrowband Snow Albedo Simulation Targeting MODIS Sensors
Based on the new integrated model, the narrowband snow albedo is simulated for the 7 reflectance bands of MOD09GA data, and the results are compared with the narrowband snow albedo retrieved from MOD09GA snow reflectance data (Figure 8). To maintain the spatiotemporal consistency between the remote sensing observations and model simulations, we choose the simulation results from the new integrated model at the same time as the remote sensing observation time. The overall accuracy and kappa coefficient between the narrowband snow albedo targeting MOD09GA snow reflectance data simulated by the integrated model and the narrowband snow albedo retrieved by utilizing the MOD09GA snow reflectance data are 0.81 and 0.5, respectively. Thus, the narrowband snow albedo simulated by the integrated model is basically consistent with the narrowband snow albedo retrieved by satellite remote sensing in each band. In the VIS band, the modeled results based on the integrated model are greater than the retrieved results based on satellite remote sensing, and both are more consistent in the NIR band than in the VIS band. The simulation results of the integrated model show that Band 1-Band 4 corresponding to the MOD09GA reflectance data are located in the VIS band with a high narrowband albedo, while Band 5-Band 7 corresponding to the MOD09GA reflectance data are located in the NIR band with a low narrowband albedo. Overall, these simulation results agree well with the snow surface albedo retrieved by the MOD09GA reflectance products.

Narrowband Snow Albedo Simulation Targeting MODIS Sensors
Based on the new integrated model, the narrowband snow albedo is simulated for the 7 reflectance bands of MOD09GA data, and the results are compared with the narrowband snow albedo retrieved from MOD09GA snow reflectance data (Figure 8). To maintain the spatiotemporal consistency between the remote sensing observations and model simulations, we choose the simulation results from the new integrated model at the same time as the remote sensing observation time. The overall accuracy and kappa coefficient between the narrowband snow albedo targeting MOD09GA snow reflectance data simulated by the integrated model and the narrowband snow albedo retrieved by utilizing the MOD09GA snow reflectance data are 0.81 and 0.5, respectively. Thus, the narrowband snow albedo simulated by the integrated model is basically consistent with the narrowband snow albedo retrieved by satellite remote sensing in each band. In the VIS band, the modeled results based on the integrated model are greater than the retrieved results based on satellite remote sensing, and both are more consistent in the NIR band than in the VIS band. The simulation results of the integrated model show that Band 1-Band 4 corresponding to the MOD09GA reflectance data are located in the VIS band with a high narrowband albedo, while Band 5-Band 7 corresponding to the MOD09GA reflectance data are located in the NIR band with a low narrowband albedo. Overall, these simulation results agree well with the snow surface albedo retrieved by the MOD09GA reflectance products. Overall, the integrated model successfully simulates the spatial distribution characteristics of the narrowband snow albedo targeting the MOD09GA reflectance data and the differences in narrowband snow albedo among different spectral bands. More importantly, the integrated model simulation process completely avoids the influence of clouds due to the controllability of the driving data and can effectively simulate continuous changes in snow surface reflectance information both spatially and temporally.
The integrated model can simulate and predict the spatially distributed narrowband snow albedo with an arbitrary time series (Figure 9). In the new integrated model, radiation monitoring images at satellite transit times can be simulated, while radiation monitoring images for periods without remote sensing data can be predicted, which compensates for the deficiencies of snow albedo observations from satellite remote sensing in the time series. As shown in Figure 9, we can utilize the new integrated model to predict the narrowband snow albedo obtained from MOD09GA reflectance data at any time. Overall, the integrated model successfully simulates the spatial distribution characteristics of the narrowband snow albedo targeting the MOD09GA reflectance data and the differences in narrowband snow albedo among different spectral bands. More importantly, the integrated model simulation process completely avoids the influence of clouds due to the controllability of the driving data and can effectively simulate continuous changes in snow surface reflectance information both spatially and temporally.
The integrated model can simulate and predict the spatially distributed narrowband snow albedo with an arbitrary time series (Figure 9). In the new integrated model, radiation monitoring images at satellite transit times can be simulated, while radiation monitoring images for periods without remote sensing data can be predicted, which compensates for the deficiencies of snow albedo observations from satellite remote sensing in the time series. As shown in Figure 9, we can utilize the new integrated model to predict the narrowband snow albedo obtained from MOD09GA reflectance data at any time. By validating the results of the integrated model and remote sensing observations in the same period ( Figure 10, Table 3), we find that the average narrowband snow albedo in 7 bands simulated by the integrated model is 0.584 and that the average narrowband snow albedo in 7 bands retrieved by utilizing the MOD09GA reflectance data is 0.564. In addition, the average MAE and RMSE between the narrowband snow albedo in 7 bands simulated by the integrated model and the narrowband snow albedo in 7 bands retrieved by utilizing the MOD09GA reflectance data are 0.024 and 0.027, respectively. The simulation errors in Band 1 and Band 7 are the smallest, and the MAE and RMSE in these two bands are 0.011 and 0.014, respectively; in contrast, Band 4 has the largest simulation error with an MAE of 0.053 and an RMSE of 0.054 between the modeled Band 4 and remotely sensed Band 4. Overall, the average narrowband snow albedo in 7 bands simulated by the integrated model is slightly larger than that retrieved by utilizing the MOD09GA reflectance data. However, the narrowband snow albedo simulated by the integrated model has a small deviation in the time series from that retrieved by using the MOD09GA reflectance data. Therefore, the integrated By validating the results of the integrated model and remote sensing observations in the same period ( Figure 10, Table 3), we find that the average narrowband snow albedo in 7 bands simulated by the integrated model is 0.584 and that the average narrowband snow albedo in 7 bands retrieved by utilizing the MOD09GA reflectance data is 0.564. In addition, the average MAE and RMSE between the narrowband snow albedo in 7 bands simulated by the integrated model and the narrowband snow albedo in 7 bands retrieved by utilizing the MOD09GA reflectance data are 0.024 and 0.027, respectively. The simulation errors in Band 1 and Band 7 are the smallest, and the MAE and RMSE in these two bands are 0.011 and 0.014, respectively; in contrast, Band 4 has the largest simulation error with an MAE of 0.053 and an RMSE of 0.054 between the modeled Band 4 and remotely sensed Band 4. Overall, the average narrowband snow albedo in 7 bands simulated by the integrated model is slightly larger than that retrieved by utilizing the MOD09GA reflectance data. However, the narrowband snow albedo simulated by the integrated model has a small deviation in the time series from that retrieved by using the MOD09GA reflectance data. Therefore, the integrated model is feasible for simulating the narrowband snow albedo targeting the MOD09GA reflectance data in the time series. Therefore, to a certain extent, the integrated model successfully predicts the observation results in 7 bands of MOD09GA reflectance data in continuous time series.
Remote Sens. 2020, 12, x FOR PEER REVIEW 20 of 31 model is feasible for simulating the narrowband snow albedo targeting the MOD09GA reflectance data in the time series. Therefore, to a certain extent, the integrated model successfully predicts the observation results in 7 bands of MOD09GA reflectance data in continuous time series.

Discussion
In a traditional land surface model, the snow radiative transfer process is not strictly considered, and even if this process is considered, the detailed transfer process of solar radiation energy in the snow layer is still not explained in detail. In addition, traditional land surface models do not take into account the variation in snow mass caused by changes in the snow grain size with meteorological factors. In view of these characteristics, combined with the basic mass-energy balance principle, we introduce a new method by which to integrate snow energy balance and snow radiative transfer processes as well as snow mass balance and snow grain size evolution processes. In addition, we also propose a simplified method for simulating narrowband snow albedo and snow spectral albedo targeting mainstream satellite remote sensing observations. The first and most important improvement is that the new integrated model extends the simulation range of snow spectral albedo

Discussion
In a traditional land surface model, the snow radiative transfer process is not strictly considered, and even if this process is considered, the detailed transfer process of solar radiation energy in the snow layer is still not explained in detail. In addition, traditional land surface models do not take into account the variation in snow mass caused by changes in the snow grain size with meteorological factors. In view of these characteristics, combined with the basic mass-energy balance principle, we introduce a new method by which to integrate snow energy balance and snow radiative transfer processes as well as snow mass balance and snow grain size evolution processes. In addition, we also propose a simplified method for simulating narrowband snow albedo and snow spectral albedo targeting mainstream satellite remote sensing observations. The first and most important improvement is that the new integrated model extends the simulation range of snow spectral albedo and can simulate the spatiotemporal evolution characteristics of snow spectral albedo in the whole shortwave band, effectively expanding the solar radiation energy simulation capacity for basin-scale snow radiation balance processes. Second, the new integrated model can be calibrated using snow spectral signature data from satellite remote sensing observations, expanding the scope of traditional hydrological and meteorological observations. Third, the new integrated model fully utilizes ground data to simulate satellite remote sensing observations and can simulate and predict narrowband snow albedo targeting observations from specific satellite sensors according to their band ranges. Fourth, this new method is a significant improvement over the direct simulation of satellite remote sensing observations based on land surface models.

Improved Direct Simulation of Snow Spectral Albedo by the New Method
In the study of land surface processes and climate change, high spectral resolution snow radiation information usually provides the following advantages. For different underlying surfaces, subtle ground object features greatly influence the albedo of thin snow. However, due to the fixed bandwidth and coarse spectral resolution of snow radiation information observed by traditional simulation methods or mainstream remote sensing satellites, it is difficult to distinguish and characterize the variations in the albedo of thin snow among different underlying surfaces, eventually resulting in incorrect snow radiation estimates in research on land surface processes and climate change. In addition, pollutants have a considerable impact on snow albedo; furthermore, different types of pollutants have remarkably varying influences on the snow spectral signature, and the effect of the same pollutant on snow radiation in different spectral ranges also varies notably. In conclusion, it is necessary to obtain snow radiation information with a high spectral resolution. High spectral resolution snow spectral albedo information can be used to not only distinguish the radiation conditions of thin snow among different underlying surfaces but also effectively describe the influences of pollutants on the snow spectral signature within a fine spectral band interval, thus improving the estimation accuracy of snow radiation information in research on land surface processes and climate change.
For many studies on land surface processes and climate change, it is often necessary to obtain snow reflectance and snow albedo data with a high spectral resolution and continuous spectral band [4,62]. However, current optical remote sensing satellites are unable to satisfy these requirements; nevertheless, they can perform snow radiation monitoring over the whole shortwave band. Additionally, satellite remote sensing observations can provide snow reflectance and snow albedo data only in the corresponding sensor bands and thus cannot obtain spectral information within a continuous spectral band. Consequently, direct observations from satellite remote sensing are difficult to employ as traditional meteorological and hydrological information (such as runoff, soil temperature and humidity, and snow depth) in the testing and calibration of land surface models. As a result, the snow radiation information observed by remote sensing is vastly wasted.
We solve this problem by coupling multiple models. That is, the new integrated model can directly predict snow reflectance and snow albedo by targeting satellite remote sensing observations with a high spectral resolution and a continuous spectral band while covering the whole shortwave band. Furthermore, the new integrated model can simulate the snow spectral albedo within any wavelength interval of the whole shortwave band (wavelength range: 300-5000 nm). In contrast, mainstream remote sensing satellites (here, we take the optical remote sensing satellites MODIS and TM/ETM+ as examples; in Figure 11, the blue band represents MODIS, while the purple band represents TM/ETM+) can obtain snow spectral albedo only within a fixed bandwidth. In addition, the SNICAR model obtains snow albedo only in the VIS and NIR bands ( Figure 11). Remote Sens. 2020, 12, x FOR PEER REVIEW 22 of 31 Due to the limitations of traditional model simulations, there are considerable uncertainties in the acquisition, simulation and calibration of snow albedo using traditional models to simulate remote sensing satellite observations. In general, several transformations need to be applied to remote sensing snow spectral albedo observations to yield the snow albedo simulated by a traditional model. However, many intermediate transformations may correspond to a large number of linear transformations, and the need for many intermediate steps generates the substantial accumulation of errors, introducing remarkable uncertainty. For example, the traditional simulation method requires three steps to transform the ground-observed snow reflectance into the remotely sensed snow reflectance.
Step one uses a radiative transfer model to simulate the land surface snow reflectance; step two transforms the land surface snow reflectance into the land surface snow albedo; and step three transforms the land surface snow albedo into the remotely sensed snow reflectance. Different from the traditional simulation method, our new integrated model requires only one step to convert ground observations into remote sensing observations, thereby reducing the conversion error ( Figure  12). In other words, the new method can directly simulate the snow reflectance of remote sensing observations. The proposed model simulation process simplifies the intermediate transformation process and reduces the impact of intermediate transformations on the simulation accuracy. Hence, the new method directly links the ground-simulated albedo with remote sensing observations and expands the application of snow radiation information from remote sensing observations in land surface models. Due to the limitations of traditional model simulations, there are considerable uncertainties in the acquisition, simulation and calibration of snow albedo using traditional models to simulate remote sensing satellite observations. In general, several transformations need to be applied to remote sensing snow spectral albedo observations to yield the snow albedo simulated by a traditional model. However, many intermediate transformations may correspond to a large number of linear transformations, and the need for many intermediate steps generates the substantial accumulation of errors, introducing remarkable uncertainty. For example, the traditional simulation method requires three steps to transform the ground-observed snow reflectance into the remotely sensed snow reflectance.
Step one uses a radiative transfer model to simulate the land surface snow reflectance; step two transforms the land surface snow reflectance into the land surface snow albedo; and step three transforms the land surface snow albedo into the remotely sensed snow reflectance. Different from the traditional simulation method, our new integrated model requires only one step to convert ground observations into remote sensing observations, thereby reducing the conversion error ( Figure 12). In other words, the new method can directly simulate the snow reflectance of remote sensing observations. The

Role of the Snow Spectral Albedo Simulation with a High Spectral Resolution
Generally, on the macroscopic scale, the accumulation, melting and refreezing of snow are all reflected in snow spectral albedo changes [26]. On the microscopic scale, the snow grain shape, snow effective grain size evolution, types of pollutants in snow, concentrations of pollutants in snow, roughness of the snow surface, and water content of snow all affect the snow spectral albedo changes [63][64][65][66]. Therefore, the simulation and prediction of snow spectral changes can help to better understand the macroscopic and microscopic changes in snow cover. Moreover, spatiotemporally distributed snow spectral albedo simulations can often indirectly reflect changes in snow accumulation, melting and refreezing and can also indirectly reflect the spatiotemporal evolution characteristics of snow cover.
In addition, high spectral resolution snow albedo information simulated by land surface models targeting satellite remote sensing observations will contribute to better describing the radiation conditions of thin snow on different underlying surfaces. Moreover, high spectral resolution snow albedo information targeting satellite remote sensing observations can better reflect the snow spectral signatures in different states, with different properties, of different kinds, and of polluted snow. Providing detailed descriptions of snow spectral signatures will contribute to accurately estimating snow radiation information in research on land surface processes and climate change. Finally, high spectral resolution snow spectral albedo information simulated based on the new integrated model can compensate for the shortcomings of satellite remote sensing observations in terms of time series and spectral range.

Limitations and Uncertainties of the New Integrated Model
In Figure 5, in comparison with the blue-sky snow albedo data (integrated, measured and parameterized), the measured albedo data show a tendency to decay beginning approximately on day 50. In contrast, the prediction results of the integrated model show an increasing trend. An analysis holds that the main causes of this error can be attributed to a data problem and a model

Role of the Snow Spectral Albedo Simulation with a High Spectral Resolution
Generally, on the macroscopic scale, the accumulation, melting and refreezing of snow are all reflected in snow spectral albedo changes [26]. On the microscopic scale, the snow grain shape, snow effective grain size evolution, types of pollutants in snow, concentrations of pollutants in snow, roughness of the snow surface, and water content of snow all affect the snow spectral albedo changes [63][64][65][66]. Therefore, the simulation and prediction of snow spectral changes can help to better understand the macroscopic and microscopic changes in snow cover. Moreover, spatiotemporally distributed snow spectral albedo simulations can often indirectly reflect changes in snow accumulation, melting and refreezing and can also indirectly reflect the spatiotemporal evolution characteristics of snow cover.
In addition, high spectral resolution snow albedo information simulated by land surface models targeting satellite remote sensing observations will contribute to better describing the radiation conditions of thin snow on different underlying surfaces. Moreover, high spectral resolution snow albedo information targeting satellite remote sensing observations can better reflect the snow spectral signatures in different states, with different properties, of different kinds, and of polluted snow. Providing detailed descriptions of snow spectral signatures will contribute to accurately estimating snow radiation information in research on land surface processes and climate change. Finally, high spectral resolution snow spectral albedo information simulated based on the new integrated model can compensate for the shortcomings of satellite remote sensing observations in terms of time series and spectral range.

Limitations and Uncertainties of the New Integrated Model
In Figure 5, in comparison with the blue-sky snow albedo data (integrated, measured and parameterized), the measured albedo data show a tendency to decay beginning approximately on day 50. In contrast, the prediction results of the integrated model show an increasing trend. An analysis holds that the main causes of this error can be attributed to a data problem and a model coupling problem. The former is reflected mainly in the WRF datasets, which have certain errors (Figure 4) [67], and thus, the data inevitably accumulate errors when driving the integrated model. The latter is reflected primarily in the differences between the descriptions of the snow cover state by different models, which inevitably lead to the accumulation of errors when integrating multiple models. Furthermore, there are some errors in the simulation of the fractional snow cover, snow cover extent and snow depth in land surface process models [68,69], which lead to deviations in the snow state and snow distribution simulations and eventually deviations in the snow albedo simulation. In addition, the scale difference between the snow station measurements and spatial resolution of the integrated model data and the input parameter uncertainties of various models both result in an overestimation of the simulation results. Generally, the above simulation errors are allowable in the simulation of snow radiative transfer processes at large spatial scales and over long time series. Although our integrated model refines most of the physical snow radiative transfer processes, the model still contains some parameterization details that are not depicted by physical models, which is one of the reasons for error.
Many studies have shown that snow albedo simulations are affected not only by the snow age, solar zenith angle, snow grain size, snow cover density, snow depth and snow water equivalent but also by light-absorbing snow impurities [70][71][72], which can reduce snow albedo [70,73]. Due to the limitations of observational data, light-absorbing snow impurities are only roughly considered in this study, which is another main reason for the inconsistency between the variation trends of the simulated and observed snow albedo. The topographic variables are an important factor affecting the snow radiative transfer simulation process in complex environments. Unfortunately, we have not yet found an effective way to solve this problem in the new integrated model. This is also a defect of the new integrated model. In future studies, combined with DEM, slope, aspect and other topographic variables to develop a snow radiative transfer model suitable for complex terrain is expected to solve this problem.
Although the simulation accuracy of the new integrated model has been comprehensively verified, due to the lack of effective observation data, the new integrated model remains unproven during the spring snowmelt period and autumn snow cover onset period. A new universal method should be tested using data from much more locations. However, in our new integrated model, driving the whole integrated model requires high-precision atmospheric forcing datasets, remote sensing data and other ground observation data. Therefore, the new integrated model lacks accuracy assessment in the typical snow-covered regions such as Antarctica and Arctic. Hence, there is still some uncertainty regarding the accuracy of the new integrated model. In future studies, we will further expand the verification area of the new integrated model, and add simulations and verification experiments in the typical snow-covered regions to prove the universality of the new integrated model.

Conclusions
In this study, a new integrated model scheme is introduced to simulate snow surface spectral reflections and narrowband snow albedo targeting satellite remote sensing observations. In the new integrated model, a snow radiative transfer module is coupled with a snow hydrological model and a land surface process model to simulate snow radiative transfer processes at the basin scale and to carefully track the transfer process of solar radiation energy in snow. Utilizing the new scheme, the snow spectral albedo simulation in the solar spectrum region is refined, and the land surface narrowband snow albedo simulation targeting satellite remote sensing observations is realized based on meteorological forcing data and snow optical data. The broadband snow albedo simulated by the integrated model agrees well with satellite remote sensing observations.
The variation trends of the spatiotemporally distributed snow spectral albedo simulated by the new integrated model and the experimentally measured snow spectral albedo exhibit good consistency in the solar spectrum region, and the simulation results can significantly reflect the temporal and spatial variation characteristics of land surface snow spectral albedo. Based on the new integrated model, we extend the simulation of snow spectral albedo to the whole solar spectrum region and refine the wavelength interval of the snow spectral albedo simulation to 1 µm. The narrowband snow albedo targeting MOD09GA reflectance data simulated by the integrated model spatially agrees well with the narrowband snow albedo retrieved by utilizing the MOD09GA reflectance data and has a smaller deviation in the time series. The results of our research prove the feasibility of simulating the snow radiation observed by satellite remote sensing while utilizing an integrated snow hydrological model and snow radiative transfer model based on ground-driven data.
Our results suggest that if real-time meteorological forcing data and snow optical data can be obtained, the new integrated model will be able to predict spatiotemporally continuous snow spectral albedo and further predict narrowband snow albedo by targeting different optical remote sensing satellite observations. We developed a new method representing an important attempt to simulate snow reflectance information directly from land surface models and ground-driven data targeting satellite remote sensing observations, effectively eliminating the impacts of multiple intermediate transformations on the accuracy of directly simulated remote sensing observations. The results of this study will be beneficial for the application of multisource remote sensing data to snow radiation monitoring in land surface integrated models. Moreover, the new land surface integrated model can compensate for the shortcomings of satellite remote sensing observations in time series.
Author Contributions: D.S. and H.L. designed the study and wrote the paper; W.X., J.W. and X.H. contributed to the discussions, edits, and revisions. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
No potential conflict of interest was reported by the authors.

Appendix A. Variation Characteristics of Snow Spectrum Albedo
As indirect evidence for the spatial distribution of snow spectral albedo, we use the control experiment method to simulate the change characteristics of the snow spectral signature under a specific solar zenith angle and effective snow grain size based on the new integrated model. region and refine the wavelength interval of the snow spectral albedo simulation to 1 μm. The narrowband snow albedo targeting MOD09GA reflectance data simulated by the integrated model spatially agrees well with the narrowband snow albedo retrieved by utilizing the MOD09GA reflectance data and has a smaller deviation in the time series. The results of our research prove the feasibility of simulating the snow radiation observed by satellite remote sensing while utilizing an integrated snow hydrological model and snow radiative transfer model based on ground-driven data.
Our results suggest that if real-time meteorological forcing data and snow optical data can be obtained, the new integrated model will be able to predict spatiotemporally continuous snow spectral albedo and further predict narrowband snow albedo by targeting different optical remote sensing satellite observations. We developed a new method representing an important attempt to simulate snow reflectance information directly from land surface models and ground-driven data targeting satellite remote sensing observations, effectively eliminating the impacts of multiple intermediate transformations on the accuracy of directly simulated remote sensing observations. The results of this study will be beneficial for the application of multisource remote sensing data to snow radiation monitoring in land surface integrated models. Moreover, the new land surface integrated model can compensate for the shortcomings of satellite remote sensing observations in time series.
Author Contributions: D.S. and H.L. designed the study and wrote the paper; W.X., J.W. and X.H. contributed to the discussions, edits, and revisions. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
No potential conflict of interest was reported by the authors.

Appendix A. Variation Characteristics of Snow Spectrum Albedo
As indirect evidence for the spatial distribution of snow spectral albedo, we use the control experiment method to simulate the change characteristics of the snow spectral signature under a specific solar zenith angle and effective snow grain size based on the new integrated model. The snow spectral albedo variation curve simulated based on the integrated model is basically consistent with that measured by a SVC (Spectra Vista Company) HR-1024 ground-based spectral The snow spectral albedo variation curve simulated based on the integrated model is basically consistent with that measured by a SVC (Spectra Vista Company) HR-1024 ground-based spectral radiometer in our practical observations [74]. Snow albedo is very sensitive to changes in the effective snow grain size and solar zenith angle, and the most sensitive bands are concentrated mainly in the VIS and NIR regions. Snow spectral albedo decreases with increasing effective snow grain size and increases with increasing solar zenith angle ( Figure A1). Our results show that the sharply decreasing stage of snow spectral albedo is concentrated mainly in the NIR band. Under parameters with different sensitivities, the change amplitude of snow albedo in the NIR band is also larger than that in the VIS band and is concentrated predominantly in the wavelength range of 0.7-1.4 µm. The main reason for this finding is that the variations in snow diffuse parameters in the NIR band with the snow grain size are more sensitive than those in the VIS band [75].

Appendix B. Variation in Snow Spectral Albedo with the Wavelength Simulated by CRREL
Remote Sens. 2020, 12, x FOR PEER REVIEW 26 of 31 radiometer in our practical observations [74]. Snow albedo is very sensitive to changes in the effective snow grain size and solar zenith angle, and the most sensitive bands are concentrated mainly in the VIS and NIR regions. Snow spectral albedo decreases with increasing effective snow grain size and increases with increasing solar zenith angle ( Figure A1). Our results show that the sharply decreasing stage of snow spectral albedo is concentrated mainly in the NIR band. Under parameters with different sensitivities, the change amplitude of snow albedo in the NIR band is also larger than that in the VIS band and is concentrated predominantly in the wavelength range of 0.7-1.4 μm. The main reason for this finding is that the variations in snow diffuse parameters in the NIR band with the snow grain size are more sensitive than those in the VIS band [75]. Figure A2. Typical spectral albedo curve for snow [60].

Appendix B. Variation in Snow Spectral Albedo with the Wavelength Simulated by CRREL
The variation in snow spectral albedo measured by CRREL is shown in Figure B1. The snow spectral albedo is relatively high between 0.6 μm and 0.7 μm, and it continues to decrease in the NIR band until a trough occurs at approximately 1.0 μm. A small peak appears at approximately 1.09-1.10 μm. Beyond 1.1-1.5 μm, the snow spectral albedo decreases rapidly, and between 1.25 μm and 1.35 μm, the snow spectral albedo decreases gradually. Other crests appear again at 1.83 μm and 2.24 μm. There is strong attenuation within 1.95-2.05 μm, and the reflection peaks at subsequent wavelengths fluctuate less.  Figure A2. Typical spectral albedo curve for snow [60].

Appendix C. Descriptions of Snow Albedo Models
The variation in snow spectral albedo measured by CRREL is shown in Figure A2. The snow spectral albedo is relatively high between 0.6 µm and 0.7 µm, and it continues to decrease in the NIR band until a trough occurs at approximately 1.0 µm. A small peak appears at approximately 1.09-1.10 µm. Beyond 1.1-1.5 µm, the snow spectral albedo decreases rapidly, and between 1.25 µm and 1.35 µm, the snow spectral albedo decreases gradually. Other crests appear again at 1.83 µm and 2.24 µm. There is strong attenuation within 1.95-2.05 µm, and the reflection peaks at subsequent wavelengths fluctuate less.