Improved Parameterization of Snow Albedo in WRF+Noah:Methodology Based on a Severe Snow Event on the Tibetan Plateau※

Lian LIU,Massimo MENENTI,Yaoming MA*,4,5,6,7,8,and Weiqiang MA,6

1Land-Atmosphere Interaction and its Climatic Effects Group, State Key Laboratory of Tibetan Plateau Earth System,Resources and Environment (TPESRE),Institute of Tibetan Plateau Research,Chinese Academy of Sciences,Beijing 100101,China

2State Key Laboratory of Remote Sensing Science,Aerospace Information Research Institute,Chinese Academy of Sciences,Beijing 100094,China

3Delft University of Technology,Delft 2600AA,Netherlands

4College of Earth and Planetary Sciences,University of Chinese Academy of Sciences,Beijing 100049,China

5College of Atmospheric Science,Lanzhou University,Lanzhou 730000, China

6National Observation and Research Station for Qomolongma Special Atmospheric Processes and Environmental Changes,Dingri 858200,China

7Kathmandu Center of Research and Education,Chinese Academy of Sciences,Beijing 100101,China

8China-Pakistan Joint Research Center on Earth Sciences, Chinese Academy of Sciences, Islamabad 45320, Pakistan

ABSTRACT Snowfall and the subsequent evolution of the snowpack have a large effect on the surface energy balance and water cycle of the Tibetan Plateau (TP).The effects of snow cover can be represented by the WRF coupled with a land surface scheme.The widely used Noah scheme is computationally efficient,but its poor representation of albedo needs considerable improvement.In this study,an improved albedo scheme is developed using a satellite-retrieved albedo that takes snow depth and age into account.Numerical experiments were then conducted to simulate a severe snow event in March 2017.The performance of the coupled WRF/Noah model,which implemented the improved albedo scheme,is compared against the model’s performance using the default Noah albedo scheme and against the coupled WRF/CLM that applied CLM albedo scheme.When the improved albedo scheme is implemented,the albedo overestimation in the southeastern TP is reduced,reducing the RMSE of the air temperature by 0.7°C.The improved albedo scheme also attains the highest correlation between the satellite-derived and the model-estimated albedo,which provides for a realistic representation of both the snow water equivalent (SWE) spatial distribution in the heavy snowbelt (SWE ＞6 mm) and the maximum SWE in the eastern TP.The underestimated albedo in the coupled WRF/CLM leads to underestimating the regional maximum SWE and a consequent failure to estimate SWE in the heavy snowbelt accurately.Our study demonstrates the feasibility of improving the Noah albedo scheme and provides a theoretical reference for researchers aiming to improve albedo schemes further.

Key words:WRF,MODIS,severe snowfall,albedo scheme,SWE,Tibetan Plateau

1.Introduction

The Tibetan Plateau (TP),also known as the“Third Pole”(Qiu,2008),constitutes a unique geographical environment with an average elevation of more than 4000 m above sea level and an area of approximately 2.57×106km2(Zhang et al.,2002).Snow frequently falls over the TP and is a vital cryosphere component,linking cryospheric and hydrospheric processes with important regional weather and climate feedback.Snow on the TP is strongly linked to the general atmospheric circulation and can affect the Asian summer monsoon (Zhang et al.,2004;Seol and Hong,2009;Xiao and Duan,2016;Li et al.,2018).

Snowfall and the subsequent evolution of the snowpack influence both the energy budget and the water cycle(Barnett et al.,1988) due to high variability in the surface net radiation flux,snow water storage,and snow meltwater release (Liu and Qian,2005).As a key factor of net radiation flux,albedo directly determines the absorption of solar irradiance at the surface and is,therefore,an important parameter for land surface processes (Sellers et al.,1996).Subtle changes in albedo affect the energy balance for the earth-atmosphere system and amplify climate variability(Bloch,1964).Albedo is affected by snow density,water content,snow particle size,and impurity content (Green et al.,2002;Painter et al.,2003;Hansen and Nazarenko,2004),and is also related to ice and snow contamination or debris deposition,temperature,and terrain,slope,and aspect (Jonsell et al.,2003).Atmospheric water content,turbidity,and cloud cover change the amount and spectral properties of surface irradiance,thus affecting surface albedo.

Many parameterization schemes of snow albedo have been developed and improved upon.Snow particle size,age,depth,density and melt rate,and air temperature have been considered (Brock et al.,2000;Greuell,2000;Klok et al.,2003;Malik et al.,2014).Solar zenith angle and cloud cover have been accounted for (Dickinson et al.,1986),and snow particle size and depth,and impurity content have been considered along with solar zenith angle in other models (Marshall and Warren,1987;Marshall et al.,1999;Marshall et al.,2003).Gardner and Sharp (2010) developed a new albedo parameterization scheme that considers the carbon concentration on the snow surface,solar zenith angle,cloud optical thickness,snow depth (SD),and snow and ice cover.Their parameterization provides an accurate representation of albedo over a wide range of snow and ice conditions.Yang et al.(2016) evaluated four different albedo parameterization schemes in various climate models.They proceeded to develop a modified scheme that considered synoptic conditions and the local land and sea-ice surface characteristics.Zhong et al.(2017) improved a snow albedo scheme in the Snow-Atmosphere-Soil Transfer model by considering the impacts of light-absorbing aerosols on snow,and better-reproduced snow albedo and SD,particularly during the period of snow ablation.

Boundary layer conditions in a coupled land-atmosphere model depend on the applied land surface schemes where albedo is an extremely significant factor.Some shortcomings in calculating land surface albedo in current land surface models (LSMs) exist.For example,the community Noah LSM employs four soil layers and a single snow layer to simulate snow accumulation,melting,sublimation,and energy exchange at snow-atmosphere and snow-soil interfaces (Chen and Dudhia 2001;Ek et al.,2003).This LSM appears to be the most readily available snow albedo scheme for long-term climate modeling research,despite its various predetermined constant parameters and simplified treatment of snow (Rai et al.,2019).However,the albedo scheme in Noah LSM only considers snow cover and age while ignoring snow thickness,known to be the most important factor influencing the albedo variation (Yang et al.,2016).This omission leads to inaccuracies in the estimated albedo during snowfall and the subsequent snowmelt processes (Liu et al.,2019).In the Community Land Model(CLM),a sophisticated albedo scheme considers multiple snow properties,e.g.,SD,snow cover,snow age,snow metamorphism,fresh snow content,liquid water content,effective ice grain size,and impurity content (Dai et al.,2003).Consequently,the advanced albedo scheme in CLM outperforms the Noah albedo scheme in modeling albedo and snow-related variables (Liu et al.,2019).Its inclusion enables the Weather Research and Forecasting (WRF) air temperature estimates to be more accurate when applying CLM than the Noah LSM (Jin et al.,2010).Nevertheless,fast computation is a major requirement for numerical land-atmosphere models with high spatial and temporal resolutions,and the CLM is much more computationally expensive than the Noah LSM.Therefore we chose to improve upon the albedo scheme in the WRF coupled with the Noah LSM to improve the model’s performance for snow event simulations.

In our study,the albedo was parameterized as a function of SD,snow age,fresh snow albedo,and snow-free albedo by adapting the scheme in Oerlemans and Knap(1998) using the observed snow-free albedo.The application of the Noah LSM to the WRF was confirmed to be able to verify not only the observed values but also the spatial pattern of SD at the meteorological stations (Liu et al.,2019).Therefore,apart from the Moderate Resolution Imaging Spectroradiometer (MODIS) albedo products and ground SD measurements,our improved albedo scheme also applies WRF estimates of SD.We repeated our earlier (Liu et al.,2019)numerical experiments for the March 2017 snow event over the whole TP to compare the performance of three model configurations:(1) WRF+Noah LSM applying our improved albedo scheme,(2) WRF+Noah LSM,and (3) WRF +CLM.The default albedo schemes are used in (2) and (3).

