Simulation of summer climate over Central Asia shows high sensitivity to different land surface schemes in WRF

Land surface processes are vital to the performance of regional climate models in dynamic downscaling application. In this study, we investigate the sensitivity of the simulation by using the weather research and forecasting (WRF) model at 10-km resolution to the land surface schemes over Central Asia. The WRF model was run for 19 summers from 2000 to 2018 configured with four different land surface schemes including CLM4, Noah-MP, Pleim-Xiu and SSiB, hereafter referred as Exp-CLM4, Exp-Noah-MP, Exp-PX and Exp-SSiB respectively. The initial and boundary conditions for the WRF model simulations were provided by the National Centers for Environmental Prediction Final (NCEP-FNL) Operational Global Analysis data. The ERA-Interim reanalysis (ERAI), the GHCN-CAMS and the CRU gridded data were used to comprehensively evaluate the WRF simulations. Compared with the reanalysis and observational data, the WRF model can reasonably reproduce the spatial patterns of summer mean 2-m temperature, precipitation, and large- scale atmospheric circulation. The simulations, however, are sensitive to the option of land surface scheme. The performance of Exp-CLM4 and Exp-SSiB are better than that of Exp-Noah-MP and Exp-PX assessed by Multivariable Integrated Evaluation (MVIE) method. To comprehensively understand the dynamic and physical mechanisms for the WRF model’s sensitivity to land surface schemes, the differences in the surface energy balance between Ave-CLM4-SSiB (the ensemble average of Exp-CLM4 and Exp-SSiB) and Ave-NoanMP-PX (the ensemble average of Exp-Noah-MP and Exp-PX) are analyzed in detail. The results demonstrate that the sensible and latent heat fluxes are respectively lower by 30.42 W·m−2 and higher by 14.86 W·m−2 in Ave-CLM4-SSiB than that in Ave-NoahMP-PX. As a result, large differences in geopotential height occur over the simulation domain. The simulated wind fields are subsequently influenced by the geostrophic adjustment process, thus the simulations of 2-m temperature, surface skin temperature and precipitation are respectively lower by about 2.08 ℃, 2.23 ℃ and 18.56 mm·month−1 in Ave-CLM4-SSiB than that in Ave-NoahMP-PX over Central Asia continent.


