Complex terrain, commonly represented by mountainous region, occupies nearly one-quarter of the Earth’s continental areas. An accurate understanding of water cycle, energy exchange, carbon cycle, and many other biogeophysical or biogeochemical processes in this area has become more and more important for climate change study. Due to the influences from complex topography and rapid variation in elevation, it is usually difficult for field measurements to capture the land-atmosphere interactions well, whereas land surface model (LSM) simulation provides a good alternative. A systematic review is introduced by pointing out the key issues for land surface processes simulation over complex terrain: (1) high spatial heterogeneity for land surface parameters in horizontal direction, (2) big variation of atmospheric forcing data in vertical direction related to elevation change, (3) scale effect on land surface parameterization in LSM, and (4) two-dimensional modelling which considers the gravity influence. Regarding these issues, it is promising for better simulation at this special region by involving higher spatial resolution atmospheric forcing data which can reflect the influences from topographic changes and making necessary improvements on model structure related to topographic factors. In addition, the incorporation of remote sensing techniques will significantly help to reduce uncertainties in model initialization, simulation, and validation.

1. Introduction

The global climate system is an interactive and coupled system consisting of five major components: the atmosphere, the hydrosphere, the cryosphere, the land surface, and the biosphere, forced or influenced by various external forcing mechanisms. As one of the major components located at the boundary between the atmosphere and the lithosphere, land surface significantly affects the energy and mass exchanges with lower boundary atmosphere. By controlling the surface energy balance and water balance, the land surface processes strongly influence the climate change from local to regional and global scale. The changes in key land surface variables (e.g., surface albedo, surface roughness length, soil moisture, land surface temperature, and land cover) will lead to variations in climate and changes in the sensitivity of climate to other disturbances. In turn, the changes in climate will provide feedbacks to land surface and govern the land surface hydrological and ecological processes. Therefore, land surface is attracting increasing attentions from the scientific society of hydrology, atmosphere, environment, remote sensing, ecology, and many others [1, 2].