The remainder of this paper is organized as follows.Section 2 describes the data,numerical experiments design and methods of improving albedo scheme,section 3 provides the results about near-surface air temperature,albedo and snow,section 4 discusses the performance of our improved albedo scheme and the influence of land surface parameters on the boundary meteorological variables estimates and section 5 concludes our findings about the improved albedo scheme and its capacity on simulating the boundary meteorological variables.

2.Data and Methodology

2.1.In situ observations and MODIS products

A severe snow event took place over a large fraction of the TP in March 2017 with a maximum SD deeper than 70 cm in Nyalam County of the Tibetan Autonomous Region,China.This extreme event featured ground observations of SD never recorded before in the eastern TP (i.e.,in Minyang),and within the thickest 4thpercentile in the eastern and southern TP -Nyalam and Songpan -by historical statistics of SD observations for more than 50 years.Variable snowfall intensity was observed,from light snow in some regions to heavy snow in the eastern,central,and southern edges of the TP.The sensitivity of this snow event to the applied LSM and initial and boundary conditions was investigated in our previous study (Liu et al.,2019).

This study focuses on the same snow event but includes an improved albedo parameterization and an updated set of static land surface parameters in WRF+Noah LSM,such as land cover and fractional vegetation cover (FVC).Observations of near-surface air temperature and snow water equivalent (SWE) were taken from 32 World Meteorological Organization surface synoptic observation (SYNOP) stations run by the China Meteorological Administration (CMA),and albedo observations were taken from six observatories operated by the Chinese Academy of Sciences (CAS).The geographic information and distribution of in situ stations are shown in Table 1 and Fig.1.

Table 1.Information of in situ China Meteorological Administration (CMA) and Chinese Academy of Sciences (CAS) stations.

Due to the sparse and uneven distribution of the few CAS observatories,the in situ observations of albedo are poorly representative of the TP.Spaceborne observations of radiance reflected by the Earth-Atmosphere system provides an effective means to monitor land surface albedo both regionally and globally (Li and Garand,1994),and MODIS albedo products are more accurate than other satellite data products when evaluated against in situ observations (Liang et al.,2002,2005).The products MOD09A1,MOD09CMG,and MYD09CMG were used to calculate surface albedo by converting narrowband reflectance to broadband albedo following Liang (2000).The product MOD09A1 is a MODIS/Terra eight-day surface reflectance product with a spatial resolution of 500 m,provided on a sinusoidal coordinate grid.The product MOD09CMG(MYD09CMG) is a MODIS/Terra (MODIS/Aqua) daily surface reflectance product with a spatial resolution of 0.05°,provided on the Climate Modeling Grid (CMG).The three surface reflectance products are L3 data products and provide surface reflectance in bands 1 to 7,from which broadband albedo can be estimated (Liang,2000):

where αshortis the surface broadband albedo; α1-α7is the surface reflectance in MODIS bands 1 to 7.The spectral coverage for MODIS bands 1 to 7 is 0.62-0.67,0.84-0.87,0.46-0.48,0.54-0.56,1.23-1.25,1.63-1.65 and 2.11-2.15 μm respectively.

The dataset MCD12Q1 is the MODIS Terra+Aqua combined land cover product.It is a yearly L3 level product with a spatial resolution of 500 m,provided on a sinusoidal coordinate grid.It contains five land cover classification schemes and is calculated using a supervised decision tree classification algorithm.The primary land cover scheme identifies 17 land cover classes defined by the International Geosphere-Biosphere Program (IGBP).Due to large differences in reflected radiation for different land cover types,the land cover 2017 product from MODIS was used to calculate land surface albedo using the IGBP classification in this study.

2.2.Noah LSM default and improved albedo parameterization schemes

The default surface albedo parameterization scheme in the Noah LSM considers snow cover and age (Ek et al.,2003;Livneh et al.,2010).The Noah LSM default albedo scheme is described by:

whereLvis a prescribed empirical parameter; αs0is a prescribed maximum albedo for each land cover type;tis the number of days since the last snowfall;AandBare constants set to 0.94 and 0.58,respectively for the snow accumulation period,i.e.,from October to February,and set to 0.82 and 0.46 otherwise; αbgis the background albedo,andfcis the fractional snow cover.Snow cover is set equal to 1.0 when SWE meets or exceeds the prescribed threshold SWE.Otherwise,snow cover is parameterized using the following equations:

whereSupis a prescribed threshold SWE in meters that implies 100 percent snow cover for each land cover type;Seqvis SWE in meters;andP0is a constant shape parameter of the distribution function of snow cover (set to 2.6 in the Noah default scheme).

The Noah LSM default albedo scheme (Livneh et al.,2010) imposes a seasonally-variable decay in albedo from its initial fresh snow value.The decay is slower (faster) during the accumulation (ablation) season [Eq.(3)],with an initial maximum snow albedo value that is usually very high[Eq.(2)].Apart from snow surface albedo,this scheme also includes the influence of ground surface albedo through the fractional snow cover [Eq.(4)],which leads to some shortcomings.For example,snow cover in the Noah LSM is calculated using an areal depletion threshold (Sup),which assumes partial snow cover when SWE is below the threshold [Eqs.(5) and (6)].Otherwise,snow cover is set equal to 1.0.The parameterSupfor non-forest land cover is fixed at 0.02 m.In summary,the areal snow cover fraction remains at 1.0,and the albedo is only determined by snow age while the snow is thicker than 0.1 m (Ek et al.,2003;Livneh et al.,2010).This approach is independent of the influence of ground albedo.Therefore,regarding the greater sensitivity to SD than snow cover for the snow accumulation and ablation periods,considering SD explicitly in the Noah LSM albedo scheme may be a better choice.

Fig.1.Land cover types and in situ stations over the Tibetan Plateau and the surrounding regions;the legend number indicates the land cover type index,with the associated specific description to the right of each index;solid red circles represent the location of in situ stations with the station number near its position.For a description of station numbers,see Table 1.

To improve the Noah LSM albedo scheme,we used the snow albedo parameterization proposed by Oerlemans and Knap (1998),which takes into account SD and snow age explicitly,and includes an application specific to glaciers:

wheresis the last day with snowfall prior to dayi,ands-i,thus giving the snow age in days.The parameters αfirn,αfreshsnow,and αiceare the albedos of firn,fresh snow,and bare ice,respectively;dis SD in meters;t*andd*are scaling parameters for the snow age and SD,respectively; α(i) is actual albedo on dayi.Equation (9) determines whether or not snowfall has occurred.

Land cover in the TP is fragmented (Fig.1) and the snow-free surface albedo depends on land cover and location.Our first change to the albedo scheme was to replace αicein Eq.(8) with the snow-free albedo appropriate to the land cover at a given location.This change accounts for the large differences in reflected radiance between different land cover types (Dong and Li,1994;Hu et al.,2019).For example,the reflectance of a water body is much lower than the reflectance of bare soil or a sparsely vegetated surface(Menenti et al.,1989).

2.3.Model Configuration and Experiment Design

The non-hydrostatic WRF model (Skamarock et al.,2008),version 3.7.1 (released August 2015),was used for this study with a horizontal resolution of 25 km.The WRF experiments were configured for a single large domain,with an upper-right boundary at 46°N and 110°E,and a lowerleft boundary at 20°N and 60°E,to fully include regions important to the Indian monsoon and westerlies.We ran the model for 10 days and 18 hours starting at 0800 Beijing Standard Time (LST) on 5 March 2017,which was at least one day before snowfall.This temporal gap allowed the model ample time to stabilize and rapidly converge against a physically dynamic balanced state.The model initial and boundary conditions were provided by the European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis dataset (ERA-Interim) provided at a 0.25° spatial resolution and six-hour temporal resolution.The model was configured to use the Noah LSM to describe all land-atmosphere interactions;the Lin scheme to represent microphysical processes;the RRTM scheme to describe longwave radiation;the Dudhia scheme to represent shortwave radiation;the YSU scheme to describe the planetary boundary layer,and the Kain-Fritsch cumulus parameterization scheme for convective clouds.