Introduction
Central Asia is an extensive arid and semiarid region located in the mid-latitude of Eurasian continent (Lioubimtseva and Henebry 2009;Peng et al. 2018;Jiang et al. 2019). Central Asia is located in the hinterland of Eurasia. It has a typical temperate continental climate of deserts and grasslands because the plateau in its southeast blocks the warm and humid air inflow from the Indian Ocean and the Pacific Ocean (Dando 2005). Central Asia is characterized by three features. First, the rainfall is scarce with annual precipitation generally below 300 mm, thus the region is extremely dry. The annual precipitation near the salt sea and the desert of Turkmenistan is only 75-100 mm, although the annual precipitation in the mountain area is 1000 mm. Second, Central Asia is located in the interior of the midlatitude continent. There are long-lasting sunny days with strong solar radiation, high temperature and large evapotranspiration. Central Asia receives about 12-15 millions W·m −2 solar radiation per year, and the amount for the Turkmenistan can be as large as 19 millions W·m −2 (Peng et al. 2018;Jiang et al. 2019). Third, the diurnal temperature range in Central Asia is large. In most regions, the difference between the daily maximum temperature and minimum temperature is 20-30 ℃. In the Pamir Plateau, the amplitude of diurnal temperature variation can be up to about 40 ℃. The monthly mean temperature is generally between 26 and 32 ℃ in the summer. Therefore, land-atmosphere interaction is crucial in shaping the characteristics of regional climate over Central Asia.
Central Asia is the key connecting region along the Silk Road Economic Belt that extends from East Asia to Europe (Peng et al. 2018). This region is facing severe water shortages and societal conflicts as a result of the dry climate (Peng et al. 2018;Jiang et al. 2019). The region "Central Asia" mentioned in our study refers to the five Asian republics of the former Soviet Union, namely, Kazakhstan, Uzbekistan, Turkmenistan, Kyrgyzstan, and Tajikistan (Lioubimtseva and Henebry 2009;Jiang et al. 2019). This region is controlled by a typical continental climate with little precipitation and high evaporation, with highly sensitive to climate change (Lioubimtseva and Henebry 2009;Schiemann et al. 2008;Jiang et al. 2019;Lai et al. 2020). Additionally, significant shrink of glaciers, reduction of snow cover, and decline of lake sizes have been observed in Central Asia due to global warming and impacts of other anthropogenic activities during recent decades (Pala 2005;Stone 2008;Li et al. 2011;Guan et al. 2018;Jiang et al. 2019). Central Asia is among the regions of highest per capita water use in the world, which apparently mismatches the condition of water deficiency (Varis 2014;Jiang et al. 2019). In-situ observations are scarce in Central Asia, and the spatial resolutions of gridded observational products and reanalysis are commonly coarse. This dampens the study of the regional climate carried out over Central Asia, especially at fine scales. As a powerful tool for downscaling studies, it's very important to investigate whether Regional Climate Model (RCM) is capable to capture more regional climate informatiom at finer scales compared to the coarse observations or reanalysis products. The Weather Research and Forecasting (WRF) model is a state-of-the-art RCM for dynamic downscaling (Liu et al. 2019). However, there are large uncertainties in the WRF model simulations. It is found that the performance of WRF model varies with regions, seasons and climatic conditions, owing to various factors such as the atmospheric circulation natural variability (Kjellström et al. 2011;Li et al. 2016), the physical parameterization schemes (Fernández et al. 2007;Mooney et al. 2013;Li et al. 2016), the size, location and resolution of simulation domain, the lateral boundary conditions (Xue et al. 2014), and the representation of land surface conditions (e.g., topography and land cover, Ge et al. 2019).
Many studies have illustrated that land surface processes have significant impacts on water vapor and momentum exchanges between land and atmosphere, which in turn affect the nature of atmospheric circulation and even causes climate feedback in the long term (Roesch et al. 2001;Li et al. 2016). Therefore, the performance of WRF should be highly dependent of the land surface scheme, which largely determines the land surface processes in RCM. Multiple land surface schemes have been implemented in the more recent versions of WRF model (Idabel et al. 2015;Liu et al. 2019). A number of studies have explored the sensitivity of RCM to land surface schemes from different perspectives and investigated the mechanisms that control land-atmosphere interactions word-wide (Meehl 1994;Douville and Royer 1996;Xue 1996;Texier et al. 2000;Douville et al. 2001;Xue et al. 2004;Wu et al. 2012;Boone et al. 2016;Li et al. 2016;García-García et al. 2020). Zeng et al. (2011) assessed the sensitivity of WRF simulated short-term hightemperature weather to different land surface parameterization schemes. Sato and Xue (2013) showed that diverse descriptions of land surface process are one major factors that contributes to uncertainties in the East Asian summer monsoon simulated by WRF. Li et al. (2016) revealed the choice of land surface schemes have important effects on the simulated East Asian summer monsoon. Liu et al. (2019) evaluated a snow event simulation using various initial and boundary conditions and different land surface schemes over the Tibetan Plateau. Chen et al. (2019) assesses the impacts of the surface coupling strength on regional climate simulation by using 4-km WRF CONUS schemes. García-García et al. (2020) illustrated the importance of the different land surface model choice in climate simulations over North America. However, the studies focusing on the sensitivity of RCM to land surface schemes over Central Asia are rare. Thus, how land surface schemes impact on climate simulation in this region remains unknown.
The resolution of RCM is another crucial factor that affects its downscaling ability. Increased resolution can more realistically represents regional land surface feature, such as topography and underlying surface conditions. Accurate underlying surface information has a significant impact on the accuracy of regional climate simulation. Previous studies (Pan and Li 2011;Shi et al. 2012) showed that RCMs with high resolutions can better simulate the distribution characteristics of the surface air temperature, particularly over complex terrain area. Giorgi and Marinuci (1996) found that the large-scale average precipitation intensity, frequency and surface flux are very sensitive to the resolution, and the RCM with higher resolution can better simulate the extreme drought and flood events. Zhao and Luo (1999) illustrated that the simulated precipitation with high vertical resolution is apparently better than that with low resolution in the regional climate simulation over East Asia. Leung and Qian (2003) found that the RCM using high resolution improves the simulation ability of rainfall in coastal mountainous areas and basins on the long-term summer simulation through two nested methods. Tang et al. (2006) proved that high resolution can enhance the ability of heavy rainfall simulation. With the increase of horizontal resolution, the effect of topography and other local forcing factors is more obvious, and the ratio of local precipitation is also increase in the Yangtze-Huaihe River Basin during Meiyu period. Lű et al. (2009) demonstrated that the increase in both horizontal and vertical resolution can produce better predictions of temperature and precipitation, and the simulation of soil temperature is better with finer horizontal resolution. However, long-term simulation with higher resolution is lack for Central Asia.
This study conducts a higher resolution simulation of summer climate over Central Asia using WRF model, and examines the sensitivity of WRF performance to land surface schemes. This study is organized as follows. In Sect. 2, we introduce the model configuration, the experimental design, and the data together with the statistical method applied to evaluate the simulation by WRF. In Sect. 3, we first analysis the simulated summer mean 2-m temperature, precipitation, and large scale circulation, then we assess the simulation results of the four experiments comprehensively by using a pretty creative method named Multivariable Integrated Evaluation. Finally, we investigate physical processes from the perspective of the energy balance that greatly affects the simulation performance of WRF with different land surface parameterization schemes. Section 4 presents the main results and discussion.

