Simulated cold bias being improved by using MODIS time-varying albedo in the Tibetan Plateau in WRF model

Systematic cold biases exist in the simulation for 2 m air temperature in the Tibetan Plateau (TP) when using regional climate models and global atmospheric general circulation models. We updated the albedo in the Weather Research and Forecasting (WRF) Model lower boundary condition using the Global LAnd Surface Satellite Moderate-Resolution Imaging Spectroradiometer albedo products and demonstrated evident improvement for cold temperature biases in the TP. It is the large overestimation of albedo in winter and spring in the WRF model that resulted in the large cold temperature biases. The overestimated albedo was caused by the simulated precipitation biases and over-parameterization of snow albedo. Furthermore, light-absorbing aerosols can result in a large reduction of albedo in snow and ice cover. The results suggest the necessity of developing snow albedo parameterization using observations in the TP, where snow cover and melting are very different from other low-elevation regions, and the influence of aerosols should be considered as well. In addition to defining snow albedo, our results show an urgent call for improving precipitation simulation in the TP.


Introduction
The Tibetan Plateau (TP), located at 27 • -45 • N, 70 • -105 • E, covering an area of 2.5 × 10 6 km 2 and with an average elevation of more than 4000 m above sea level, is the largest and most supreme highland in the world (Gao et al 2015b, Xu et al 2014, Yao et al 2012, You et al 2016. With high altitude, strong solar radiation, extremely varied topography and complex landscapes, the TP has a complicated impact on the regional climate, the Asian monsoon system, inland aridification and even global climate change (Duan et al 2013, Duan et al 2012, Guo and Wang 2013, Ma et al 2014, Wu et al 2012, You et al 2013. With global warming, the TP has been illustrated as warming rapidly, with the warming rates being almost two times that shown by a global average in the past 3030a (Liu and Chen 2000). As a result, the TP is regarded as an amplifier and driver of environmental change in the adjacent regions (Gao et al 2014, Yang et al 2014, Yang et al 2010, Yi et al 2014. In addition to its global effects, climate change in the TP has gained growing attention in recent years (Duan et al 2011, Yang et al 2011. However, under harsh environmental conditions, observations of the TP are very limited, especially in the vast northwestern TP, where there is almost no observation. In fact, numerical modeling is becoming a main method to study climate change quantitatively. However, due to the insufficient presence of topography and model resolution, both regional and global climate models (GCMs) produce large cold biases, particularly in winter (see details in table 1). Recently,  Gao et al (2015a) used the Weather Research and Forecasting (WRF) Model to simulate regional climate in the TP, showing that WRF can realistically capture the trend of a 2 m air temperature change in summers, while there is no pronounced improvement in the dry seasons. As a result, they called on more reliable observation combined with remote sensing data and higher-resolution modeling to improve the simulations.
Due to the limitation of the harsh environment of the TP, a finescale observation is very hard to perform. Another better way to improve the regional climate simulation is to integrate satellite remote sensing data. Many studies have illustrated that a better description or data assimilation of land surface parameters can improve regional climate modeling, land atmosphere interaction, and even air quality modeling (Meng et al 2009, Meng et al 2014a, 2014b, Ran et al 2016, Yin et al 2016. As one of the basic forcing parameters of boundary conditions in climate models, surface albedo quantifies solar energy absorbed by the earth, balancing the energy that is transferred by the atmosphere. In the TP, albedo not only influences the regional climate, but also affects the climate change of East Asia through sensible heating (Ghatak et al 2014, Wang et al 2014, Wu et al 2012. In addition, the evaluation of albedo from the Coupled Model Intercomparison Project Phase 5 (CMIP5) GCMs indicated that the seasonal cycle of albedo was overestimated, particularly in winter in the Northern Hemisphere (Li et al 2016). Wang and Zeng (2010) assessed snow albedo in four land models using in situ data over grass and forests, and suggesting that a probable update of snow albedo using the Moderate-Resolution Imaging Spectroradiometer (MODIS) albedo product may improve the model albedo presence over forest in the NCEP global model. As a result, the accurate presence of albedo is important in regional climate modeling and also for land-atmosphere interaction studies.
In this study, albedo from the Global Land Cover Facility, i.e. the Global LAnd Surface Satellite (GLASS) MODIS albedo, was integrated into the WRF model to improve 2 m air temperature simulation in the TP. Then, the reasonfor the overestimation of albedo in the model was investigated. We also discussed the applicability and further potential updates for the MODIS albedo product, as well as possible improvement when it is used in the models.

Albedo datasets
In the WRF model, the default albedo (referred to as CTL albedo) was derived from the Advanced Very High Resolution Radiometer (AVHRR) sensor onboard the National Oceanographic and Atmospheric Administration (NOAA) satellites, and this was based on monthly averaged clear-sky, snow free surface broadband albedo data globally. The time period was between April 1985 to December 1987 and January 1989 to March 1991. It was merged to a monthly climatology data with the spatial resolution of 0.25 • . Further details can be found in Csiszar (2009). In this study, albedo for updating WRF lower boundary conditions was derived from the GLASS surface albedo product (Liang and Liu 2012, Qu et al 2013. GLASS is a series of remote sensing products produced by the College of Global Change and Earth System Sciences, Beijing Normal University, China. The GLASS surface albedo product has an eight day temporal resolution, and is available from 1981-2010. The GLASS surface albedo product in 2000-2010 was derived from MODIS data. In this work, the GLASS surface albedo from 2000-2009 produced by MODIS data in a resolution of 0.05 is used. Chen et al (2015) performed a detailed validation for the accuracy of the GLASS surface albedo product with ground measurements from seven site observations in the TP. The results showed that the variation tendency of GLASS surface albedo products is consistent with ground measurements, with the biases and root mean square errors (RMSEs) of 0.0184 and 0.0124 respectively. They also compared the GLASS surface albedo product, the MODIS MCD43B3 product and the GlobAlbedo product, indicating that the proportion of GLASS high quality data number is higher than that of MODIS, which is concentrated in more than 80%, while it is 50%-80% for the MCD43B3 surface albedo product. Therefore, it is relatively reliable to use the GLASS surface albedo product in the WRF model in the TP region.

Meteorological observations
The observed monthly climatology of 2 m air temperature and precipitation from 2000-2009 in a spatial resolution of 0.5 • were used to assess the simulation for spatial and seasonal distribution. In addition, observed 2 m air temperature from 83 sites was used to assess the simulated time series. All the observed datasets are produced by the Climatic Data Center, the National Meteorological Information Center, China Meteorological Administration (http://data.cma.cn/site/index.html).
We also used observed surface albedo to compare with the two datasets. Observation was set up at a grassland land surface scape (33 • 53'N; 102 • 8'E) in the east of the TP . Vegetation is composed by Cyperaceae and Gramineae, which is a typical alpine meadow in the TP. The soil is silt clay loam, including 66.7% silt, 29.8% sand and 3.5% clay in the top 40 cm. The annual mean 2 m air temperature is 1.2 • C. The annual mean precipitation is 620 mm. The terrain is flat, homogeneous, and 3423 m a.s.l. Radiation sensors were installed at a height of 3.15 m above the ground.

Model setup
The WRF model is used in this work. WRF was developed for the needs of both weather forecasting and atmosphere related research. It was generated by a collaborative partnership between the National Oceanic and Atmospheric Administration (the National Centres for Environmental Prediction (NCEP), the National Center for Atmospheric Research (NCAR), the Air Force Weather Agency (AFWA), the Forecast Systems Laboratory (FSL)), the University of Oklahoma, the Naval Research Laboratory, and the Federal Aviation Administration (FAA). In this work, we used Version 3.6.1 (Skamarock et al 2008) by using the physics schemes as follows: the Rapid Radiative Transfer Model (RRTM) longwave radiation scheme, the Dudhia shortwave radiation scheme, WRF Single Moment 5c lass microphysics scheme, the Kain-Fritsch cumulus physics scheme, Monin-Obukhov surface layer similarity, the Yonsei University boundary layer scheme and Noah landsurface scheme, with snow albedo parameterized using the Livneh formulation. The model region centered We performed two simulations. The first is the control simulation (hereafter referred to as CTL) in which the default surface albedo was used. Another is the MODIS surface albedo simulation (hereafter referred to as MOD). In the CTL simulation, we used the default monthly surface albedo from AVHRR described above.
In the MOD simulation, as snow cover albedo was parameterized in the WRF_Noah land surface model, which can result in a snow-dependent albedo when and where there is snow, we closed the snow albedo parameterization in MOD and updated the lower boundary condition with the 'real' MODIS surface albedo product.
The MODIS surface albedo product was used to update the WRF lower boundary condition following the works of Meng et al (2014aMeng et al ( , 2014b and Evans et al (2017). Firstly, MODIS surface products were missing, filled by the interpolations between the two nearest valid values. Then MODIS values were reprojected to the WRF Lambert conformal projection to be consistent with the WRF model setup. Thirdly, the MODIS values were linearly interpolated to 6 h intervals and 30 km in spatial resolution to match with the WRF lower boundary condition inputs. At last, the prepared MODIS surface albedo file was used to replace the WRF lower boundary condition inputs in each year from 2000-2009 to update the surface albedo in the MOD simulation.

Comparison between MODIS and WRF default albedo
First, we compared MODIS surface albedo with that initialed by WRF default run. The albedo biases were apparent in boreal winter (December, January, and February (DJF)) and smaller in boreal summer (June, July, and August (JJA)) (figure 1). In the west TP, the WRF default shows much higher values than those of MODIS surface albedo in almost all seasons. In contrast, the biases in the central and east TP are relatively smaller. We calculate an average of the two albedo datasets; the average of the WRF default is 0.36, while it is 0.22 for the MODIS datasets, showing a high overestimation in the WRF default run. This can create a widespread influence on energy budget in the model. It is noted that in the Qaidam basin in the north TP, the biases were shown to be small in all seasons.
Similarly, time traces show a high overestimation of surface albedo in winter; the maximum of the WRF default surface albedo was greater than 0.50 in winter, while it was about 0.28 in the MODIS products (figure 2). Both products can consistently reflect interannual and annual variations, with WRF default surface albedo exhibiting far larger amplitudes than MODIS traces.   figure 3. There is systematic overestimation of surface albedo in the WRF default compared to in the MODIS product, except for the days with snow fall. Biases on snow days in both surface albedo datasets implied that better description for snow albedo is necessary for snow-albedo retrieval or snow-albedo parameterizing, although there is a scale-match issue when using point measurements compared with the gridded model or satellite data.

Comparison of 2 m air temperature simulated by the two runs
We updated the lower boundary condition using the two surface albedo products in the WRF simulations. Seasonal comparisons indicate definite improvements for 2 m air temperature in all seasons (figure 4). The marked region is in the west TP, with prominent improvements in boreal winter, substantially in boreal spring (March, April, and May (MAM)) and in boreal autumn (September, October, and November (SON)) last.
In the time series (figure 5), both simulations can prescribe the interannual and annual variations of 2 m air temperature in the TP, with both simulations underestimating the 2 m air temperature. There was significant improvement for cold biases in the MOD run, indicating the influence of the overestimation of energy. The biases of the CTL and MOD runs are −4.31 and −2.51 • C respectively compared with the observations. The RMSEs of the CTL and MOD runs are 5.25 and 3.24 • C respectively, showing large improvements for cold biases corrections.

Comparison of net radiation simulated by the two runs
Surface albedo differences can directly result in changes in the radiation budget. The simulated patterns of net radiation differences between the two runs (figure 6) were clearly consistent with those in surface albedo. The surface albedo biases led to more than 40 W m −2 maximum in the TP in spring and in winter, which is almost 50% of the annual average net radiation and with the fact that the magnitude in spring was somewhat larger than in winter.

The reason of surface albedo overestimation in the CTL run
In the WRF model, surface albedo was parameterized using the background albedo with snow cover. The parameterization scheme is as follows: where 0 is snow-free albedo, i.e. background albedo, f snow is fractional snow cover, and snow is snow cover albedo; this is calculated by where snow0 is the maximum snow albedo from new snow which is considered as 0.85 in Robinson and Kukla (1985). LV coef is a user-definable coefficient for adjusting snow albedo, which is set in file GENPARM.TBL. The Livneh scheme boosts the snow albedo toward 85%, then reduces it according to the age of the snow. LV coef should range between 0 and 1. Values lower than 0.5 will tune this more toward the incoming snow albedo (and ultimately lower surface albedo). Values greater than 0.5 will tune this more toward 85% (and ultimately higher albedos). In this simulation, it is used as 0.5. In figure 1, two surface albedo products are relatively comparable in summer in which there is little snow. This implies the key influences of snow on surface albedo values in spring, autumn and winter. Figure 7 shows the differences between simulated snow and observed precipitation in different seasons. It displays large positive biases in snow simulation in the TP in spring, autumn and winter. According to the snow parameterization scheme, surface albedo is closely related to both snow-free albedo and simulated snow cover. Figures 8 and 9 show the simulated snow albedo and snow cover fraction respectively. Overestimation of snow cover corresponds with overestimation of snow, which results in substantiallyoverestimated snow albedo and surface albedo.
In table 1, it is evidently shown that there is a large overestimation of precipitation in all models in the TP. In our simulation, large positive biases exist in  both snow and precipitation simulations in the TP. As was shown by previous studies, the bias for precipitation simulation in the TP was 0∼2 mm day −1 in average in WRF simulations in the work of Gao et al 2015a and Yu et al 2015. In our study, however, we found larger biases for precipitation simulation, which is 0∼2.37 mm day −1 . This might be because the time period of our simulation was shorter than previous simulations. After we updated surface albedo, the MOD performed even worse for precipitation simulation with a bias of 0∼2.74 mm day −1 . As we can expect, lower surface albedo can result in more absorption of energy in land surface, which could enhance moist static energy, and then more precipitation in advance (Meng et al 2014a).
In addition, snow albedo was parameterized according to the age of the snow. Due to the colder temperature simulated by the WRF model, snow tends to exist for a longer time. This results in an older age for snow, leading to larger and long-existence snow albedo over the TP. For example, in spring and early winter, large-area snow falls at night or in the morning, but melts quickly in the afternoon due to strong solar radiation, with only very small parts left in the valley or roadsides, etc., which were overestimated in the model almost all the time.    Another important factor that affects surface albedo in the TP is from the direct and indirect influence of light-absorbing aerosols. Based on observations, black soot aerosols deposited on Tibetan glaciers and snow are found to be responsible for surface albedo reduction, and then have a significant contribution to observed rapid glacier retreat and surface energy budget on snow cover (Xu et al 2009, Zhang et al 2017. The influence can exceed 50% in maximum to the surface albedo reduction, which is an nonnegligible aspect in surface albedo parameterization in winter and spring (Zhang et al 2017). A statistical analysis based on remote sensing data indicated that the influence of light-absorbing aerosols on snow albedo might be 0.75 when considering a 1 K increase in land surface temperature and a 0.1 increase in aerosol optical depth (Lee et al 2016). In addition, a previous study showed that the annual mean snow albedo forcing from black carbon is 2.9 Wm −2 averaged over snow-covered regions on the TP, which is a factor of three larger than the value over global land snowpack (He et al 2014). These results suggested ineligible influences on snow albedo from light-absorbing aerosols on the TP.

Discussions
We demonstrated that the overestimated snow causes overestimation of surface albedo and then cold temperature biases in the TP. Firstly, the overestimation of precipitation generally leads to large snow albedo in winter. Secondly, snow-albedo parameterization overestimates surface albedo in the TP in advance due to the age of snow and the simulated cold 2 m air temperature. Lastly, the influence responded by light-absorbing aerosols could be large, which has not been considered in the model. Previous studies have concluded that large uncertainty in snow cover fractions at grid-cell scale in land surface models (Frei and Gong 2005) leads to errors in modeling snow and its interactions with the atmosphere. However, even when the snow cover extent is captured well by MODIS surface albedo datasets, large biases of 2 m air temperature still persist in winter in the MOD simulation. Part of the reason might be the influence of the overestimation for precipitation and the substantial overcalculation for evaporation, which may result in cold air temperature. However, cause of large cold biases remains unclear and requires more comprehensive investigation in the TP in WRF and other models.
Although global numerical forecasting models have been improved when progressing to more realistic weather and climate predictions, those processes related to water vapor, particularly for cloud and precipitation, are still poorly simulated . Xu et al (2010) showed that the ensemble of CMIP3 models overestimates precipitation over the TP by up to 100%. As to CMIP5, Su et al (2013) showed that the models overestimate the precipitation over the eastern TP by up to 62%-183%. Even based on regional climate models, precipitation overestimation is still more than 50%, with the maximum exceeding two times in the TP as was investigated in table 1. The influence of the overestimated precipitation remains a large space for improvement in all models and may have complicated effects on both cold biases and surface albedo estimation in models.
In addition to the influence of snow albedo overparameterization, large cold biases on the TP in the simulations could also be caused by the overestimation of precipitation as was discussed above. The maximum of 2 m air temperature in summer was still underestimated, which was shown in figure 5. Our simulation shows one of the problems that resulted in the cold biases in the simulations, and it is also a feasible way to improve 2 m air temperature simulation. However, much more work still needs to be done in the future to improve both precipitation and air temperature simulation for the TP region.

Conclusions
We updated the surface albedo in the WRF lower boundary condition based on MODIS surface albedo products and demonstrated an evident improvement for cold temperature biases. The large surface albedo overestimation in winter and spring in the WRF model tended to result in large cold temperature biases, which were caused by simulated precipitation biases and overparameterization for the surface albedo calculation. It highlighted the necessity of correctly representing snow albedo in the TP, where snow is not as large as was simulated in the model, and snow melting is quick most of the time due to large solar radiation compared with low-elevation regions. The influence of light-absorbing aerosols should be involved as well, especially with the enhancement of air pollution in recent years. In addition to the snow albedo defining, our results show an urgent call for improving precipitation simulation in the TP.