Four experiments were carried out with WRF,each implementing a different albedo parameterization scheme (Table 2).The default Noah LSM albedo scheme [Eqs.(2)-(4)] was implemented in the control experiment (EXP1),with hourly model output.In the second experiment (EXP2),the improved albedo scheme [Eqs.(7)-(9)] was implemented in WRF+Noah LSM,using remote sensing products(MOD09A1,MOD09CMG,and MYD09CMG) and sparse in situ observations of SD to accurately estimate αfirn,αfreshsnow,t*,d*,and bare ground albedo in the albedo parameterization.Although the WRF model simulation has model errors including systematic deviation and mode integration error accumulation,the WRF+Noah LSM appears to be able to verify not only the spatial pattern but also theobserved values of SD at the sparse meteorological stations,confirming the model’s performance of regional snow estimates (Liu et al.,2019).Therefore,to improve WRF model albedo estimates during the snow event,the modeled SD was considered to generate the improved albedo scheme,which was implemented in the third experiment (EXP3).The third experiment (EXP3) was similar to EXP2 but used hourly estimates of SD from EXP1.Model results for EXP2 and EXP3 were output at six-hour intervals.To evaluate the performance of the WRF+Noah LSM using the improved albedo scheme,i.e.,EXP2 and EXP3,we also considered results from our previous numerical experiments with WRF +CLM,in which the more complex WRF+CLM albedo parameterization was used (Liu et al.,2019),referred to here as EXP4.The first day was used for model spin-up.

Table 2.Overview of the design of numerical experiments with the WRF.

Land cover and vegetation are important factors in the snow event simulation through their influences on the estimation of surface albedo and the interception of snowfall.In this study,the four experiments used the same default static land surface products to define land cover and green vegetation coverage.The default land cover dataset in the model was produced by the United States Geological Survey(USGS),using multispectral image data acquired by the Advanced Very High Resolution Radiometer (AVHRR)from April 1992 to March 1993,and adopted the 24 classification categories from the USGS.The default green vegetation coverage dataset in the model has a spatial resolution 0.144°,and was derived by Gutman and Ignatov (1998)using the AVHRR Normalized Difference Vegetation Index(NDVI) from 1985 to 1990.

2.4.Estimation of parameters for the improved albedo scheme

The albedo parameterization scheme is expanded to represent albedo as a function of SD,snow age,snow-free albedo,fresh snow albedo,firn albedo [Eqs.(7)-(9)].Snow depth (SD) and snow age were taken from ground observations and WRF outputs.Fresh snow albedo and snow-free albedo were taken from data retrieved from MODIS.Firn albedo and the scales for snow age and SD were estimated from nonlinear fitting.Based on the WRF+Noah LSM,the expanded albedo scheme was used in EXP2 and EXP3.

2.4.1.Estimation of snow free and fresh snow albedo

The MOD09A1 and MCD12Q1 products from 2017 were used to determine snow-free albedo and fresh snow albedo.More specifically,these data were used:

(1) Broadband albedo was calculated from the MOD09A1 dataset and Eq.(1) for each pixel.

(2) Pixels with MODIS retrieved albedo greater than 0.7 and lower than 0.3 were classed as fresh snow and snow free,respectively.Thus distributions of fresh snow and snow-free albedo were assigned to each specific land cover type,assuming land cover types as mapped in the MCD12Q1 product.

(3) The third quartile of the distribution of MODIS retrieved fresh snow albedo,selected as described in (2),was used as fresh snow albedo [αfreshsnowin Eq.(7)] for each specific land cover type.Similarly,the third quartile of the distributions of snow-free albedo found for each land cover type in (b) was used for land cover specific αsnowfree,which we replaced with αicein Eq.(8).

Taking all land cover types together,the first and the third quartiles of the fresh snow albedo are 0.71 and 0.84,respectively (Fig.2),and the first and third quartiles of the snow-free albedo are 0.02 and 0.27,respectively.For simplicity,the third quartile values for the two subsets of different land cover types were averaged to give the fresh snow albedo (0.79) and snow-free albedo (0.19),which were then used in the improved albedo parameterization when it was implemented in the Noah LSM scheme.

2.4.2.Dependence of snow albedo on snow age

Fig.2.Box plot of broadband albedo higher than 0.7 and lower than 0.3 in various land cover indexes (LU index),the first and third quartile values are marked beside each box.

Surface spectral reflectance from the MOD09CMG product and hourly SD estimates from EXP1 were used to calculate the parameters needed for the improved albedo scheme [Eqs.(7) and (8)] in EXP3.The in-band reflectances from the MODIS product were combined to estimate the broadband albedo using Eq.(1).We combined the WRF SD values with the albedo values calculated from the retrievals corresponding to the time closest to the snowfall.The snowfall time was extracted from the hourly WRF output for EXP1.As SD increases,the albedo increases for shallow snow but remains constant for deep snow (where the albedo is already high).We define deep snow as snow where SD is at least 20 cm.Therefore,snow albedo is considered equal to the albedo retrieved from the MODIS data,where the SD is at least 20 cm.The MODIS albedo for locations where the SD values from WRF (EXP1) are at least 20 cm was used to determine and,for use in the improved albedo scheme in EXP3 [Eq.(7)].

To determine the relationship between the ground observations of SD and albedo,MODIS surface spectral reflectance products from the Terra (MOD09CMG) and Aqua(MYD09CMG) platforms were used in the improved albedo scheme for EXP2.Similar to EXP3,where ground observations of SD are at least 20 cm,the albedo from MODIS was used for αfirnandt*in the improved albedo scheme [Eq.(7)]for EXP2.

2.4.3.Dependence of albedo on SD

Spatial variability in SD and fragmentation of snow cover immediately after snowfall are both high for this snow event.The modeled SD and MODIS retrieved albedo during the modeled snowfall are shown in Fig.3.It is clear that albedo is less than 0.7 for a large range of SD estimates (values from WRF),with a wide range of values.For example,albedo corresponding to a SD greater than 100 cm varies between 0.1 and 0.75 at higher elevations in the Himalayas,which illustrates that the modeled snowfall is actually a snowmelt process in this complex terrain region(Fig.3).In these conditions,the difference may be attributable to the difference in spatial resolution between the MODIS (500 m×500 m) and WRF (25 km×25 km) data,from which the albedo and SD were taken.The different resolutions mean that the highly variable surface types and elevation are captured to different degrees in the two datasets.It is concluded that the coarse grid resolution WRF generally overestimates SD for the high elevations in the Himalayas.We,therefore,used SD values from WRF only where the estimates were less than 100 cm to determined*for use in the improved albedo scheme [Eq.(8)] in EXP3.We considered snow albedo equal to the fresh snow albedo during the snowfall event.The ground observed SD values were used to estimated*in the improved albedo scheme [Eq.(8)]in EXP2 during snowfall.

The steps to estimate the albedo-dependent parameters(i.e.,αfreshsnow,αsnowfree,αfirn,t*,andd*) in the improved albedo parameterization scheme are shown in Fig.4.The values used for these parameters in EXP2 and EXP3 are listed in Table 3.Our value for firn albedo in EXP2 is 0.51,which is similar to values from other studies,e.g.,0.53 in Oerlemans and Knap (1998) and 0.5 in Yang et al.(2013),but is different than that used in EXP3 (Table 3).

Table 3.Optimal values for parameters in the improved albedo parameterization scheme.

2.5.Evaluation of the Model’s Performance