The WRF model configuration
The WRF model version 4.0 is used here. To perform higher resolution simulation, a nested run simulation is applied in this study. The spatial resolutions of the outer and inner domain are 30 km and 10 km, respectively (Fig. 1). The inner domain of the model is centered at 40° N and 60° E, covering an area of 3990 km (west-east) × 5280 km (south-north). There are 28 eta levels in vertical extending from the surface to 50-hPa. We adopted the one-way feedback option without considering the feedback from the inner domain to the outer domain. The NCEP-FNL product is one of the most commonly used data to drive WRF (Kalnay et al. 1996;Li et al. 2009;Liu et al. 2019), and precipitation simulations can be well reproduced (Maussion et al. 2014). The NCEP-FNL analysis data is derived from the Global Data Assimilation System (GDAS), which continuously collects observational data from the Global Telecommunications System (GTS), and other sources, for many analyses. The NCEP-FNL data used in the present study at 6-h intervals

Experimental design
The simulation covers 19 summers from 2000 to 2018 over Central Asia. For each summer, the model integration starts on 23th May and ends on 1st September. The first nine days (23th to 31th, May) are regarded as the spin-up period and discarded. The remaining 3-month simulations (June-July-August) are used for analysis. To investigate the impact of land surface scheme options on simulation, experiments were conducted separately with four different land surface schemes, including the Community Land Model version 4 (CLM4; Lawrence et al. 2011), the Noah Multiple Parametrization land surface model (Noah-MP; Niu et al. 2011), the Pleim-Xiu land surface model (PX; Pleim and Xiu 1995;Xiu and Pleim 2001), and the Simplified Simple Biosphere (SSiB) land surface model (Xue et al. 1991;Sun and Xue 2001). The four experiments are referred as Exp-CLM4, Exp-Noah-MP, Exp-PX and Exp-SSiB respectively hereafter. Other options of physics parametrization schemes are identical in all four experiments (Table 1). Therefore, the simulation differences among these experiments should be attributed to the land surface schemes.
The CLM4 model is composed of complex biogeophysical, biogeochemical, hydrology processes with dynamic vegetation considered. CLM4 contains ten-layer soil column, up to five-layer snow-pack, and single-layer vegetation canopy in the vertical direction. The CLM4 model represents the surface with five primary sub-grid land cover types including the glacier, lake, wetland, urban, and vegetation in every grid unit (Lawrence et al. 2011). The vegetated portion of a grid unit is further divided into 17 plant functional types (PFTs, including the bare soil) and each type is described with distinct physiological parameters.
The Noah-MP land surface model contains a vegetation canopy layer to compute the canopy and the ground surface temperatures separately on the basis of the Noah-V3 model structure (Niu et al. 2011). Noah-MP model adds the two stream radiation transfer scheme considering canopy gaps to compute fractions of sunlit and shaded leaves and their absorbed solar radiation. The Noah-MP model not only has a Jarvis scheme that relates stomatal resistance, but also adds a Ball-Berry type stomatal resistance scheme to photosynthesis of sunlit and shaded leaves. It is significant that Noah-MP includes a short-term dynamic vegetation model, a simple groundwater model with a TOPMODEL-based runoff scheme, a physically based three-layer snow model, and a frozen soil scheme that produces a greater soil permeability. The Noah-MP model with multiple parameterization options has a great potential to facilitate physically based ensemble climate predictions, identification of the optimal combinations of schemes and explanation of model differences, and identification of critical processes controlling the coupling strength between the land surface and the atmosphere.
The Pleim-Xiu model includes a two-layer soil with a 1-cm surface layer and a 1-m root zone layer, surface fluxes including parameterization of vegetation, and a non-local closure planetary boundary layer model (Pleim and Xiu 1995;Xiu and Pleim 2001). In Pleim-Xiu model, the vegetation and soil parameters of each gird cell are derived from fractional coverage of soil texture types and land use categories. The ground surface (1 cm) temperature is computed based on the surface energy balance using a force-restore algorithm for heat exchange within the soil (Xiu and Pleim 2001). Stomatal conductance is parameterized according to root zone soil moisture, air temperature and humidity, photosynthetically active radiation (PAR), and several vegetation parameters such as vegetation leaf index (LAI) and minimum stomatal resistance (Xiu and Pleim 2001).
The SSiB3 model is simplified version of the SiB (Simple Biosphere Model) developed by Xue and Sun (Xue et al. 1991;Sun and Xue 2001). SSiB3 is based on biophysical process, which is designed for global and regional models to describe the interaction between the land and the atmosphere. The parameters of the simplified model SSiB3 are reduced approximately by half (from 44 to 21), contributing to more efficient simulation of the land-atmosphere coupling model. The SSiB3 model contains one-layer vegetation, three-layer snow-pack to accurately and reasonably describe the snow processes such as destructive metamorphism and densification processes resulting from the snow load and melting respectively, two-layer soil temperature (the surface and deep layer), and three-layer soil moisture (0-0.02 m for the surface layer, 0.2-0.5 m for the root layer, and 0.3-2 m for the water gravity leakage layer). In SSiB3 model, the soil temperatures are calculated by the force-restore method. The aerodynamic resistances that controls the momentum, heat and water exchange between the land and the atmosphere is determined by the vegetation characteristics, underlying surface conditions and bulk Richardson number based on the revised Monin-Obukhov similarity theory.