Obviously, land surface plays a distinctive role in climate system. Many organizations or communities have made many efforts to establish land surface observations system [3, 4]. For example, surface observations are one of the key observation components of the Global Observing System (GOS) (http://www.wmo.int/pages/prog/www/OSY/GOS.html) of World Meteorological Organization (WMO). Composed of about 11,000 stations at or near the Earth’s surface, it provides atmospheric parameters such as atmospheric pressure, wind speed and direction, air temperature, and relative humidity with 3-hourly or even hourly frequency. A regional and global network of micrometeorological tower sites (FLUXNET) is also established around the world (http://fluxnet.ornl.gov/) [5]. Currently, more than 650 tower sites from about 30 regional networks across five continents are operated on a long-term and continuous basis to measure the exchanges of carbon dioxide (CO2), water vapor, and energy between terrestrial ecosystems and the atmosphere by using Eddy covariance method [6] for different landscapes. Except the discrete site observations on land surface, the rapid development of satellite technology enhances the Earth observation at regional or global scale. The National Aeronautics and Space Administration (NASA) proposed the program named Earth Observation System (EOS) for long-term global observations of the land surface, biosphere, atmosphere, and oceans of the Earth by a series of artificial satellite missions and scientific instruments in Earth orbit. The observations include quantitative derivation of land surface key parameters (e.g., land surface temperature, land surface water, energy fluxes, and carbon flux) and monitoring of their spatiotemporal variation. Various land products from different remote sensing observations [710] greatly accelerate land surface studies related to land resource, environment, ecology, climate, natural disaster, and so forth. In addition, numerous land surface experiments, focusing on different scopes, were conducted by combining field measurements with airborne or spaceborne observations to explore the land surface-atmosphere interactions [1113].

Generally, the established land surface observation system from site measurements to remote sensing methods enables us to better understand the conditions of the land surface radiative, ecological, hydrological, and biological processes. Because land surface processes are coupled system with close connection, it is important to capture the interaction among the processes and the feedback effects of land surface perturbations on climate. Many evidences have found that the changes in land surface conditions will bring changes in climate from regional to global scale [14]. The development of land surface model (LSM) since the 1960s provides one way to help us understand the complex interactions (including biophysical, hydrological, and biogeochemical interactions) between land surface and the atmosphere at micro- and mesoscales. Through numerical parameterization of land surface, it can provide a simple and realistic way to show the transfer of mass, energy, and momentum between land surface and atmosphere. The basic task of LSM is to accurately simulate the partitioning of surface net radiation at the land surface into the emission of thermal infrared radiation to the atmosphere, latent heat flux associated with evaporation and transpiration, sensible heat flux, and ground heat flux. Furthermore, it provides the relevant information on land surface and climate data. Because LSM shows increasing importance for both weather forecasting and climate models, physical modeling of land surface processes using LSM plays a key role, not only in large-scale atmospheric models (including general circulation models (GCMs)), but also in regional and mesoscale atmospheric models [1517].

Although LSM simulation usually concerns the land surface processes in regional or mesoscale to couple with climate models, few researches have been done for local scale with big spatial heterogeneity. A typical area is the complex terrain, commonly represented by mountainous region with great topographic changes within a short distance, whose distinct characteristics are of high relative relief, steep slopes, and remarkably biological diversity. In addition, it is usually described as the water tower and climate regulator of the planet [18]. It plays a huge part in the maintenance of hydrological cycle and helps to protect the fresh water sources and climatic stability for the lowlands and all global ecosystems. The ecosystem in this area is essential to the survival of the global ecosystem but at high fragile level. The global change and intensive human activities increasingly influence the capacity of the ecosystem to provide ecological goods and services, maintain the ecoenvironmental balance, and make it more susceptible to accelerated soil erosion, landslides, and rapid loss of habitat and genetic diversity. In view of the fact, it is required to conduct researches for a better understanding of the function of complex terrain ecosystem and its feedback from the external environmental changes. However, the violent elevation change of this region induces high levels of spatial heterogeneity in major land surface properties, such as vegetation, surface roughness, albedo, emissivity, and soil properties, which in turn leads to strong spatial variation of surface fluxes of heat, moisture, and momentum through the atmospheric boundary layer [19]. This change results in the large spatial variation in surface net radiation distribution, surface land cover condition, surface water condition, and atmospheric environment. These situations make the land surface processes over complex terrain present higher complexity than those of flat area. Therefore, it is necessary for LSM to accurately present the interaction between land surface and atmosphere over complex terrain, address the influences of climate change and land cover change on surface-atmosphere energy-water-mass exchanges, and monitor the dynamic change related to the ecosystem environment and mountain hazards.

In this paper, we concentrate on land surface processes simulation by reviewing the development of land surface model and the applications of LSM over complex terrain. Key issues of complex terrain land surface modelling are addressed, and the suggestions for better simulation are proposed finally.

2. Land Surface Model

2.1. Structure

In general, LSM usually consists of four major components as shown in Figure 1. Parameters include the major information about land surface characteristics which are important for the model initialization and operation. Atmospheric forcing data are the driving forces for the energy and water exchange between land surface and atmosphere. For model structure, the physical, chemical, and biological processes connected with different cycles (surface energy, water, and carbon) are usually included with different process definitions and connectivity. Under the support of the above components, a range of processes connected with land surface dynamic information at a range of time scales from hourly to yearly scales (Figure 2) is simulated including the diurnal cycle of the biogeochemical fluxes, pools, vegetation/soil temperature, water content, daily or weekly dynamics of ecosystem information (phenology and nutrient cycle), yearly variation of vegetation dynamics, and land use/cover change.

2.2. Development

In the history of LSM development, LSMs have been improved substantially both in mechanism structures and numerical techniques. LSMs can differ in many ways: the representation of subgrid spatial variability in soil and vegetation characteristics, topography, and water storage, treatment of surface runoff, and parameterization of snow melt, and so forth. Generally, there are three typical generations [14, 21].

The first generation corresponds to the traditional “bucket model” [22]. The bucket model is the simplest LSM based on a simple energy balance equation where the soil heat flux is ignored. Soil layer in this model is assumed to have the same soil depth and fixed water capacity globally, and the soil properties are spatially constant. The evaporation only occurs under the condition that the soil water content is above a threshold value and the runoff will be produced when the soil moisture exceeds the saturation limitation and further precipitation still continues. The transpiration of vegetation is not considered in the model, and the uptake of soil water by roots from the deep soil is also ignored. As concluded by Pitman [14], the common simplification used to simulate evaporation is the major conceptual limitation of the first generation model. The role of turbulence is represented through a bulk aerodynamic resistance. The surface resistance is omitted and the evaporation is controlled by one special variable () which represents the soil moisture availability. It is a simple function of soil water content and changes from zero for dry soil to one for saturated soil.

To overcome the limitations in the “bucket model,” Deardorff [23] implemented an efficient time-dependent equation to predict ground surface temperature based on the “force-restore” method of Bhumralkar [24]. The involvement of two soil layers in the model enables the simulation of heat conduction inside soil. Meanwhile, a single vegetation layer which is regarded as a pervious layer in the “bucket model” is also considered. The explicit vegetation treatment becomes a milestone for the LSM development. The physically based land surface parameterization represents the second generation LSM. It is the basis of numerous LSMs developed later with kinds of improvements on scheme components, data input, land surface parameterization, and mathematical method, such as the original Simple Biosphere (SiB1) model [25], Biosphere-Atmosphere Transfer Scheme (BATS) [15], Noah land surface model [16], and the land surface model Interactions between Soil, Biosphere, and Atmosphere (ISBA) [26]. For example, SiB1 has three soil layers and two vegetation layers (one canopy layer and one ground cover layer), while Noah LSM separates the vertical profile of soil into four soil layers with lower boundaries at 10 cm, 40 cm, 100 cm, and 200 cm below the surface and has one canopy layer and snow layer above the surface. Additionally, a simplified version of SiB1 (SSiB) was developed by Xue et al. [27] to reduce physical parameters and improved computational efficiency. In the second generation LSM, the interactions in the soil-vegetation-atmosphere system are well illustrated. The different properties of soil and vegetation and the 3D structure of the canopy induce big variations in key land surface processes, such as radiation absorption, momentum transfer, biophysical control of evapotranspiration, precipitation interception and interception loss, soil moisture availability, and insulation [21]. Compared with the first generation, the second generation is more physical. The inclusion of hydrologic cycle and vegetation layer enable the second generation model to explicitly represent radiative transfer, turbulent processes above and within the canopy, and the physical and biological controls on evapotranspiration.

In the above two generations, the simulations mainly focus on the energy and water exchanges between land surface and lower atmosphere interface. Since the late 1980s, global climate change has become a more and more popular issue. Many scientists realize that the carbon cycle and exchange with the atmosphere play a pivotal role in the dynamics of the Earth system [28, 29]. However, the previous models are far from satisfying the need of further studies. The progresses in plant physiological research advance the development of biochemical models on leaf photosynthesis by combining photosynthesis and transpiration through a model for stomatal resistance such as the semiempirical model of leaf conductance [30]. The physiological control of evapotranspiration by plants seems to act as an optimization mechanism that seeks to maximize carbon fluxes by photosynthesis and reduce the water loss from the plant by closing stomata. Therefore, the third generation LSM incorporates the biological control of evapotranspiration. In these models, the biochemistry of photosynthesis is linked with the biophysics of stomatal conductance. Photosynthesis and conductance in leaf level are scaled to the plant canopy level based on the optimal allocation of nitrogen and photosynthetic capacity in relation to light availability [31]. These biochemical models improve the LSMs to describe the biological and chemical processes, enhance the Bowen ratio simulation, and simulate the fluxes of carbon, water, and heat simultaneously. As concluded by Sellers et al. [21], the advantages of the third generation models are as follows: the models are more realistic by using photosynthesis-conductance model to couple the fluxes of energy, water, and carbon simultaneously; an important parameter governing the surface fluxes can be obtained; (3) the models can represent realistic feedback to changes in atmospheric CO2.

2.3. LSM Running Mode

The parameterization of land surface processes in LSMs is very important for numerical weather prediction or climate simulation. The implementation of climate model requires the accurate determination of the radiation and energy and water fluxes and momentum at the land surface boundary. In contrast, land surface processes simulation needs the near surface atmospheric parameters such as near surface air temperature, humidity, wind speed, wind direction, air pressure, and precipitation.

According to the connection between climate model and LSM, LSM is initiated to couple with climate models at different spatial scales to study land-atmosphere interactions and feedback, evaluate sensitivity to perturbations, and sort out cause-effects. The coupling is typically done in this way: the atmospheric model provides incident solar and terrestrial radiation and meteorological variables (temperature, wind, humidity, precipitation, and pressure) from the lowest model level to LSM; and the LSM feeds the atmospheric model with heat and moisture fluxes, surface stress, and radiation at the land surface. Currently, more complex planetary boundary layer (PBL) schemes are sensitive to surface fluxes, while cloud/cumulus schemes are sensitive to the PBL structures. Numerical weather prediction (NWP) models increase their grid spacing (1 km and sub-1 km). There is a great need to capture mesoscale circulations forced by surface variability in albedo, soil moisture/temperature, land use, and snow [32]. As the importance of interactions between land and atmosphere has been increasingly recognized, LSM becomes a key component with an increasingly important role for climate models not only for numerical weather prediction but also for regional climate simulations. In previous coupling examples, the implementation of SiB1 for use within GCMs showed that there are big improvements in the simulation of the continental fields of absorbed radiation, evapotranspiration, and precipitation [33]. Chen and Dudhia [32] coupled an advanced land surface model (previous version of Noah LSM) with NCAR MM5 mesoscale model (MM5) to provide the reasonable diurnal variations of surface energy fluxes and soil moisture information. The coupling system was extended to the coupled Noah/WRF model. In addition, LSM also plays an important role in the RegCM regional climate modeling system where the land surface process was described by BATS in early versions (Versions 1, 2, and 3) and a more advanced model named CLM3.5 in the latest Version 4. The LSM development promotes the incorporation of more physically based LSM to replace old ones or their old versions in the coupled mode.

In addition to the coupled mode, offline mode is also an important aspect of LSM simulation. The offline mode is standalone, and the LSMs are detached from the host atmospheric models. Driven by provided atmospheric forcing data, the offline LSMs simulations are useful to assess the realism of LSMs (evaluation, calibration, and validation), for analyzing the sensitivity to forcing parameters and land use/cover, for helping to improve the parameterizations of land surface, for providing initial (soil moisture) data for coupled runs, and for developing new methods for LSMs development. Under this mode, a big program for LSMs simulation should be referred to the Project for Intercomparison of Land Surface Parameterization Schemes (PILPS) [34]. It was conducted to perform model intercomparison to evaluate the parameterization of energy and water fluxes to/from the land-atmosphere interface and finally to improve the parameterization of the continental surface, especially hydrological, energy, momentum, and carbon exchanges with the atmosphere. The results revealed substantial differences among simulated surface energy, water vapor flux, and surface runoff, driven by identical atmospheric forcing data. Apart from this project, there are also many studies connected with LSM modification, evaluation, and land surface parameter estimation [3537].

3. Applications over Complex Terrain

In most LSMs, the subgrid composition of different land cover types is usually used to represent the spatial variability of land surface. In different LSMs, different treatments for land surface representation are applied. In Community Land Model (CLM), a three-level subgrid system is used to represent land surface with the vegetated surface specified by 15 plant functional types (PFTs) [38]. In Noah LSM, the land surface characteristics of each grid are predefined by using a vegetation lookup table for static vegetation parameters according to its land cover type in the different land cover system such as USGS, MODIFIED_IGBP_MODIS_NOAH, and USGS-RUC. However, it has only one PFT for a defined model grid. For Common Land Model (CoLM), each surface grid cell of CoLM is composed of up to five tiles: dominant vegetation, secondary vegetation, bare soil, wetland, and inland water. Each tile contains a single land cover type. Energy and water balance calculations are performed for each tile at each time step, and each tile maintains its own prognostic variables [39].

It is clear that most LSMs have tried to represent the spatial heterogeneity of the land surface characteristics, and many improvements [40, 41] have been done to better represent land surface by developing fine resolution land cover datasets or using different subgrid parameterization scheme. The changes in land surface elevation also have been considered to analyze the topographic effects on land surface thermal and hydrological processes. The Variable Infiltration Capacity (VIC) model divides the model grid cell into multiple subgrid elevation bands with 500 m elevation interval for better simulations of surface hydrological processes [42]. Ke et al. [43] developed a new subgrid classification (SGC) method based on the method of Leung and Ghan [44] to account for the effect of spatial variation of topography and vegetation cover. The model grid is divided into different elevation classes, and a variable number of vegetation types are defined for each elevation class. Therefore, the model grid is composed of a variable number of subgrid classes or land response units (LRUs). The structure is shown in Figure 3.

Despite the fact that more and more sophisticated and advanced LSMs are used to increase the representation ability of the spatially heterogenetic land surfaces, it is still difficult for the simulation results in offline modes, coupling modes, or land data assimilation systems to meet the demand of environmental, ecological, hydrological, and atmospheric research over complex terrain. The high spatial heterogeneity exists not only in land surface but also in atmospheric environment which tremendously limits the application of current LSMs. Additionally, it is also suggested from field observations that the intrinsic parameter values vary over space and environmental conditions, even for a certain PFT [45].

To address the critical issues concerning hydrological or environmental problems over complex terrain, many efforts have been made either by developing new LSMs or improving existing LSMs. Numbers of LSMs have been developed to perform mountain hydrological and ecological modelling to understand the high spatial and temporal variability of hydrometeorological conditions over complex terrain. A Regional Hydroecological Simulation System was developed by Band et al. [46] to compute and map forest evapotranspiration and net primary productivity over complex mountainous terrain. The methodology is based on the interface of geographic information processing and remote sensing with FOREST-BGC [47], a nonlinear deterministic model designed to simulate carbon and water and nitrogen cycles in a forest ecosystem. Wigmosta et al. [48] set up a physically based, distributed hydrology-vegetation model for complex terrain. The model explicitly accounts for canopy interception, evaporation, transpiration, and snow accumulation and melt, as well as runoff generation via the saturation excess mechanisms. It can be applied over a range of scales, from plot to large watershed at subdaily to daily time scales. Digital elevation data were used to model topographic controls on incoming solar radiation, air temperature, precipitation, and downslope water movement. A physical-mathematical model was also proposed to describe hourly dynamics of water and heat transfers in a forested watershed over complex terrain [49]. Interconnecting energy budget and water balance processes are integrated by simultaneously solving related equations using a specific feedback-iteration method.

One important aspect of the topographic influence owes itself to its impact on soil water drainage and soil water distribution. The nonuniform distribution of soil moisture significantly affects surface runoff, surface evapotranspiration, and surface hydrological cycle. Because most LSMs are one-dimensional model including a detailed description of the vegetation and root zone, the interactions between groundwater, root zone, and surface water, as well as the lateral surface and subsurface flows, are normally neglected [50]. The ignorance of subsurface lateral transport induced by topography or moisture gradients causes error in surface flow calculations and soil water budget which in turn influences the spatiotemporal distribution of surface, subsurface water, and terrestrial water and energy budget finally. A common way is to couple a hydrological model with LSM to represent the influence of topography on surface hydrological components and to consider the lateral surface and subsurface flows. Choi et al. [50] developed a conjunctive surface-subsurface flow (CSSF) model and coupled it with CoLM to get a hydrologically enhanced model (CoLM + CSSF). The offline evaluation results showed that the simulated interaction between surface and subsurface flow significantly improves the stream discharge prediction. Moreover, the hydrological model TOPMODEL [51, 52] was used with a high frequency because it employs digital terrain model data to simulate the lateral flow. Famiglietti et al. [53] coupled TOPMODEL with a surface energy balance model to compute the spatial variability of evaporation at the catchment scale. Band et al. [54] used a coupled model consisting of FOREST-BGC and TOPMODEL to account for the effects of lateral soil moisture redistribution on ecological processes. Similar applications of this approach can also be found in Stieglitz et al. [55], Deng and Sun [56], and Bouilloud et al. [57].

In the complex terrain hydrological modelling, snowpack is another important aspect. Snow cover has an effect on local meteorological condition, soil moisture condition, the distribution and growth season of vegetation, and the timing and availability of runoff. For this reason, many physically based snow models [58, 59] have been developed to describe snow processes and the dynamics of snow water equivalent (SWE) accumulation and ablation by solving energy and mass balance equations for the snowpack over a grid. For instance, Alpine3D developed by Lehning et al. [58] is a surface energy balance model that has been used to simulate fine scale snow processes in mountainous regions. It has been recognized as a valuable tool to investigate surface dynamics in mountains. Another model, Generate Earth Systems Science (GENESYS), was also demonstrated to have good performance in simulating daily snowpack and the spatial extent of snow cover over complex mountainous terrain [59]. Except for the snow models, the snowpack simulation with Noah LSM faces the problems of an early depletion due to excessive sublimation and untimely onset of snowmelt. To overcome the deficiencies, Barlage et al. [35] made five model modifications to Noah LSM for the simulation in the headwater region of the Colorado Rocky Mountains. The solar radiation adjustment for terrain slope and orientation is one of the modifications. The slope and orientation of the surface made substantial differences to the amount of incident shortwave radiation through change of the incident solar zenith angle. Although the results showed the terrain slope and orientation adjustment for radiation has little effect on average snow water equivalent conditions for the entire modeling domain, the addition of the terrain slope and orientation effects does have an effect on local surface energy flux components depending on the cell slope and orientation. The same evaluation is also conducted with Noah LSM in cold alpine area. Chen et al. [60] conducted simulations in the upper reaches of Heihe River at two typical sites. The comparison between simulated and observed soil temperature and liquid contents indicates that the Noah LSM is capable of monitoring the variation trend, but with some unignorable differences. The insufficient consideration of soil profile heterogeneity in Noah LSM is one of the major problems.

4. Key Issues of Land Surface Modelling over Complex Terrain

Over complex terrain, the most important phenomenon indicates that land surface characteristics such as topography, vegetation, and soils vary greatly in space combined with surface temperature, soil temperature, and moisture. Atmospheric forcing data, such as near surface air temperature, air humidity, precipitation, and wind speed, exhibit similar spatial variation characteristics. The absorption of solar radiation is highly dependent on local topography such as slope and aspect. These small scale heterogeneities impact soil moisture redistribution, runoff generation, and the land surface energy balance. It can be concluded that the complex terrain shows much greater spatial variation of land surface variables than flat areas. Therefore, the modelling work in this area is much more complex than that in the flat area.

According to the study of Duan et al. [61], several factors show significant influences on the ability of a LSM in representing the land surface processes accurately and reliably. The first factor is the equations or parameterization schemes in the model structure; the second factor is the quality of external atmospheric forcing data and the initial and boundary conditions; and the third one is the appropriateness of the model parameter specification. Based on these points, the key issues faced by current LSMs for complex terrain simulation can be concluded as below.

4.1. Horizontal Variation of Land Surface Parameters

The horizontal variation is mainly reflected by the spatial heterogeneity of land surface characteristics. Realistic and high spatial resolution representation of land surface characteristics is significant for accurate estimation of surface hydrology, heat fluxes, and surface CO2 exchanges in land surface models for applications across global, regional, and subregional scales.

Among various land surface characteristics, land use/cover is an important one which affects the biophysics, biogeochemistry, and biogeography of the Earth’s surface and thus brings about far-reaching consequences to human well-being. Due to various impacts of the big gradients of heat and water, the land cover condition of complex terrain exhibits violent spatial variations. The development of remote sensing technology also makes it easy to map global or regional land cover with a relatively high accuracy at different spatial resolution [62, 63]. Figure 4 presents the land cover map of Mount Gongga area located in southwestern China in 2010. Two land cover maps are shown in Figure 4 with one from EOS-MODIS land cover product at 500 m resolution and the other acquired from land cover classification using Landsat TM/ETM data at 30 m resolution [64]. The big vertical gradient (the highest peak of 7,556 meters) in topography induces the significant variation of land cover types in horizontal direction. Regardless of the resolution differences in these data, the derived land use/cover information from satellite data enables the LSMs to incorporate fine resolution land cover dataset to conduct land surface processes simulation at different spatial scales. However, the current LSMs are far from satisfactory. CLM has the base spatial resolution of 0.05°. The surface grid size of CoLM is 2.8°  × 2.8° [65]. The Noah LSM simulation studies were conducted at the grid size in kilometers level [66, 67]. Although tiles or subgrids are used to present the heterogeneity within one grid cell, the coarse grid cell sizes are not suitable to resolve the spatial heterogeneity of land surface processes over complex terrain at local scale.

Additionally, as an important parameter in LSM parameter components, soil property information determines the physical and chemical properties of soil, including the representation of the soil column and the infiltration parameterization. It is particularly important due to its impacts on the vertical distribution (evaporation, transpiration, and infiltration), simulation of soil moisture, and the influence on the energy balance by controlling the parameter albedo. The sensitivity study [68] shows that hydraulic conductivity at saturation, parameter for soil water potential calculation in Clapp and Hornberger [69], and wilting point have a profound impact on the simulation of soil moisture. In current LSMs, soil physical parameters are usually predefined for each soil texture type. For example, four parameters (soil wetness exponent, soil tension at saturation, hydraulic conductivity at saturation, and soil porosity) are defined for seven types in SiB2 [70]. Noah LSM also provides soil parameters associated with 19 soil types. Although the soil parameterization has predefined the parameter values, soil texture type distribution varies greatly and usually has big change even within a short distance [71]. For complex terrain area, the elevation change and the aspect difference influence the development of soil and finally affect the soil texture type spatial distribution. This phenomenon is even more serious in high mountain area [72]. It is difficult for current soil texture distribution dataset applied in CoLM or other LSMs to describe the high spatial variation over complex terrain. Therefore, the accurate description of soil property horizontal variation is one important issue that should be considered in land surface modelling in complex terrain area.

4.2. Vertical Variation of Atmospheric Forcing Data

In vertical direction related to elevation change, as a vital issue, the variation of atmospheric forcing data should highly vary with the impact of the significant topographic changes over complex terrain. Among the forcing data, air temperature is one of the key factors driving LSMs. It is essential for land surface energy flux and water flux modelling to strictly control the energy exchange between land surface and near surface atmosphere. Generally, air temperature decreases linearly with the increase of elevation. However, for current LSMs, the air temperature in the gridded analysis is considered to be homogeneous in each coarse grid. The neglect of the impact of the variation of the inner elevation in each grid will definitely bring a big uncertainty to the energy fluxes simulation results in complex terrain area when the grid has a very coarse spatial resolution.

Besides air temperature, the other forcing data such as solar radiation, wind speed, precipitation, and humidity also face the same problems. The high spatial variation of air temperature and precipitation observations in the Langtang valley in the Nepalese Himalayas also certified the large influence of elevation or topography on spatial change of atmospheric parameters [73]. As shown in Figure 5, seven sites are located at different elevations in a small watershed, and the observed mean air temperature of different time periods for these sites exhibits great vertical variation with elevation change.

Therefore, considering the impact of topography, the applications of LSM at different sites over complex terrain are usually driven by the forcing data from their corresponding meteorological stations in previous simulation study [74, 75]. Dense networks of meteorological sites with high temporal resolution are required to characterize the patterns of atmospheric forcing data that occur over complex terrain. However, such dense field observations are seldom available especially for rugged regions. The typically sparse distribution of meteorological sites inadequately resolves the spatial and temporal pattern of atmospheric forcing data. It is not realistic to complete the simulation over complex terrain with site observations to establish the accurate description of surface atmospheric environment. Additionally, the comparison study [76] for simulations with process-based hydrologic model forced by two different types of forcing data illustrated that the differences (such as shortwave radiation) greatly impacted hydrologic states and fluxes, particularly at high elevation. The differences altered the timing of snowmelt and runoff and the partitioning of precipitation between evapotranspiration and runoff. Therefore, the studies directly proved the importance of the accuracy of atmospheric forcing data in land surface modelling over complex terrain. It is a big challenge to generate reliable atmospheric forcing data to capture the spatial heterogeneity of atmospheric environment.

4.3. Scale of Land Surface Parameterization

Land surface parameterization is a vital step to conduct land surface modelling. For large-scale application, the surface data are provided with coarse spatial resolution. For example, the surface dataset required by CLM, including percent glacier, percent lake, percent wetland, percent urban, percent sand/clay, soil organic matter density, soil color, and PFTs, is parameterized with the grid cell size changing from 5 minutes to 1 degree. However, for CoLM, the surface dataset, including terrain elevation, land-water mask, 25-category land cover, and 17-category soil, is parameterized with the grid cell size of 30 seconds. The coarse resolution enables most LSMs to couple with atmospheric model for climate simulation at continental or global scale [7779]. Although the above parameterization of land surface is enough for the global studies, the large-scale simulation is not useful to analyze the impact of topography on the spatial variation of the radiation environment, precipitation, temperature, and soil water drainage. Consequently, the application of current LSMs with coarse land surface parameterization dataset cannot effectively monitor the spatial distribution of soil-vegetation-atmosphere interactions over complex terrain at local scale. The application of LSM for complex terrain land surface modelling will cause spurious drifts in simulation results because the LSM parameters are of high uncertainty and of high sensitivity to the simulated states and fluxes. Although the application of LSMs over complex terrain has been performed for point-based simulation, regional simulation is hardly conducted. Therefore, the parameterization of land surface characteristics with high spatial resolution should be developed to present the land surface processes over complex terrain.

Furthermore, it is important to realize that the increase in spatial resolution implies increased complexity but not necessarily a higher degree of accuracy in simulation [80]. As indicated in Figure 4, two land cover maps at different spatial scales are shown for the same area. It is obvious that the high resolution map (30 m) from Landsat data can present more detailed information (not only pixel size but also the number of land cover types) about the spatial distribution of land cover types in the mountainous area than the MODIS land cover map (500 m). However, it is still uncertain which scale associated with the land cover map is more suitable to perform the land surface parameterization and conduct land surface processes simulation. In general, increased model resolution is able to depict land surface characteristics with a relatively high accuracy. In practice, this is not always true because of the limitations in physical parameterizations and input data, as well as other numerical errors that can be introduced in model operation. In addition, a very fine spatial resolution may violate the assumptions of some physical parameterizations or numerical solutions. Wood et al. [81] discussed the challenges for monitoring terrestrial water with hyperresolution global land surface modelling. The study showed that the hyperresolution land surface modelling will be particularly challenging because the basis for the current estimation of the exchange of mass, energy, and momentum from the land surface to the atmosphere will not be suitable at hyperresolution. Although high resolution simulation is useful to capture the spatial patterns of land-atmosphere interactions at local scale, the improvement in the accuracy of the simulation may be relatively little. Therefore, the determination of an appropriate spatial scale in the simulation is the primary and important work for accurately modelling land surface processes for complex terrain area.

4.4. Two-Dimensional Modelling

Soil moisture is a link between the surface energy, water, and biogeochemical cycles. It influences the partitioning of precipitation into infiltration and runoff, controls the biogeochemical processes and surface evapotranspiration by limiting water availability to plants, and affects the partitioning of energy into latent and sensible heat. Over complex terrain, soil water distribution is mainly controlled by vertical and horizontal water divergence and convergence, infiltration recharge, and evapotranspiration. The topography plays a key role in water resources forming and runoff generation. Generally, regions of local concavity (lowlands) tend to be zones of convergent flow (surface and baseflow) and therefore zones of high soil moisture content. In comparison, uplands soils tend to be progressively drier [55].

However, most LSMs usually use one-dimensional mass conservation equations based on Darcy’s law to simulate the hydrological components (soil water and surface runoff). Soil water is simulated from a multilayer model, in which the vertical soil moisture transport is governed by infiltration, surface and subsurface runoff, gradient diffusion, gravity, canopy transpiration through root extraction, and interactions with groundwater. Taking CLM as an example, the water flow in soil layers can be stated by the conservation of mass equation:where is the volumetric soil water content; t is time; is the height above some datum in the soil column; is soil water flux; and is a soil moisture sink term. Only soil layers existing in the vertical direction are used to model the water transfer. The horizontal interactions of land surface processes between adjacent grids such as the lateral flow redistributing water via both surface and subsurface flow pathways are commonly ignored. Each grid is an independent soil-vegetation-snow-atmosphere system, and they are not affected by what happens in systems of the neighboring grid squares. The major reason can be attributed to the relatively coarse spatial resolution of LSMs used to couple with global or regional climate model that does not allow accurate representation of surface and subsurface water movement.

Therefore, the use of current one-dimensional LSMs over complex terrain is unable to describe the horizontal movement of soil water within soil layers which is seriously impacted by topographic factors. Great uncertainty is brought to the prediction of soil moisture distribution and the generation of surface runoff.

5. How to Model Land Surface Processes over Complex Terrain

Currently, land surface processes simulation operated by LSMs has made much progress. Much work has been done to address the issues related to surface energy, water, and carbon cycles over complex terrain.

5.1. Generating Atmospheric Forcing Data at Fine Resolution Scale

For the land surface boundary information, atmospheric forcing data are one of the key inputs for LSM simulation. Theoretically, atmospheric forcing data determine the external atmospheric environment for land-atmosphere interaction and have strong spatial heterogeneity over complex terrain. Therefore, high accuracy atmospheric forcing data at fine resolution scale are very important driving factors for land surface processes simulation with LSMs over complex terrain.

In the forcing data, the incident radiation controls the radiation received by land surface and plays a major role in LSMs’ surface energy balance. Land surface topography directly alters incident shortwave radiation due to the shading effect and multiple terrain reflections [82]. It then influences the downwelling longwave radiation because of the effects from the longwave irradiance emitted from the surrounding terrain [83]. In addition, slope irradiance also exerts influence on incident radiation. Usually, it is neglected when the spatial resolution is very coarse (tens of kilometers) or the slope is relatively low or moderate. However, the contribution of slope irradiance to surface radiation will be enlarged when the spatial resolution is higher (1 km or less) with steep slope.

Taking the shortwave radiation as an example, Figure 6 illustrates the radiation accepted by land surfaces over complex terrain. For sunlit area (), surface absorbed radiation includes direct solar radiation (), sky radiation (), and radiation reflected by adjacent terrain (). However, for shaded area (), solar irradiation is totally shaded. Therefore, the topographic influence mechanism should be imported into LSMs to discriminate the spatial variation in incident radiation. The geometrical factors such as slope and aspect should be emphasized in the study. Hauge and Hole [84] proposed a surface radiation flux estimation method by taking into account the effect of slope and aspect. The total irradiance on a surface inclined by an angle towards an azimuth angle can be expressed aswhere is solar elevation; is ground slope (calculated using forward differences); and is the solar beam angle of incidence. is the diffuse sky irradiance; is the direct radiation (beam); and is ground reflected irradiance. Helbig and Löwe [82] also proposed a complete radiation parameterization scheme for subgrid topography accounting for shading, limited sky view, and terrain reflections.

In addition to the incident radiation, the influence of topography on the other atmospheric parameters (air temperature, wind speed, humidity, precipitation, and air pressure) should also be seriously considered. The homogeneous distribution treatment in the global simulation is not appropriate for the land surface simulation over complex terrain. With the differences in altitude, latitude, and terrain environment, the atmospheric parameters show significant differences between the valley area and the top of the mountain. However, because there are insufficient meteorological stations located over complex terrain especially in the high mountain area, the atmospheric parameters could not be provided accurately at local scale.

The traditional way to get these small scale data is to apply the spatial interpolation method by using sparse site measurements combined with geolocation and topographic information [85, 86]. The interpolation methods include the Inverse Distance Weighting (IDW) approach, the spline method, the kriging method, regression analysis, and geostatistical methods. Therefore, they are used to generate fine resolution atmospheric forcing data to enhance the simulation over complex terrain. During the interpolation, to account for the physiographical controls on spatial and temporal distribution of atmospheric forcing data, numerous methods have been proposed to parameterize the relationship. For example, precipitation-elevation relationships were built to extrapolate precipitation from station observations to spatial distribution [59]. For air temperature, temperature lapse rates were used to account for temperature changes as a function of elevation. Accordingly, high resolution gridding of atmospheric data often relies on assumptions such as a constant surface temperature lapse rate of 6.5°C km−1 or 6.0°C km−1 [87, 88]. However, Minder et al. [89] found that it should not be a uniform or constant value in terms of the impacts from air mass characteristics (e.g., temperature, humidity), solar radiation, and other factors. Immerzeel et al. [73] analyzed the air temperature lapse rates and precipitation gradients in the Langtang valley. The results indicated that the lapse rates were seasonal differences connected with water vapor in the atmosphere, and the precipitation gradients showed big horizontal variation.

Another source of generating atmospheric forcing data is from the simulation results of climate model. However, as climate models are still computationally limited to relatively coarse spatial resolutions of  km at its best, the high resolution data are obtained by applying downscaling method to disaggregate low spatial resolution atmospheric parameters derived from climate model. A quasi-physically based, high resolution meteorological distribution model (MicroMet) was developed by Liston and Elder [90] to produce high resolution meteorological forcing distributions. The relationships between atmospheric parameters and the surrounding landscape (primarily topography) are introduced in the downscaling procedure. Zabel et al. [91] used another scaling tool SCALMET (Scaling Meteorological variables) developed by Marke et al. [92] to downscale the coarse meteorological data provided by MM5 climate model to higher resolution for the surface fluxes simulation with PROMET. This tool has proven to be robust under a variety of climatological and hydrological conditions and has already been successfully applied in various regions of the world. The downscaling methods provide a new direction for overcoming the limitation of field observations to run spatially distributed LSMs over complex terrain.

Therefore, based on the previous studies, it is expected that we have more knowledge about the topographic control on atmospheric parameters variations and the mechanism behind the differences over complex terrain and that we know more about how to develop a systematic approach to generate the reliable forcing data.

5.2. Solving Two-Dimensional Movement of Soil Water

It is clear that the parameterizations of spatial heterogeneity in land surface processes should accurately evaluate the distributions of water movements in both vertical and horizontal directions that subsequently affect the water and energy exchange between land surface and atmosphere. As the issue indicated in Section 4.4, most LSMs focus on the water movement in the vertical direction. Multilevel soil layers, vegetation layers, and snow layers are usually separated to describe the processes. Under this parameterization, the effects of heterogeneous distribution of soil moisture on surface runoff generation, subsurface flow, and surface evapotranspiration are totally ignored. With this deficiency, it will obviously influence the surface energy partition, the relationship among precipitation, runoff and evapotranspiration, and the calculation of surface runoff and underground flow. These effects will cause more serious problems with the function of gravity over complex terrain. Therefore, it demands LSMs to be developed to describe the soil water and hydrological processes in two dimensions.

To solve the downslope water movement over complex terrain, many hydrological models, such as the Distributed Hydrology-Soil-Vegetation Model (DHSVM) of Wigmosta et al. [48], its later version the GIS-based modeling system for watershed analysis (GISWA), Variable Infiltration Capacity (VIC) model [93], and TOPMODEL [51], have been developed. In general, topography influences land surface hydrological processes such as precipitation and runoff and controls the spatial distribution of air temperature, precipitation, soil, and vegetation within a catchment.

Among these hydrological models, topography based TOPMODEL has the biggest advantage as concluded by Deng and Sun [56]: TOPMODEL is built on the clear physical basis. The solutions of hydrological processes in TOPMODEL are simple analytical ones with few parameters, and (3) TOPMODEL is able to reflect the effects of topography on soil moisture and groundwater movements. The subsurface flow treatment in TOPMODEL attracts many attentions from researchers of LSMs study, and many studies have been conducted by coupling it with LSMs [56, 57]. It suggests a general direction for the model development by utilizing the advantage of hydrological model in the consideration of soil water spatial distribution and topography effect.

In addition to coupling with hydrological model, LSM modification is another direction. Wang and Chen [94] imported Muskingum method into Noah LSM to improve the simulation accuracy of the surface runoff by combining flow concentration values for considering the redistribution of soil moisture in a two-dimensional level. Results indicate that the soil moisture and temperature simulated by the improved model are closer to the observed values. Compared with the division of land surface into a number of regular grids, Tesfa et al. [95] suggested a novel way for land surface treatment by dividing the surface into irregular subbasins. This method has the following advantages. First, the subbasin based representation follows the natural topographic division and river network structure. Second, the redistribution of water from neighboring subbasins should not be considered. Third, the representation is more suitable for coupling hydrologic component in LSM and agrees with the conceptual basis of most hydrological schemes.

In general, the studies concerning coupling with hydrological models or modifying the hydrologic component in LSM to resolve the two-dimensional movement of soil water in land surface modelling indicate a general direction for the simulation of land surface processes over complex terrain.

5.3. Coupling LSM with Remote Sensing Observations

In recent years, the development of remote sensing techniques accelerates the development of the quantitative remote sensing methods. Numerous land surface parameters are able to be provided by remote sensing methods at different spatial and temporal scales. It offers a big chance to improve LSM simulation over complex terrain.

5.3.1. Dynamic Information from Remote Sensing

One important aspect of coupling LSM with remote sensing observations is to use the land surface dynamic information derived from remote sensing data. The initial state of land surface is an important input of LSM simulation. For complex terrain, remote sensing techniques can provide time series of observations for land surface parameters such as land cover, fraction of vegetation cover (FVC), leaf area index (LAI), and surface albedo with various spatial and temporal resolution. It can get rid of the difficulty to parameterize surface condition over complex terrain due to the lack of field observations in this area. The dynamic surface information guarantees that the LSM initialization obeys the right condition at the corresponding resolution. Yuan et al. [96] used the vegetation information from AVHRR data to simulate streamflow for Hanjiang River basin with VIC-3L land surface model. Choi [97] assessed the applicability of high-resolution simulations for terrestrial processes in a small study basin using CoLM in 1 km surface boundary conditions (SBCs) with remote sensing products from MODIS and SPOR-VGT observations. Ghilain et al. [98] proposed a practical methodology using the biophysical products LSA-SAF LAI and FVC from MSG SEVIRI data in H-TESSEL land surface models. Compared to the use of the semistatic database ECOCLIMAP-I, the proposed method showed better ability in monitoring surface evapotranspiration.

The applications suggest that the remote sensing method has become an important way for LSM initialization and calibration, and frequent observation of land surface parameters by remote sensing can specify the dynamic change of surface conditions in land surface modelling over complex terrain. However, with the increase of remote sensing observations, multiple observations will participate in LSM simulation. The application of multiple observational constraints should be concerned about. Model-data fusion will be one option to improve the observational constraint on the model over that afforded by any single data set on its own [99]. Renzullo et al. [100] demonstrate an application of a “multiple constraints” model-data fusion (MCMDF) scheme to integrate AMSR-E soil moisture product and MODIS LST data with a coupled biophysical model of surface moisture and energy budgets for savannas of northern Australia. The result suggested that the approach is able to improve the predictive capability for a range of surface variables and fluxes.

5.3.2. Model Evaluation

The evaluation of LSM simulation is another aspect for coupling LSM with remote sensing. Evaluation of the performance of LSMs is a complex problem, yet it is necessary to determine the efficacy of hydroclimatological predictions. For streamflow validation, the traditional methods are to compare the predicted and observed streamflow over an entire watershed level [101]. For energy and water flux, field measurements from field stations are usually used to validate the simulated fluxes under the assumption of homogeneous distribution [102]. It is obvious that these validation methods face difficulty in determining whether the model can describe the spatial variation of energy and water fluxes well due to the effect from the inconsistency in spatial scale. The evaluation work is constrained by the number of field stations; and the spatial representativeness of the field measurements is a big issue especially for complex terrain.

Compared with the above validation methods, remote sensing techniques offer an alternative data source for LSM evaluation. Currently, the quantitative estimation of land surface parameters, such as land surface temperature (LST), land surface emissivity (LSE), evapotranspiration (ET), soil moisture, LAI, FVC, net primary productivity (NPP), and energy fluxes, has made huge improvements in accuracy by using remote sensing data. The estimation results have been applied a lot to evaluate the performance of LSM simulation. Rhoads et al. [103] compared surface temperature simulated by the Variable Infiltration Capacity (VIC) hydrologic model with surface temperature retrievals from TOVS and GOES over the Arkansas-Red River basin. The results showed that modeled and satellite-derived surface temperatures agree well when aggregated in space or time. In order to evaluate the performance of the Noah LSM after some parameterization, Wei et al. [104] used satellite-retrieved land surface skin temperature combined with site measurements to validate the simulation results and found that there are significant reductions in biases in latent heat, total runoff, and land skin temperature simulation.

Although the evaluations quantitatively showed the performance of LSM, the comparisons often focus just on one aspect of the LSM. There are more and more recognitions from land modelling community to use a benchmarking system to make more comprehensive model evaluation. Many studies have been done [105108], and the remote sensing dataset becomes a key component in the benchmarking system. Therefore, to effectively evaluate the LSMs simulation over complex terrain, it is important to build a benchmarking system based on a collection of remote sensing datasets and site-based observations for key land surface processes evaluation.

5.3.3. Data Assimilation

LSM driven by atmospheric forcing data has shown great potential in generating estimates of land surface states at different temporal resolution from a regional to global scale. However, such estimates will be heavily dependent on the characteristics of the chosen LSM and the quality of the atmospheric data used. In addition, because LSMs are simplifications of land surface processes, there are unavoidable model errors and uncertainties in the initialization and running process. Even though they are driven by identical atmospheric forcing data, different LSMs may give very different estimates of surface fluxes [109]. Therefore, it is desirable to combine model simulations with remote sensing observations to update land surface simulation states so as to minimize the effects of errors in the model physics, model parameters, and atmospheric forcing.

Data assimilation is an alternative to integrate the data from variety of sources with different resolutions and accuracies with model prediction to improve deterministic model accuracy [110]. The popular data assimilation methods are Variational Assimilation (VA), the Ensemble Kalman filter (EnKF), and Particle Filter (PF).

In recent years, the applications of data assimilation in land surface modeling have expanded with the increased availability of remotely sensed land surface variables. Under the support of satellite and ground based observations, a Global Land Data Assimilation System (GLDAS) has been developed by using a series of advanced LSMs and data assimilation techniques to generate optimal fields of land surface states and fluxes at global and continental scale [111]. GLDAS is unique in that it is an uncoupled land surface modeling system that drives multiple models, integrates a huge quantity of observation-based data, and produces the hourly 1/8th degree product for North America (NLDAS) or the 3-hourly 0.25 and 1.0 degree products globally. Similar systems also exist in Europe (ELDAS) and China (CLDAS). Except for the large scale, site specific or watershed scale applications have been conducted by many researchers [112, 113]. These studies are performed by using different data assimilation algorithms coupled with different LSMs with various satellite based observations. Meanwhile, the challenges in current land surface data assimilation are also concluded from previous studies. Because many data assimilation techniques assume bias-free models as well as observations, it is important to deal with the biases of model predictions or remote sensing observations, the characterization of second order error, and the nonclosure or imbalance problems [114]. In general, the studies help to establish the theoretical basis of similar applications over complex terrain. Data assimilation is becoming a promising way to overcome the limits in the uncertainties of land surface processes simulations in this special area.

In addition to the above aspects addressing the key issues related to land surface modeling over complex terrain, there are also some issues related to LSM development that need to be considered in further study: the importance of diffusional limitations on carbon fluxes (e.g., stomatal and mesophyll conductance) and the representation of land use change or dynamic vegetation. As pointed out by Claussen et al. [115], the changes in land cover and land use have a huge impact on Earth’s systems including the biogeophysical and biogeochemical effects on the atmosphere. The vegetation dynamics are greatly correlated with the carbon assimilation and allocation. The importance analysis of vegetation dynamics for future terrestrial carbon cycling by Ahlström et al. [116] showed that the vegetation dynamics represent a significant fraction on the uncertainties in modelled global C-uptake. It is realized that increasing complexity in LSM may increase realism of LSM simulation result, but the reliability and robustness of LSM will decrease because of the addition of poorly known model parameters. Therefore, against this fact, Prentice et al. [117] suggested simplifying complex processes by stochastic parameterisation (the representation of unresolved processes by statistical distributions of values) for next-generation complex LSMs.

6. Conclusion

The scientific importance of complex terrain has been widely recognized by many researchers in recent years. Complex terrain plays a critical role in global water cycle and carbon cycle, and its ecosystem is also globally important as the center of biological diversity. However, the big spatial variation in surface parameters occurring with sparse field observations constrains the understanding of land surface processes and their dynamic changes in this area. Therefore, it is urgent to develop an effective way to monitor and capture the dynamic information for global change study. Although field observations are able to measure land surface parameters with a high accuracy, it is inappropriate to deploy sufficient sites over complex terrain. The physical-based LSM provides a valuable tool for mountain ecological and environmental research and sharpens our understanding of the land surface processes over complex terrain.

The development of LSM and applications over complex terrain are reviewed. Due to the impact of topography change, the atmospheric environment, surface energy partition, and energy/water/carbon cycle in this area have significant difference compared with those of flat area. Key issues faced by LSM application over complex terrain are concluded as horizontal variation of land surface parameters, vertical variation of atmospheric forcing data, scale problem of land surface parameterization, and two-dimensional modelling of soil water movement.

To resolve these issues, researches are needed to determine the appropriate scale which is suitable for LSM simulation over complex terrain. The simulation at this scale should be able to present the spatial heterogeneity of land surface processes within the accuracy requirement. Then, efforts should be made to generate atmospheric forcing data at this fine resolution scale by considering the topographic impacts on key forcing data. Moreover, the two-dimensional movement of soil water should be solved with the help of coupling hydrological model or LSM modification. It is also important to couple LSM with remote sensing observations by integrating remote sensing observations into LSM initialization, simulation, and evaluation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This study is partially supported by International Cooperation Partner Program of Innovative Team (no. KZZD-EW-TZ-06) of Chinese Academy of Sciences, the National Natural Science Foundation of China (nos. 41271433 and 41401425), the External Cooperation Program of BIC of Chinese Academy of Sciences (GJHZ201320), the West Light Foundation of Chinese Academy of Sciences, and the “Hundred Talents” Project.