To assess the performance of the WRF when applying various albedo parameterization schemes,we compared field observations of near-surface air temperatures and SWE from 32 CMA stations,albedo values from 6 CAS stations,and MODIS retrieval products with the model’s estimated values.At local solar noon in Lhasa (1400 LST,LST=UTC+8 h),the observed albedo value is closer to the Lambertian albedo described by the WRF model when coupled with LSMs.Thus,we used albedo observations at 1400 LST to evaluate the model estimated albedo.The root mean square error (RMSE),mean absolute deviation (MAD),temporal correlation coefficient (CC),spatial correlation coefficient (SCC),and relative difference (δ) were used to evaluate the model’s performance [Eqs.(10)-(13)].The calculation of the SCC follows that of the CC,but areas are weighted by the sine of the latitude.Weighting is also performed appropriately for unequally spaced grids.The SCC was calculated in the region from 65°E to 106.8°E,and from 25°N to 40.8°N,for 10752 grids in total.Specifically:

Fig.3.(a) Scatterplot of snow depth (SD) from EXP1 estimates and Terra MODIS albedo and (b)spatial distribution of SD estimates from EXP1 during snowfall.

whereNis the total number of observed or simulated data(N=10 when evaluating the performance of albedo estimates at one station,andN=40 when evaluating the performance of near-surface air temperature estimates at one station);are the modeled and observed values at timestepj,respectively;are the mean of the modeled and observed values,respectively; δ(j) is the relative difference at timestepjin units of percentage.

3.Results

3.1.Near-surface air temperature

Fig.4.Flowchart illustrating the steps to estimate the parameters for the improved albedo parameterization scheme in EXP2 and EXP3 using MODIS albedo products and snow depth (SD).

We assessed the air temperature calculated by the WRF using the different albedo parameterization schemes by comparing it against CMA observations (Fig.5).Compared to the ground observations,the air temperature in EXP1 has a high RMSE.The highest RMSE for air temperature is from EXP2,0.6°C higher than the RMSE for EXP1,and the CC between the model calculated air temperature and the ground observations is also lower for EXP2 than for EXP1(CC=0.66 for EXP2).The RMSE for EXP3 is 0.7°C lower than for EXP1,with a CC value of 0.71.The lowest RMSE corresponds to EXP4,0.2°C lower than for EXP3,which also has the highest CC (0.75).The MAD for all experiments (Fig.5) shows similar results.For EXP1,the MAD is 8.2°C,while the MAD is 0.4°C higher for EXP2 than for EXP1.The MAD for EXP3 and EXP4 is lower than EXP1 by 0.7°C and 0.9°C,respectively.

Fig.5.(a) RMSE (units:°C),temporal correlation coefficient (CC) and (b) mean absolute deviation (MAD,units:°C) for near-surface air temperature,comparing model estimates from EXP1,EXP2,EXP3,EXP4,and ground observations.

Albedo is the main determining factor for the net radiation flux;thus,an overestimated albedo leads to an underestimated net radiation flux,which may contribute to the large cold bias in modeled near-surface air temperatures for the TP.The relatively simple representation of albedo in Noah LSM leads to a strong cold air temperature bias of -8°C in EXP1,while the improved albedo scheme in the Noah LSM in EXP3 and the CLM (applying the very complex and advanced albedo scheme in EXP4) reduces the cold air temperature bias by 1°C.These results illustrate that the bestmodeled air temperature estimates may correspond to realistic albedo estimates,in agreement with Meng et al.(2018),who reduced a cold near-surface air temperature bias by 1.8°C in the WRF by updating the lower boundary conditions using the MODIS time-varying albedo.The improved albedo parameterization based on the Noah LSM and MODIS albedo products was implemented in EXP2 and EXP3.In EXP2,in situ observations of SD were used,and the hourly SD calculated by EXP1 was used in EXP3.The accuracy of the near-surface air temperature calculations is weakest in EXP2 because of the sparse and uneven distribution of in situ stations,which are especially sparse in the northwestern TP,and also because the spatial locations of the ground observations do not coincide with the WRF grid locations for the SD calculations.The spatial matching between SD estimates in EXP1 and the albedo calculations in EXP3 results in reasonably accurate model performance.However,the method assumes that the SD calculated in EXP1 is accurate.This improved performance in EXP3,relative to EXP1,is evident by the much lower RMSE and MAD;additionally,the cold bias in EXP3 is greatly reduced,indicating a comparable performance with the CLM.

3.2.Albedo

The albedo estimated at 1400 LST in our four numerical experiments was evaluated against the MODIS albedo product by calculating the SCC (Fig.6) and the relative difference (Fig.7) between them.The MODIS albedo data (Fig.6a)shows that albedo exceeded 0.5 for a large fraction of the TP,i.e.,most of the TP is covered with snow with a snowbelt evident in the southern TP.Fresh snow with higher albedo coexists with melting snow in the eastern and northwestern TP.Old snow with lower albedo is melting along the northern boundary of the TP,with sparse snow cover in the southeastern TP.The difference in albedo between the EXP1,EXP2,and EXP3 estimates and the MODIS retrievals show these experiments (EXP1,EXP2,and EXP3) overestimate the relatively high (low) albedo in the northwestern(southeastern) TP,with local relative differences exceeding 80%.At the same time,they underestimate the relatively high albedo in the snowbelt in the southern TP (Figs.6b-d,Figs.7a-c).Experiments EXP1 and EXP2 overestimate albedo at the northern boundary of the TP,but this overestimation is greatly reduced in EXP3.The albedo for a 250 km long belt centered on Namco Lake is overestimated in EXP2,with the mean relative difference being less than 40%.The albedo in EXP4 is underestimated over nearly the whole TP by 60% but is overestimated by more than 80% in small regions in the southeastern TP and the southern Himalayas (Figs.6e and 7d).

During the severe snow event,the SCC between the EXP1 albedo estimates and the MODIS data is between 0.31 and 0.46,significantly higher than EXP4.During the continued snowfall after 11 March,EXP2 and EXP3 featured a higher SCC than EXP1.The maximum SCC values occur on 13 March and are 0.53 and 0.5 for EXP2 and EXP3,respectively.On the same day,the SCC for EXP1 is 0.46,and for EXP4,the SCC is less than 0.3.The SCC for EXP1,EXP2,and EXP3 increases slightly through the snow event,while it decreases for EXP4 (Fig.6f).Thus,WRF simulations using the improved albedo scheme lead to improved estimates of albedo over the whole TP.In contrast,the configuration of WRF+CLM significantly underestimates the albedo evidenced by an SCC less than 0.4.

The observed albedo at 1400 LST at six CAS stations was used to evaluate the modeled albedo at specific grid points (Fig.8).The model performs differently at CAS stations characterized by different land cover and terrain (Fig.8a).At NASDE and SETS,where no snowfall occurred,EXP1,EXP2,and EXP3 result in a high RMSE (＞0.5),while for EXP4,the RMSE=0.25.However,the CC for EXP4 is comparable with EXP3 and was much lower than EXP1 and EXP2 at NASDE.All experiments result in a large,negative CC at SETS.Snowfall occurred at MASWE,QOMS,SHSEX,and MAQU,where the RMSE is low for all experiments,except for MASWE.A relatively high CC is achieved for all experiments at SHSEX and MAQU.At MASWE,the RMSE for the control experiment,EXP1,is slightly lower than the other experiments.The model configurations using the improved albedo scheme (EXP2,EXP3) result in a slightly lower RMSE and a significantly higher CC,showing a positive correlation.The lowest RMSE at MASWE is observed in EXP 4 but with a large,negative CC.At QOMS,EXP3 results in the most accurate albedo,having the lowest RMSE,and is the only experiment that corresponds to a positive CC.At SHSEX,EXP2 and EXP3 result in an RMSE that is slightly lower than the control experiment,EXP1,and all three share similarly high CC values(0.76-0.83).The least accurate albedo at SHSEX is calculated using WRF+CLM (EXP4),for which the RMSE is high,and the correlation is negative.At MAQU,EXP3 performs well and gives albedo estimates comparable to the control experiment,EXP1 (with a similar RMSE and CC),while EXP2 results in a higher RMSE than in EXP1 along with a slightly higher CC.In general,EXP3 performs well,i.e.,with a lower RMSE and a higher CC than the control experiment,EXP1.The highest RMSE (0.37) and CC (0.32)are from EXP2.The lowest RMSE (0.25) is obtained for EXP4 but with a CC of nearly zero (Fig.8b).