Validation data
To evaluate the WRF model performance, the downscaled large scale circulations of wind, geopotential height and vertical velocity are validated against the ERA-Interim (ERAI) reanalysis dataset. EARI is produced by a data assimilation system based on the Integrated Forecast System developed by the European Centre for Medium-Range Weather Forecasts (ECMWF). EARI covers the data-rich period ranging from 1979 to 2019. Berrisford et al. (2009) and Dee et al. (2011) detailed the ERAI reanalysis product archive, describing its atmospheric forecast model and sequential data assimilation schemes, as well as the performance of this systems. The data we used here has a horizontal resolution of 0.1° latitude by 0.1° longitude and temporal resolution of 6-h. Note that the native spatial resolution of ERAI is 0.75°, and the 0.1° product is produced by interpolation. The data is downloaded from https:// apps. ecmwf. int/ datas ets/ data/ inter imfull-daily/. The National Oceanic and Atmospheric Administration (NOAA) reconstructed dataset of monthly mean 2-m temperature over the land at a horizontal resolution of 0.5° latitude by 0.5° longitude, which covers the period from January 1948 to the present. They are available at https:// www. esrl. noaa. gov/ psd/ data/ gridd ed/ data. temp. html. TEMP/. The data is produced by the observations collected at stations in the Global Historical Climatology Network version 2.0 (GHCN V2.0) and the Climate Anomaly Monitoring System (GHCN-CAMS) dataset, which are interpolated to regular latitude-longitude grids using the cubic spline method. Standard quality controls mainly include the removal of outliers, the interpolation at points where data is missing, and the check for temporal consistency in time series. Digital terrestrial height data (Global 30 arc-second elevation) are used to eliminate the impact of terrain on the accuracy of interpolation. Xue et al. (2004) indicate that the GHCN-CAMS data can well represent the 2-m temperature. It is the best available large scale near surface temperature with high resolution at this point.
The downscaled monthly mean precipitation is evaluated based on the data from the University of East Anglia Climate Research Unit (CRU) with version 4.01, which covers the period from 1901 to 2016. The CRU data are produced by interpolating the observations collected from a number of stations worldwide to 0.5° latitude by 0.5° longitude grids. The dataset is obtained at https:// dx. doi. org/ 10. 5072/ edf8f ebfda ad48a bb2cb af7d7 e846a 86/. The precipitation data from CRU is homogeneity adjusted because spatially interpolating observational precipitation at stations to grids may induce large bias, especially in the regions with sparse stations. Harris and Jones (2017) show that this gridded rainfall data is highly reliable as a result of a height correction over Central Asia. Therefore, we verify WRF model simulation of precipitation against CRU data in this study.