Fig.6.(a) MODIS surface albedo product and albedo difference between numerical estimates and the MODIS product on 11 March.(b) EXP1 -MODIS,(c) EXP2 -MODIS,(d) EXP3 -MODIS,and (e) EXP4 -MODIS;(f)the spatial correlation coefficient (SCC) between model estimates and MODIS albedo during the severe snow event (X-axis denotes the day in March 2017).

Fig.7.The relative difference between the albedo calculated in the model experiments and the MODIS albedo product on 11 March for (a) EXP1;(b) EXP2;(c) EXP3;(d) EXP4.

3.3.Snow water equivalent (SWE)

The SD is a good indicator of the magnitude of snow events,but it is dependent on snow density and is easily affected by snow compaction.The SWE is the product of snow density multiplied by SD and is a more accurate indicator of snow quantity.Snow water equivalent (SWE) observations at the CMA stations on the TP (Fig.9) and in surrounding regions were used to evaluate the accuracy of model SWE estimates from 11-13 March 2017.

The snow event was captured in a large fraction of the TP,with a heavy snowbelt (SWE ＞6 mm) clearly evident and oriented in the NE -SW direction,with the most severe snowfall occurring in the southern TP on 11 March.The snowfall moved north on 12 March and was observed over a large area.On 12 March,daily SWE observations show that the intensity and extent of snowfall decreased sharply in the southern TP,remained reasonably unchanged in the central TP,and increased in the eastern TP,where the maximum SWE increased by 8 mm.On 13 March,the snowfall continued to move north,and the intensity decreased in the eastern and northeastern TP,where the observed daily SWE was less than 3 mm.

Fig.8.The RMSE and temporal correlation coefficient (CC) between model estimated albedo from EXP1,EXP2,EXP3,EXP4 and observations at each station (a) NASDE;(b) SETS;(c) MASWE;(d) QOMS;(e) SHSEX;(f) MAQU),and (g) its average.

Fig.9.Daily ground observations of SWE.

The SWE estimates from EXP1 to EXP4 (Fig.10)show that all four experiments successfully simulated the characteristics of the snowfall event,particularly the northward movement from 11-13 March.The extent of the snowfall in the eastern,central,and southern TP was estimated successfully over this period.However,all experiments overestimate the intensity and extent of the snowfall in the southeastern TP.Compared with SWE estimates from the control experiment (EXP1) (Fig.10a),EXP2 simulates a smaller extent for the snowfall on 12-13 March,although larger than that calculated for 11 March,with even heavier snowfall,SWE ＞6 mm,on 12 March in the southeastern TP(Fig.10b).Experiment three (EXP3) calculates less snowfall than EXP1,except in the northwestern TP,and calculates a higher snowfall intensity in the heavy snowbelt,a smaller extent for the snowfall in the southeastern TP,and a higher snowfall intensity in the eastern TP on 12 March(Fig.10c).A snowfall extent and SWE,similar to EXP3,is calculated for EXP4 (see section 2.3),but the snowfall extent in the southeastern TP is overestimated,and the SWE is underestimated in the heavy snowbelt area on 11 March(Fig.10d).Compared with the SWE observations,EXP3 results in the most accurate estimates of SWE and snowfall extent for 11 March,including in the heavy snowbelt,possibly due to the improved albedo scheme (Figs.9 and 10).

A deeper analysis of the four experiments requires a more detailed analysis of SWE observations from 11-13 March 2017 (Fig.11a).The maximum observed SWE was 53.3 mm in the southern TP and 48.5 mm in the eastern TP,while the SWE was only 8.8 mm in the central TP.The stations collecting SWE observations on the TP were sparse and unevenly distributed.To mitigate the impact of this spatial distribution on our analyses,we defined three subregions:East (E),with 105 stations,Middle (M),with 19 stations,and South (S),with 20 stations (Fig.11a).Heavy snowfall occurred on 11-12 March,and we evaluated our model experiments based on the maximum and mean SWE values by subregion (Figs.11b-c).

As shown in Fig.11b,the maximum observed SWE is 51 mm on 11 March in the southern region and 26.9 mm on 12 March in the eastern region.In the middle region,the maximum observed SWE was less than 7 mm.The EXP1 estimates of maximum SWE in the M and E regions were reasonably close to observations.The relative differences were 34.6% on 11 March and 12.3% on 12 March in the M region,22 % on 11 March,and 22.6% on 12 March in the E region.Differences were lower in the E region on 12 March for EXP2 (4.5%) and EXP3 (12.5%) but were much larger in the M region.Due to the improved albedo scheme,EXP3 gave the best estimate of the maximum SWE in the E region on 11 March.The maximum SWE was underestimated by only 0.16 mm,corresponding to a relative difference of less than 1%.These results provide additional evidence of the advantages of implementing the improved albedo parameterization.In the E region,EXP4 significantly underestimated the maximum SWE on 11 March but overestimated it on 12 March,while it slightly underestimated the maximum SWE(by less than 1.5 mm) in the M region.All four experiments largely underestimated the maximum SWE in the S region.One possible explanation is that the modeled snowfall area is displaced relative to the observations (Figs.9 and 10).

The mean SWE was overestimated in all four experiments in the E and M regions (Fig.11c),although EXP4 underestimated the mean SWE by～50% on 12 March in the M region.Compared with the simple albedo scheme in the control experiment,the improved albedo scheme has the potential to give higher estimates of daily snowfall amount(see the discussion section for further details) in the E and M regions.In the S region,the mean SWE was underestimated by more than 50% on 12 March in all four experiments when the observed mean SWE was less than 1 mm.Snowfall in the S region is a rare event for this time of year,but the model successfully captured the general spatial pattern.There are differences in the spatial distribution of the snowfall in comparison with observations from the 20 stations,which likely explains the difference in SWE between the model and observations.Compared with the other experiments,EXP4 gave the largest underestimation of mean SWE (by 2.2 mm),with a relative difference that exceeded 35% on 11 March.

Fig.10.Daily SWE estimates from numerical experiments (a) EXP1;(b) EXP2;(c) EXP3;(d) EXP4.

Overall,the improved albedo scheme in EXP3 calculates the spatial distribution of SWE similar to that in EXP4,where the CLM advanced albedo scheme was applied.In addition,the calculations of SWE by EXP3 in the heavy snowbelt,and the maximum SWE values for the E region,were more accurate than those from EXP4.All four experiments generally underestimated the maximum and mean SWE in the S region where the snowfall was severe.This result may be attributed to the combined effect of differences between the observed and modeled snowfall pattern and the scarcity and sparseness of the 20 stations in the large and rugged S region.

4.Discussion

4.1.The performance of the improved albedo scheme

Fig.11.(a) Observations of total SWE from 11 through 13 March 2017 in E (East),M (Middle),and S (South)regions,(b) regional maximum SWE,and (c) regional mean SWE from 11-12 March 2017 in East,Middle,and South regions.

Near-surface air temperature,which drives snow accumulation and melt,is a diagnostic in WRF that is largely controlled by the land surface temperature,which is determined by the net radiation flux.It has been shown that air temperatures from both model simulations,as well as from reanalysis datasets,are much lower than ground observed air temperatures in the TP (Frauenfeld et al.,2005;Gao et al.,2011;Wang and Zeng,2012;Ji and Kang,2013;Su et al.,2013;Chen and Frauenfeld,2014;Hua et al.,2014;Chen et al.,2017).According to Chen and Frauenfeld (2014),the significant cold bias in near-surface air temperature in the TP could be due to poor model estimates of snow cover and,therefore,poor representation of the albedo feedback.To improve the accuracy of land-atmosphere interaction estimates,Bao and Lyu (2009) added solar zenith angles to the albedo scheme,which led to a temperature increase of 1.2°C,which considerably reduced the cold bias and improved the representation of diurnal ground temperature changes.Malik et al.(2014) pointed out that the Noah LSM overestimated snow albedo during springtime and that an alternative albedo scheme improved the estimates of albedo,SD,and snowmelt rate further by considering the shape of the variogram for optically thick snowpacks.This new scheme from Malik et al.(2014) yields a 0.105 albedo RMSE when applying the default Noah albedo scheme,which is reduced to 0.088 when applying the alternative albedo scheme (a relative improvement of 16%).Park and Park (2016) improved the winter albedo estimates by including vegetation information,i.e.,leaf and stem indices in the albedo scheme,which decreased the albedo RMSE by 69%.However,in our work,the optimally improved albedo scheme reduces the RMSE for the albedo from 0.35 with the default Noah scheme to 0.33 with the improved scheme,a relative improvement of only 6%,leading to a temperature RMSE decrease of 0.7°C.Such a large albedo RMSE and a small improvement in albedo and temperature estimates in this study are closely related to the WRF model configuration’s coarse grid resolution (25 km),leading to large uncertainties compared with in situ observations on the topographically complex TP.In addition,apart from systematic model deviations,the modeled temperature bias is related to initial temperature errors and the integrated accumulation of errors from the associated physics schemes.In the WRF,the near-surface air temperature is diagnosed from ground temperature,which,in turn,is determined by the surface energy budget.This budget is influenced by atmospheric conditions and weather conditions generated from shortwave and longwave radiation,convective,microphysics,and land surface parameterization schemes.The substantial temperature bias is only partially reduced by the LSM changes.The remaining temperature bias may be associated with the model initialization,systematic deviation,and integrated error accumulation from shortwave and longwave radiation,convective,and microphysics schemes that cannot be excluded from coupled experiments.Therefore,in future studies,offline Noah LSM experiments forced by observations are needed to accurately quantify the simulation error caused by the incomplete snow albedo parameterization scheme.Nevertheless,our improved albedo scheme outperforms the Noah default albedo scheme and presents comparable performance with the CLM in terms of air temperature,albedo,and snow estimates.These results demonstrate the improved albedo scheme’s potential for modeling land-atmosphere interaction during snow events.

Computational efficiency is an important consideration in operational numerical weather prediction research.The calculation of surface parameters in Noah LSM is aided by simple parameterization schemes,enabling fast calculations and low computational cost,justifying the widespread use of the Noah LSM in current mainstream operational numerical weather prediction models such as WRF (Liu et al.,2020;Thiruvengadam et al.,2020).In our study,we were surprised to find that the computational efficiency of WRF with the default Noah LSM (EXP1) was 5.5 times that of the WRF with CLM (EXP4),despite the desirable outcome that the CLM performs much better in land-atmosphere interaction estimates.The high computational efficiency is a strong motivating factor to pursue improving the performance of WRF by applying Noah LSM.The performance of the WRF coupled with the Noah LSM is greatly improved by applying an improved albedo scheme which additionally considers snow depth and the satellite-retrieved albedo;this approach shares similar high computational efficiency with the default Noah LSM and comparable modeling accuracy to that of the CLM.

This study focuses on the methodology needed to develop an advanced albedo parameterization scheme,to improve the performance of the WRF for severe snow estimation.The development of an advanced albedo scheme cannot be separated from the input of snow-related variables,e.g.,snow age and depth,acquired from WRF’s land-atmosphere coupling during just one snow event.Therefore,the methodology is based on only one severe snow event over the TP.Our preliminary results demonstrate that the improved albedo scheme greatly outperforms the default Noah LSM scheme and show its strong application potential across the TP region.However,the performance improvements from the modified albedo scheme are not universal across the TP,and should be studied further.The intensity of snowfall and snowmelt rate over the TP widely varies depending on the complex terrain and heterogeneous underlying surfaces.These results justify the need for more coarse(25 km) and nested fine (1-5 km) resolution trials,applied to more snow events over the TP,to assess the potential of the modified scheme to downscale simulations characterized by different snowfall intensities and snowmelt processes.

4.2.Impact of real-time land and vegetation cover on the model’s performance

4.2.1.Updating land cover and fractional vegetation coverage

Land cover has significant and diverse impacts on model performance.Some land surface properties,e.g.,leaf area index,surface albedo,roughness length,and FVC,are closely related to land cover and are assigned in the WRF based on land cover type.Changes in land cover bring about changes in surface albedo,thus driving changes in the surface energy budget that can affect the general circulation at local and synoptic scales.Charney (1975) and Charney et al.(1977) proposed an important biosphere-geophysics feedback mechanism.They demonstrated that overgrazing in the sub-Saharan belt destroyed surface vegetation,leading to higher surface albedo,ultimately leading to changes in the surface energy balance and an apparent radiative heat sink.The large and rapid changes in albedo during snowfall (snowmelt cycles may have similar effects) are investigated here.Kumar et al.(2013) noted that realistic satellite-derived FVC might improve forecast model skills for air temperature and precipitation.Yin et al.(2016) used near-real-time FVC and surface albedo to improve the accuracy of the Noah LSM.Zhang et al.(2017) noted that actual FVC affected air temperature calculations and suggested that the actual land cover and FVC should be used in WRF.

The default land cover dataset in WRF version 3.7.1 adopts the 24 classification categories of the USGS.This dataset was produced using AVHRR multispectral images from April 1992 to March 1993.To facilitate the comparison between the 24 categories of the USGS and the 17 categories of the MODIS IGBP,the default land cover dataset was reclassified to 17 categories according to the land cover classification of IGBP (Fig.12a).The default FVC dataset in the WRF is produced using the AVHRR NDVI from 1985 to 1990 (Fig.12b).In recent decades,the spatial distribution of green vegetation has changed over the TP in response to climate variability.The TP goes through serious desertification in the north and northwest;in contrast,green vegetation expanded to the southeast due to the advection of warm and moist air.These trends can be seen by comparing Fig.1 and Fig.12.Therefore,the WRF default land cover and FVC dataset are out of date for this region and do not capture the current conditions.To investigate the effect of land cover and FVC on the accuracy of the snow cover extent calculated in WRF for the TP,three additional numerical experiments,EXP5,EXP6,and EXP7,were carried out using EXP3 as a starting point.An overview of parameters for all experiments is listed in Table 2.In EXP5,we updated the land cover product in WRF using MCD12Q1 in 2017 (Fig.1).In EXP6,we used the same land cover as EXP5 and updated the FVC using the high spatial resolution MODIS NDVI product (MOD13A3) in March 2017.Fractional vegetation coverage (FVC) was estimated by applying a simple radiative transfer model and a semi-empirical formula [CR algorithm,Eq.(14);Carlson and Ripley,1997].The framework of EXP7 was similar to EXP6,but the FVC was estimated using a pixel dichotomy model,assuming a linear relationship between NDVI and FVC [GI algorithm,Eq.(15);Gutman and Ignatov,1998].The two algorithms that estimate FVC are:

where σf1is the FVC estimated by CR algorithm;σf2is the FVC estimated by GI algorithm; NDVImaxand NDVIminare the regional maximum and minimum NDVI,respectively.