Model assessment
To better evaluate the downscaling ability of the WRF, three statistical metrics including mean bias (BIAS), spatial correlation coefficient (SCC) and root-mean-square error (RMSE) are applied. Meanwhile, we calculate the statistical significance of the BIAS, SCC and RMSE through a bootstrap method (Efron 1979;Li et al. 2016). The bootstrap method estimates the distribution properties of the quantity a Monte Carlo simulation, in which a distribution of the quantity is created by sampling the set of measurements with replacement many times. The main advantage of using the bootstrap resampling approach is that good estimates can be obtained, regardless of the complexity of the data processing. In present study, resampling with replacement is applied from 1000 bootstrap to estimate the probability distribution of the multi-year averaged fields. When the difference is considered between two data sets and if the significance level of the difference is of interest, the variance for two data sets have to be combined to obtain the standard error of the difference. The method of combination depends on whether the variances of the two data sets are equal. Statistical tests exist for determining the equality of the two data sets (Cheng and Yeager 2007). Hence, the BIAS, SCC and RMSE three statistical measures were calculated for each bootstrap. Then the results are sorted. The 5% significance intervals are determined from the 2.5% and 97.5% of the values obtained by 1000 repetitions. The detailed steps to design the significance test and determine the confidence interval of RMSE (also BIAS and SCC) between two experiments can be found in Cheng and Yeager (2007), Zhang (2011) and Peng et al. (2015).
In addition, a Multivariable Integrated Evaluation (MVIE) method is applied to evaluate the overall model performance in simulating multiple fields. The core concept of MVIE is to group various scalar fields into a vector field and compare the constructed vector or scalar field with the observations. The MVIE method can be flexibly applied to full fields (including both the mean and the anomaly) or anomaly fields depending on the application. The detailed introduction and application of this method can be found in Xu et al. (2016Xu et al. ( , 2017 and Huang et al. (2019).

2-m temperature and precipitation
The summer mean of observed and simulated 2-m temperature and differences between the four simulation experiments and GHCN-CAMS respectively in the period of 2000-2018 are presented in Fig. 2. All the four experiments can reasonably capture the spatial pattern of 2-m air temperature, while there are large differences among the simulated results. In the northern Kazakhstan over Central Asia, Exp-Noah-MP and Exp-PX obviously overestimate the 2-m temperature, whereas simulations of Exp-CLM4 and Exp-SSiB are in closer agreement with the GHCN-CAMS data relatively. The BIASs, SCCs and RMSEs of simulated 2-m temperature from simulations verified against GHCN-CAMS data over Central Asia are displayed in Table 2. Averaged over the study region, Exp-Noah-MP and Exp-PX exhibit warm biases, while Exp-CLM4 and Exp-SSiB exist cold biases. It is clear that 2-m temperatures simulated by Exp-Noah-MP and Exp-PX are respectively about 0.26 ℃ and 0.21 ℃ greater than the GHCN-CAMS data, with 2-m temperature simulated by Exp-CLM4 and Exp-SSiB are only respectively 0.13 ℃ and 0.14 ℃ lower than the GHCN-CAMS data. The SCCs for Exp-Noah-MP and Exp-PX are smaller than that of Exp-CLM4 and Exp-SSiB, with the SCCs for Exp-Noah-MP and Exp-PX against GHCN-CAMS data being respectively 0.96 and 0.94, whereas the SCCs for Exp-CLM4 and Exp-SSiB against GHCN-CAMS data being respectively 0.97 and 0.96. At the same time, the RMSEs of 2-m temperature for Exp-Noah-MP (2.30 ℃) and Exp-PX (3.07 ℃) simulations are higher than that of Exp-CLM4 (2.06 ℃) and Exp-SSiB (2.09 ℃). Obviously, the BIASs, SCCs and RMSEs for Exp-CLM4 and Exp-SSiB differ from that for Exp-Noah-MP and Exp-PX at the 0.05 statistical significance level.
Central Asia is facing severe water deficiency and conflicts (Peng et al. 2018). Precipitation is one of the most important components of the water cycle, and its changes can substantially affect sustainability economic and social developments (Gillett et al. 2004;Peng et al. 2018).   Figure 4 displays the summer mean wind vector (UV850) and geopotential height (H850) during 2000-2018. It is shown that the long term mean structures of UV850 and H850 of the ERAI and the NCEP-FNL reanalysis are similar. The WRF can also capture well the spatial pattern of 850-hPa circulation. The simulated UV850 and H850 of the four experiments agree well with that of the ERAI and NCEP-FNL reanalysis, although the simulated winds are slightly weaker in intensity. Table 3 lists the BIASs, SCCs and RMSEs of U850 and H850 over the simulation domain. The BIASs, SCCs and RMSEs in Exp-CLM4 and Exp-SSiB are significantly different from that in Exp-Noah-MP and Exp-PX at the 0.05 significance level. The statistical evaluation shows that U850 and H850 in Exp-CLM4 and Exp-SSiB simulations closely match reanalysis than that in Exp-Noah-MP and Exp-PX. As can be seen from Table 3, U850 simulated by Exp-CLM4 and Exp-SSiB are respectively 0.22 m·s −1 and 0.12 m·s −1 lower than the ERAI reanalysis data, while U850 simulated by Exp-Noah-MP and Exp-PX are separately 0.28 m·s −1 and 0.27 m·s −1 lower than the ERAI data. The SCCs of U850 and H850 for the In order to show more clearly the deviations between the simulation results and the reanalysis data, we analyze the differences of lower-level summer mean geopotential height and wind vector at 850-hPa between the four simulation experiments and ERAI reanalysis during 2000-2018 respectively (Fig. 5). All the four simulations also well capture the spatial pattern of UV850 and H850 overall. However, large positive biases of H850 are found near southwestern Kazakhstan and the northeastern Aral sea in Exp-CLM4 and Exp-SSiB simulations. Meanwhile, H850 is overestimated over most Central Asia in the four experiments. Results verify that the simulations of H850 in all the four experiments have significant positive biases. However, note that the biases of the geopotential height over Central Asia are small in Exp-CLM4 and Exp-SSiB (Fig. 5a, d) than that in Exp-Noah-MP and Exp-PX (Fig. 5b, c). As can also be seen from Table 3, the absolute values of BIASs of Exp-Noah-MP and Exp-PX simulation results are higher than that of Exp-CLM4 and Exp-SSiB at the significance level of 0.05. Both Exp-Noah-MP and Exp-PX produce large biases of H850 around the southwestern Kazakhstan between the Caspian sea and the Aral sea. Same as the 850 hPa, the summer mean wind vector (UV500) and geopotential height (H500) simulated by Exp-Noah-MP are similar to that of Exp-PX, while results of Exp-CLM4 resemble that of Exp-SSiB.

Large-scale circulation
The above analyses and comparisons of downscaled summer mean 2-m temperature, precipitation and large scale atmospheric circulation among simulation results of the four experiments and that from the validation data demonstrate that the performances of WRF is sensitive to land surface schemes over Central Asia. Moreover, the WRF model simulations with higher resolution actually add more detailed information to Central Asia compared to reanalysis data that provides the initial and lateral boundary condition for the WRF simulations (not show). The simulation result of 2-m temperature is more similar to the GHCN-CAMS data in Exp-CLM4 and Exp-SSiB than in Exp-Noah-MP and Exp-PX. The similar performance of Exp-Noah-MP and Exp-PX is probably because both of them use Jarvis scheme to calculate the leaf stomatal resistance, which is a dominant factor of the evapotranspiration. Figure 6 presents the MVIE diagram to quantitatively and comprehensively evaluate the model performance in simulating summer mean 2-m temperature, precipitation, vector winds, and geopotential height at each level in the four experiments. Figure 6a is a generalized Taylor diagram, which shows how much the overall root mean square deviation of different fields can be attributed to the differences in the root mean square value and how much is due to the poor pattern similarity. The statistics include similarity coefficient (SC) and normalized root-mean-square vector length (RMSL) between the ERAI reanalysis/GHCN-CAMS/CRU data and simulations. As seen from Fig. 6a, the SCs of various experiments range from 0.86 to 0.94, indicating different simulation capability of reproducing the spatial pattern of multivariable over the simulation domain in the four experiments. The normalized RMSLs are smaller than 1 for each variable, which indicates that the four experiments systematically underestimate the magnitude of vector winds and geopotential heights in the four pressure levels over Central Asia. Results of the four experiments generally show better performance at the 850 hPa than at other levels (Fig. 6a).

Multivariable integrated evaluation (MVIE)
In addition, the model performs very well in the simulations of 2-m temperature and precipitation in all the four experiments. Meanwhile, the results of SCs and normalized RMSLs in Exp-CLM4 and Exp-SSiB are more consistent with observations compared to that of Exp-Noah-MP and Exp-PX.
Although larger SCs and smaller normalized RMSLs generally represent better correspondence in the model, they do not increase or decrease monotonically following the improvement of model performance (Xu et al. 2017;Huang et al. 2018). We therefore computed the model skill scores, which satisfy the monotonic relationship with the model performance in the four experiments. Figure 6b summarizes  Table 3 Mean bias (BIAS, unit: m·s −1 or m), spatial correlation coefficient (SCC) and root-mean-square error (RMSE, unit: standard deviation) of summer mean zonal wind at 850-hPa (U850) and geopotential height at 850-hPa (H850)

Land surface heating
The downscaling capability of the WRF model for Central Asia simulation is sensitive to different land surface schemes as shown above. The representation of land surface parameters, such as the surface roughness, stomatal resistance, and leaf area index as well as the photosynthesis process, plays a crucial role in the energy partitioning and balance in the land surface schemes. In the meantime, the description of physical processes and different land surface model components are also critical for the energy partitioning and balance. What is more, the surface roughness is one of the most important parameters to determines the surface wind flow, especially for the surface dynamic roughness. Central Asia region includes three underlying surface types (bare soil, desert and semi-desert, and grassland) at least, the surface roughness of each land cover is different in four land surface schemes so that there are slight differences of the surface wind among four experiments in this study. To further understand the possible cause for this sensitivity, physical processes associated with the land surface schemes are discussed. The land surface schemes influence the surface energy balance, which consequently influences the atmospheric conditions via the land-atmosphere interactions. In the WRF model, surface sensible heat flux and latent heat flux are obtained in the land surface scheme, which are then propagated to the planetary boundary layer scheme and affect atmospheric circulation over Central Asia. The discrepancies of the surface energy budgets and their impacts on atmospheric circulation are compared between Ave-CLM4-SSiB (the ensemble average of Exp-CLM4 and Exp-SSiB) and Ave-NoahMP-PX (the ensemble average of Exp-Noah-MP and Exp-PX) because Exp-CLM4 resembles Exp-SSiB and Exp-Noah-MP resembles Exp-PX as discussed above. Figure 7 and Table 4 present the land surface energy balance differences between Ave-CLM4-SSiB and Ave-NoahMP-PX simulations. Table 4 only shows the values averaged over the land grid points. There are obvious difference in sensible heat flux and latent heat flux (Fig. 7a, b) between Ave-CLM4-SSiB and Ave-NoahMP-PX. The difference can be up to 35 W·m −2 in some areas. In general, the simulation result of sensible heat flux in Ave-NoahMP-PX is much greater than that in Ave-CLM4-SSiB over the land, especially in the Uzbekistan area. As shown in Table 4, .86 W·m −2 higher respectively in Ave-CLM4-SSiB than that in Ave-NoahMP-PX. There are smaller differences in the ground heat flux compared to that in SH and LH (Fig. 7c). In the southern part of Central Asia, the differences of the latent heat flux are substantially small while the differences of the sensible heat flux are relatively large. The Bowen ratio, representing energy partitioning between the SH and LH, is smaller in Ave-CLM4-SSiB compared with Ave-NoahMP-PX (Fig. 7d). Bowen ratio has important impacts on the atmospheric circulation as illustrated in many previous studies (Xue et al. 2004(Xue et al. , 2006(Xue et al. , 2010a(Xue et al. , 2010bFischer et al. 2007;Berg et al. 2014).
The energy source of SH and LH is the net solar radiation absorbed by the land surface. The net radiation flux is composed of net shortwave radiation (NSR) and net longwave radiation (NLR) fluxes. Compared to NSR (− 15.63 W·m −2 ), the difference in NLR (− 4.82 W·m −2 ) is relatively small (Fig. 7e, f; Table 4). Therefore, the difference in net radiation (NR) mainly depends on the difference in NSR (Fig. 7e, g; Table 4). The NSR is strongly affected by surface albedo (Fig. 7h; Table 4). The lower albedo (0.05) in Ave-NoahMP-PX results in a higher NSR and thus a higher NR on the surface compared to that in Ave-CLM4-SSiB, which partially leads to the higher SH. As seen from Table 4, the NR difference (− 20.45 W·m −2 ) is balanced by reduced SH (− 30.42 W·m −2 ) and enhanced LH (14.86 W·m −2 ). We have compared two ensemble average results (SH; LH; ground heat flux; Bowen ratio; NSR; NLR; NR; albedo) with Global Land Data Assimilation System (GLDAS) data to evaluate the performance of two groups of WRF simulations (Figures not shown). It can be seen that both group simulations can capture well the spatial pattern of land surface energy components, with the Ave-CLM4-SSiB more close to the GLDAS data than the Ave-NoahMP-PX. This suggests the land surface process influences the surface energy balance by regulating the energy partitioning between SH and LH fluxes in addition to the budget of surface radiation.
Land surface processes have direct and significant impacts on 2-m temperature, which can result in high temperature weather under special conditions in summer (Fischer et al. 2007;Berg et al. 2014;Li et al. 2016). Ave-NoahMP-PX overestimates 2-m temperature and surface skin temperature in the most areas of the land, especially in the northern Tajikistan, while there is no such overestimation in Ave-CLM4-SSiB simulation results. The simulations of 2-m temperature and surface skin temperature in domain averaged are respectively higher about 2.08 ℃ and 2.23 ℃ in Ave-NoahMP-PX than in Ave-CLM4-SSiB over Central Asia continent (Fig. 8c, f). The surface energy balance and energy partitioning directly influences 2-m temperature and surface skin temperature. Furthermore, large scale atmospheric circulations can also contribute to the variations of 2-m temperature, which will be analyzed in next section.

Influence on atmospheric circulation
In Sect. 3.4, we have discussed the differences in the SH and LH between Ave-CLM4-SSiB and Ave-NoahMP-PX simulation results. The LH is released with the phase changes of atmospheric moist, mostly vaporization and condensation (Xue et al. 2014;Li et al. 2016), and the temperature of the atmosphere is directly influenced by the SH. The absolute difference value of the SH between Ave-CLM4-SSiB and Ave-NoahMP-PX is approximately twice that of the LH (Table 4), suggesting that the difference in total upward heat flux is mainly dominated by the difference in the SH, and have an important impact on the air temperature difference between simulation result Ave-CLM4-SSiB and Ave-NoahMP-PX. The simulated temperature difference at 850 hPa (T850) over Central Asia is shown in Fig. 8a. T850 is lower in Ave-CLM4-SSiB than that in Ave-NoahMP-PX, with the difference in domain averaged T850 is 2.06 °C over land (Table 4). Such a large difference over the land can significantly affect the thermal contrast between the land and the sea because the difference in the sea is negligible (Xue et al. 2004;Li et al. 2016). Li and Yanai (1996) has indicated that the thermal contrast between the land and the sea influences the intensity of monsoon. Compared with Ave-CLM4-SSiB result (Fig. 8a), Ave-NoahMP-PX overestimates T850 over land in Central Asia, resulting in a abnormal thermal difference between the land and sea. The thermal contrasts contribute to the difference of the H850 covering the entire Central Asia between Ave-CLM4-SSiB and Ave-NoahMP-PX (Fig. 8b). The differences in geopotential height and resulting differences in pressure gradient between Ave-CLM4-SSiB and Ave-NoahMP-PX lead to wide-spread differences in the geostrophic winds (Fig. 8b). It is also obvious that the intensified southwesterly winds in Ave-CLM4-SSiB brings warm air from the tropics to Central Asia in Fig. 8b, suggesting that the atmospheric water vapor transport in Ave-CLM4-SSiB is enhanced. Consequently, precipitation and surface soil moisture are increased in Ave-CLM4-SSiB over Central Asia continent. The simulations of domain averaged precipitation and surface soil moisture are respectively lower 18.56 mm·month −1 and 0.046 m 3 ·m −3 in Ave-NoahMP-PX than in Ave-CLM4-SSiB (Fig. 8d, e). The stronger intensity of SH in Ave-NoahMP-PX produces higher tropospheric temperature Fig. 8 Same as Fig. 7, but for a temperature at 850-hPa (T850; unit: ℃); b geopotential height (H850; shaded, unit: m) and wind (UV850; vectors, unit: m·s −1 ) at 850-hPa; c 2-m temperature (unit: ℃); d pre-cipitation (unit: mm·month −1 ); e surface soil moisture (unit: m 3 ·m −3 ); f surface skin temperature (unit: ℃). The green dots in a-f denote the areas passing the 95% significance test over the land, thus the thermal gradient between twogroup simulations is enhanced. As a result, the gradients of geopotential height over Central Asia is also stronger in Ave-NoahMP-PX than in Ave-CLM4-SSiB, consistent with the large heating gradient, which produces stronger southwesterly in Ave-NoahMP-PX simulation (Fig. 9a, b). Ave-CLM4-SSiB produces relatively lower pressure in the north and relatively higher pressure in the south (Fig. 9c). An anomalous southwestward and southward winds in case Ave-CLM4-SSiB would therefore be produced while Coriolis forcing balanced the pressure gradient force based on the adaptive modulation of geopotential height gradients, particularly at lower levels. Consequently, the simulations of downscaled large scale atmospheric circulation, 2-m temperature, precipitation, surface soil moisture and surface skin temperature are all affected (Xue et al. 2004;Li et al. 2016;Sugimoto and Takahashi 2017).

Summary and discussion
Our study proves that WRF is a useful and promising tool for dynamic downscaling. However, it does show certain uncertainties in the simulation with different land surface schemes over Central Asia. There are few studies paying attention to the sensitivity of the WRF model with higher resolution to land surface schemes and evaluation of long Fig. 9 Summer mean sensible heat flux (SH; shaded, unit: W·m −2 ) and 850-hPa geopotential height (H850; contours, unit: gpm): a Ave-CLM4-SSiB; b Ave-NoahMP-PX; c differences between Ave-CLM4-SSiB and Ave-NoahMP-PX. The yellow dots in c denotes the areas passing the 95% significance test term simulations over Central Asia comprehensively. The effects of land surface schemes on surface energy balance are unclear, and the long-term simulation with higher resolution is still lack over Central Asia. Taking advantage of the recently available higher resolution reanalysis data and a pretty creative MVIE method, the sensitivity of higher resolution WRF model to land surface schemes in simulating summer climate over Central Asia is investigated in this study. We design four experiments applying WRF coupled with different land surface schemes. The comparisons of downscaled 2-m temperature and precipitation, as well as large scale atmospheric circulation between all the four experiments simulation results and that from ERAI reanalysis data, GHCN-CAMS and CRU data averaged cover the 19 summers from 2000 to 2018 suggest that the downscaling capability of WRF model over Central Asia is sensitive to the land surface schemes. Additionally, the WRF model simulations actually gain more detailed information over Central Asia compared to reanalysis data that provides the initial and boundary conditions for the WRF simulations. The 2-m temperature and precipitation in Exp-CLM4 and Exp-SSiB simulations are more similar to GHCN-CAMS and CRU data than the simulations of Exp-Noah-MP and Exp-PX results. By using different evaluation methods, it is found that the performance of Exp-CLM4 is consistent with that of Exp-SSiB, while the Exp-Noah-MP performance is consistent with that of Exp-PX.
The present study also aims to deepen our understanding of the physical and dynamical mechanisms for the influence of land surface schemes on the WRF model for the dynamic downscaling over Central Asia. The differences in the surface energy balance simulations between Ave-CLM4-SSiB and Ave-NoanMP-PX simulations and their subsequent impacts on atmospheric circulations are investigated. The simulation results illustrate that the intensity of the simulated SH over Central Asia is pretty stronger in Ave-NoahMP-PX than in Ave-CLM4-SSiB. The thermal contrast is increased because of the higher temperature gradient between two-group simulations over the land. Wind fields are subsequently affected because of the geostrophic adjustment process. Consequently, the simulations of large scale atmospheric circulation, 2-m temperature, precipitation, surface soil moisture and surface skin temperature are all affected. Meanwhile, the analysis of the ensemble mean results indicates that the Ave-CLM4-SSiB simulation is better than the single-model simulations of Exp-CLM4 and Exp-SSiB, which indicates that the ensemble study may be a useful method for the studies using WRF model or other regional climate models as shown in previous studies (Xue et al. 2010a, b;Liang et al. 2012;Mearns et al. 2012;Li et al. 2016).
The surface water and energy fluxes play a significant role in regional climate model dynamic downscaling simulation over Central Asia, especially the radiation and the sensible heat flux. Hence, more studies are required to answer the question whether a realistic simulation of land surface fluxes would lead to a better dynamic downscaling over Central Asia. It's difficult to validate the land surface flux simulations directly in the WRF model or other regional climate models because of lacking of in-situ observations. We expect more station observational data will be available to verify regional climate model simulations over Central Asia. In addition, the initial and lateral boundary conditions are prerequisite for reliable dynamic downscaling. Therefore, more realistic high-resolution reanalysis data and simulation experiments are necessary.