Land cover and FVC affect the surface albedo and,therefore,affect the surface energy budget,influencing the quantity and intensity of snow events.Large differences in the FVC between the outdated data currently implemented in the default WRF configuration and the updated CR and GI retrievals can be seen in the TP and surrounding regions(Fig.12).The mean FVC in the outdated dataset is about 10% in the E and M regions and 15% in the S region.Using the GI algorithm,the mean FVC is increased by more than 11% in the E and S regions and 7% in the M region.Using the CR algorithm,the mean FVC is decreased by about 3%-5% in the E,M,and S regions.It should be noted that the CR algorithm is considered to be more accurate than the GI algorithm for FVC retrieval (Jiang et al.,2006).Hong et al.(2009) believed that it was still controversial which FVC algorithm was most suitable for WRF.

Fig.12.(a) The distribution of default land cover types from April 1992 to March 1993 and (b) fractional vegetation coverage from 1985 to 1990 in WRF,(c) real-time fractional vegetation coverage calculated using the CR algorithm,and (d) using the GI algorithm.

4.2.2.Near-surface air temperature

Due to the significance of land use and FVC in coupled land-atmosphere interactions,significant increases in the accuracy of near-surface air temperatures calculated in the WRF were obtained when updated land use and FVC products were used (Sertel et al.,2010;De Meij and Vinuesa,2014;Schicker et al.,2016;Li et al.,2020;Yan et al.,2020).The RMSE,temporal CC,and MAD between ground observations and model estimations that applied the real-time land and vegetation cover data are shown in Fig.13.Compared with EXP3,EXP5 results in air temperature estimates of similar accuracy,and EXP6 results in slightly increased model skill,as evidenced by the lower air temperature,RMSE,MAD,and similar CC values.The results of EXP7 demonstrate a lower model skill than EXP3,which is attributed to higher RMSEs of air temperature and MAD and a lower CC value (Fig.13).In recent decades,changes in land cover have been limited.Still,changes in FVC have been large in the east and southeast of the TP where CMA stations are denser,causing differences in the accuracy of air temperature calculated by the different model configurations in EXP3,EXP6,and EXP7,and also explaining the similar performance of EXP3 and EXP5.These differences are also related to the different algorithms:the CR algorithm appears to be the most appropriate to WRF with a potential for accurate estimates of air temperature in the eastern and southeastern TP,possibly because the FVC error caused by the CR algorithm is less than that introduced by GI algorithm (Jiang et al.,2006).

4.2.3.Albedo

The experiments using the updated land surface parameters result in spatial albedo patterns with similar accuracy to those simulated in EXP3 and calculate reasonable albedo values for the heavy snowbelt.These experiments generally overestimate albedo in the snow-free case relative to the albedo retrieved from satellite data,which is very low,for example,less than 0.4 in the southeastern TP.The SCC for the model calculated and satellite-derived albedo is about 0.5 during snowfall events,slightly higher than during periods of no snowfall.

Fig.13.(a) RMSE (units:°C),temporal correlation coefficient (CC) and (b) mean absolute deviation (MAD,units:°C) for near-surface air temperature,comparing model estimates from EXP3,EXP5,EXP6,EXP7 and ground observations.

The default and updated land cover types are similar at the six CAS stations,except at NASDE,where the default land cover type is mixed shrubs and grassland,which was updated to the bare soil or sparse vegetation type.The default FVC at the six CAS stations is very different from that in the updated datasets.The performance of WRF using the real-time updated land and vegetation cover data was investigated by calculating the albedo RMSE and temporal CC between ground observations and experimental estimates (Fig.14).Compared with EXP3,the model experiments using real-time updated land and vegetation cover data (i.e.,EXP5,EXP6,and EXP7) result in no improvement in the albedo calculated at the NASDE,SETS,and MASWE stations (Fig.14a-c),where poor model performance was discussed earlier regarding our previous experiments.Also,no improvement in albedo estimates is produced at the QOMS station by using the updated land and vegetation cover data in the model,although EXP5 and EXP6 reduced the albedo RMSE (Fig.14d).Coarse model resolution,i.e.,a 25 km simulation,is a potential reason for misrepresenting the underlying heterogeneous surface conditions.The results of EXP7,which used the updated land surface type and FVC retrieved with the GI algorithm,have albedo estimates that are slightly more accurate than EXP3,EXP5,and EXP6 at the SHSEX station (Fig.14e).This means that the real-time land and vegetation cover is equally important in the WRF model to estimate albedo at the SHSEX station.The experiments EXP5,EXP6,and EXP7,reduce the RMSE for albedo by 4%-13% and increase the positive correlation between the model albedo estimates and observations at the MAQU station.The CC increases by 0.16 with a relative improvement of 36% in EXP7 and by 0.24 with a relative improvement of 53% in EXP5 and EXP6 at the MAQU station (Fig.14f).In general,the improved land surface type and FVC data are necessary input data for accurate albedo estimates during the snow event using the WRF model.The updated land surface parameters,i.e.,land cover and FVC,improve the model skill for albedo estimation,reducing the RMSE by 1%-4%.In general,the updated land cover data contributes more to these improvements than the updated FVC through CR or GI algorithm (Fig.14g).

4.2.4.Snow water equivalent (SWE)

The updated land cover and FVC substantially mitigate the overestimation of SWE in the southeastern TP seen in the experiments using the default land cover types and FVC(Figs.10c and 15).The narrow snowbelt is successfully simulated in EXP7,while a relatively wide and sparsely covered snowbelt is estimated in EXP5 and EXP6.Compared with EXP3,the experiments using the updated land cover and FVC data perform better and mitigate the overestimation of regional mean SWE in the E region.The fifth experiment,(EXP5),which applied the updated land cover data,had the most accurately modeled SWE results in the E region,reducing the overestimation of mean SWE by 35% on 11 March and by 4% on 12 March,relative to the experiments using the default land cover data (Figs.10c and 15).The model skill for estimating mean SWE in the M region increases by 4%-6% in EXP7 but decreases by 25%-43% in EXP5 and EXP6,compared to EXP3.Using the updated land cover data and FVC does not result in more accurate estimates of mean SWE in the S region during the severe snowfall event.EXP7 performs better than EXP3,EXP5,and EXP6 in estimating the maximum SWE in the E and S regions during periods of heavy snowfall.In the M region where the observed maximum SWE is 6.3 mm,EXP5 and EXP6 reduce overestimation of the maximum SWE by 36%-41% on 11 March,and EXP7 similarly reduces the overestimation by 4.1 mm(33% relative improvement) on 12 March,relative to EXP3(Figs.10c and 15).In general,the updated land cover and FVC data contribute to improvements in model performance for SWE estimates.Still,it is unclear which is the best approach and which is the most appropriate dataset to best optimize the simulation of land surface properties in the WRF (see also Hong et al.,2009).

4.3.The link between land surface parameters and estimation of boundary meteorological variables

Fig.14.The RMSE and temporal correlation coefficient (CC) between model estimated albedo from EXP3,EXP5,EXP6,EXP7,and observations at each station (a) NASDE;(b) SETS;(c) MASWE;(d) QOMS;(e) SHSEX;(f) MAQU and (g) its average.

Land surface parameters in the WRF+Noah LSM affect the model simulation of snow events since they affect the net radiation flux and the land-atmosphere interaction(Fig.16).Albedo,an important factor for calculating net radiation flux,dramatically changes during snowfall and snowmelt.Albedo is more sensitive to SD than to snow cover.Therefore,considering SD explicitly in the improved albedo scheme is a better choice for improving the model’s performance in estimating albedo and air temperature.The improved albedo scheme [Eqs.(7)-(9)] imposes a variation in albedo from fresh snow to bare ground,revealing that albedo decreases with snow age because snow meltwater absorbs a proportion of the shortwave radiation.In contrast,albedo strongly increases with SD over shallow snow but varies only slightly over deep snow.The merit of the improved albedo scheme for albedo estimates has been demonstrated in the previous analysis.

Albedo is more accurately estimated in numerical experiments where the improved albedo parameterization scheme was applied.As a significant driver of land-atmosphere interactions,the net radiation flux determines the energy and water vapor exchange at the land-atmosphere interface.The net radiation flux drives snowmelt and changes the snowmelt rate,affecting snow age and SD estimates.However,the net radiation flux also influences ground temperature,thus impacting near-surface air temperature.The near-surface air temperature is diagnosed from ground temperature and is one of the rain and snow discrimination factors in the Noah LSM.Consequently,numerical experiments where the improved albedo parameterization is implemented have the potential to provide improved estimates of albedo,nearsurface air temperature,and solid precipitation.Changes in solid precipitation rates lead to changes in SD,thus contributing to improved albedo parameterization.This represents potential feedback between albedo estimates made using the improved albedo scheme and the SD estimates.

The snowmelt rate affects snow cover,and both snow cover and depth (SD) affect surface roughness and surface emissivity,leading to changes in the effectiveness of latent heat transport through air turbulence.This interaction impacts the net radiation flux and the energy partitioning between latent and sensible heat fluxes (Bowen ratio),further affecting energy and water vapor exchanges at the landatmosphere interface,thus affecting the accuracy of near-surface air temperature and solid precipitation simulated in the model.

Fig.15.SWE estimates in EXP5 (a,d),EXP6 (b,e),and EXP7 (c,f);The left column shows SWE estimates on 11 March,and the right column shows SWE estimates on 12 March.

Updating the land cover type and green vegetation fraction leads to changes in surface characteristics such as surface roughness and emissivity,and background surface albedo (Gao et al.,2008),which affect the estimation of the land surface energy budget and the partitioning of energy into latent and sensible heat fluxes (Bounoua et al.,2002;Matsui et al.,2005;Bonan,2008).This effect further impacts the surface energy balance and water cycle.For example,empirical evidence shows that a change from pasture land to greenhouse farming leads to increased albedo,resulting in strong negative radiation forcing and cooling for the changed area throughout the year,while surrounding areas may be warming (Campra et al.,2008).Fishman et al.(1994) investigated the cooling effect of high albedo sandy surfaces and found a difference in daytime air temperature between sandy surfaces and their surroundings of 1°C-2°C,associated with an albedo difference of approximately 0.4.Nair et al.(2007) further noted that shortwave radiation fluxes are higher after clearing native vegetation for agricultural purposes.

Green vegetation is a key factor for calculating snow redistribution because it can alter the near-surface wind speed (Li et al.,2000).The characteristics and spatial distribution of green vegetation determine the evolution of the SD spatial distribution (Essery and Pomeroy,2004;Yan et al.,2019).Near real-time changes in FVC affect the evolution of radiation transfer between the vegetation and the ground,changing the estimated net radiation flux.Therefore,updating the land cover type and FVC affects the accuracy of the near-surface air temperature,albedo,and solid precipitation simulated in the model during the snow event.

5.Conclusions

Fig.16.The effect of implementing the improved albedo scheme and using the updated land cover type and fractional vegetation coverage on model estimates of near-surface air temperature,albedo,and solid precipitation during the snow event.

The albedo parameterization scheme implemented in WRF+Noah LSM was investigated by simulating a severe snow event in March 2017 in the TP.The default albedo scheme in Noah LSM considers snow cover and snow age.We developed and evaluated an alternate parameterization that considers SD,representing the improved albedo scheme presented in this study.We used the observed and modeled SD,MODIS surface reflectance,and albedo products to retrieve the parameters required for the improved albedo scheme.In our study,the performance of the WRF+Noah LSM,with the improved albedo scheme was implemented and was compared with that of WRF +CLM and that of WRF+Noah LSM,with the default albedo scheme implemented.The default land surface parameters (land cover and FVC) are out of date in the WRF,and these parameters have a significant effect on the model calculations of albedo.Additional experiments were therefore carried out,which updated the land cover and FVC data that were used,and the impact on the accuracy of the modeled air temperature,albedo,and SWE was discussed.

The spatial pattern of air temperature is successfully simulated,although air temperature is generally underestimated,potentially attributable to the simplified representation of albedo in the model.Using the WRF with the default Noah LSM albedo parameterization,the accuracy of the modeled air temperature is much lower,and the mean albedo is greatly overestimated in the E and SE regions of the TP.The complex advanced albedo parameterization in the CLM,which considers solar zenith angle,effective ice grain size,and pollutants,increases the model skill in estimating albedo,resulting in the lowest RMSE and MAD for air temperature across all experiments.The improved albedo scheme,using model estimates for SD,mitigates the overestimation of mean albedo in the E and SE region and increases the albedo CC and SCC,thereby reducing the RMSE and MAD for air temperature by 0.7°C.The scheme strongly reduces the cold bias (by 1°C) in air temperature estimates,achieving comparable accuracy to that of the WRF+CLM.The improved scheme,using modeled SD,results in a higher CC and SCC for albedo than was achieved using the WRF+CLM.The sparse and uneven distribution of ground stations and differences between the observed SD spatial patterns and those simulated by the model means that the improved albedo scheme,when implemented using the observed SD,results in inaccurate estimations for air temperature and albedo.

The WRF+Noah LSM,with the improved albedo scheme implemented,results in a higher CC and SCC between model estimates and satellite retrievals of albedo,leading to an accurate spatial distribution of model SWE in the heavy snowbelt (SWE ＞ 6 mm) and also to accurate estimates of maximum SWE in the E region.The WRF+CLM underestimates albedo for a large fraction of the TP,including the heavy snowbelt,and underestimates the regional maximum SWE but reasonably estimates regional mean SWE.The maximum SWE in the S region during periods of severe snowfall is generally underestimated in all experiments because of the sparse and uneven distribution of ground stations in the southern TP and differences between the spatial patterns in the observed and modeled snow cover.

The default land cover and FVC data in WRF are out of date,but these data significantly impact model estimates of air temperature,albedo,and SWE.Using updated land cover and FVC data improves model performance in simulating the severe snow event,reducing the albedo RMSE by 1%-4%.Using the updated land cover data has a larger impact on albedo estimates than using updated FVC data,and the choice of FVC retrieval algorithm has a large impact on the accuracy of the retrieved FVC data.The CR algorithm outperforms the GI algorithm in the WRF model on air temperature estimates.The optimum choice of an FVC retrieval algorithm to retrieve data for inclusion in simulations is unclear in terms of the accuracy of air temperature,albedo,and SWE estimates in the TP.

Acknowledgements.This research was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA20060101),the Second Tibetan Plateau Scientific Expedition and Research program (STEP) (2019QZKK0103),the National Natural Science Foundation of China (Grant Nos.91837208,91637312,41830650,and 91737205),MOST High-Level Talent Grant No.G20190161018,the Chinese Academy of Sciences President’s International Fellowship Initiative Grant No.2020VTA0001,the Key Research Program of Frontier Sciences of Chinese Academy of Sciences (QYZDJ-SSW-DQC019),and Key Research and Development Projects of the Ministry of Science and Technology (2018YFC1505701).The authors express thanks to ECMWF for sharing the atmospheric reanalysis data set (ERAInterim dataset is available from http://apps.ecmwf.int/datasets/data/interim-full-daily/),to NASA for offering MODIS reflectance,land cover,and NDVI products (https://modis.gsfc.nasa.gov/),and to staff from CMA and CAS stations for very hard work in meteorological observations and for offering the data (CMA meteorological data is available from http://data.cma.cn/en;CAS albedo observation is available from https://data.tpdc.ac.cn/en/).The authors would like to acknowledge all anonymous reviewers for reviewing this paper and providing constructive comments.