• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Evaluation of the New Dynamic Global Vegetation Model in CAS-ESM

    2018-04-08 10:59:07JiawenZHUXiaodongZENGinghuaZHANGYongjiuDAIDuoyingJIFangLIQianZHANGHeZHANGandXiangSONG
    Advances in Atmospheric Sciences 2018年6期

    Jiawen ZHU,Xiaodong ZENG,M inghua ZHANG,Yongjiu DAI,Duoying JI,Fang LI,Qian ZHANG,He ZHANG,and Xiang SONG

    1International Center for Climate and Environment Sciences,Institute of Atmospheric Physics,Chinese Academy of Sciences,Beijing 100029,China

    2Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters,

    Nanjing University of Information Science and Technology,Nanjing 210044,China

    3School of Marine and Atmospheric Sciences,Stony Brook University,NY 11790,USA

    4School of Atmospheric Sciences,Sun Yat-Sen University,Guangzhou 510275,China

    5College of Global Change and Earth System Science,Beijing Normal University,Beijing 100875,China

    1. Introduction

    Land vegetation plays a pivotal role in regulating the exchange of heat,water,and carbon fluxes between the land surface and the atmosphere(Li and Xue,2010;Xue et al.,2010).The seasonal grow th of vegetation significantly impacts surface latent heat,downward and reflected solar radiation(Zhu and Zeng,2015,2017),and the interannual variability of vegetation can significantly regulate surface energy budgets via evapotranspiration(Guillevic et al.,2002;Zhu and Zeng,2016).Vegetation dynamics are considered to be as important for climate as atmospheric dynamics,ocean circulation(Pielke et al.,1998),and have gained much attention in recent years(Cramer et al.,2001;Seddon et al.,2016).

    In the past several decades,rapid climate change has resulted in considerable changes in terrestrial ecosystems.One w idespread dynam ic change is the so-called“greening”of arctic ecosystems—a feature considered to be majorly contributed by boreal shrubs(Fraser et al.,2011;M yers-Sm ith et al.,2011).Both observation-and model-based studies have shown an increase in shrub biomass,coverage and abundance in arctic ecosystems in recent decades(Myers-Smith et al.,2011).Another key region is the tropics,which is ecologically sensitive to climate variability and shows ampli fi ed responses compared to other regions(Seddon et al.,2016).Many tropical ecosystems are more sensitive than other regions to environmentalperturbations and externaldisturbances,and are highly likely to cross a threshold to an alternative state(Holling,1973;Scheffer et al.,2009).

    These vegetation changes can in turn significantly affect the climate.Over high latitudes,higher shrub abundance warms the winter soil temperature by trapping more snow via their branches(Sturm et al.,2001).During early spring,taller shrubs that expand into arctic tundra ecosystems also tend to systematically warm the soil because of their lower albedo than snow(Bonfils et al.,2012).In boreal summer,however,the shading imposed by shrubs is known to reduce the soil temperature and consequently the thickness of the active layer(Blok et al.,2010).Many studies have also reported that shrub canopies can alter arctic nutrient cycling,biodiversity and ecosystem services(Myers-Smith et al.,2011).Conversely,in the tropics,the complete deforestation of Amazonia may result in warmer,drier conditions at the local scale and lead to extratropical changes in temperature and rainfall through teleconnections(Lawrence and Vandecar,2015).Reduced plant cover in the Sahara,meanwhile,may cause a decline in precipitation because of the resultant increase in albedo(Charney et al.,1975),and the vegetation dynamics of West Africa have been shown to play a crucial role in regulating present-day climate,and probably future climate,via influenceson precipitation,evapotranspiration and energy balance(Erfanian et al.,2016;Yu et al.,2016).

    Atpresent,Dynamic Global Vegetation Models(DGVMs)are the best available tool to represent vegetation dynamics at the global scale(Quillet et al.,2010).They can simulate transient and potential responses of vegetation to past and future climate change via parameterizing physical and biogeochemical processes of vegetation(Peng,2000).With the help of DGVM s,it is possible for global climate models(GCM s)to include the bidirectional interactions between vegetation and climate(Quillet et al.,2010).Therefore,the coupling of DGVMs and GCMs is an important approach to assess the influences of climate change on vegetation dynamics and their feedbacks to climate change.

    However,there are still many uncertainties in reproducing and forecasting vegetation dynamics with DGVMs because of the complicated physical and biogeochemical processes.In response to climate change,most DGVMs simulate an increase in woody coverage over high latitudes,while others suggest a gain in herbaceous vegetation or no changes(Falloon et al.,2012).In terms of the Amazon ecosystem,the predictions from different models also vary widely,which is dominated by the differences in large-scale forest dieback and forest resilience(Betts et al.,2004;Friedlingstein et al.,2006;Baker et al.,2008;Restrepo-Coupe et al.,2017).Falloon et al.(2012)suggested that these responses of DGVM s to climate change are strongly associated with their simulation of present-day vegetation cover.Therefore,an important step to reduce these uncertainties in DGVMs is to systematically evaluate and understand their present-day performance,which is the main focus of the present study.

    CoLM-DGVM is the default DGVM of the Common Land Model(CoLM).It was developed from the Lund–Potsdam–Jena DGVM(LPJ-DGVM;Sitch et al.,2003),and combined an early version of a temperate shrub sub-model of the DGVM developed at the Institute of Atmospheric Physics(IAP-DGVM;Zeng et al.,2008).In recent years,a boreal shrub sub-model,a process-based fire parameterization and a new establishment parameterization scheme have been further developed in IAP-DGVM.These developments have resulted in considerably improved reproductions of the presentday vegetation distribution and carbon cycle by IAP-DGVM(Zeng et al.,2014).At present,IAP-DGVM has been coupled to CoLM,and both are important components of the current version of the Chinese Academy of Sciences’Earth System Model(CAS-ESM).This work assesses the performance of IAP-DGVM,through comparison with that of CoLMDGVM,within the framework of CAS-ESM,which is a necessary step to using CAS-ESM for investigating vegetation–climate interactions.In the next section,the model,experimental design and observational data are described.The results are presented and discussed in section 3,followed by conclusions in section 4.

    2. M odel description,experimental design and observational data

    2.1. Model description

    IAP-DGVM is a DGVM developed at the Institute of Atmospheric Physics,Chinese Academy of Sciences.It originates from LPJ-DGVM(Sitch et al.,2003)and the Community Land Model’s DGVM(CLM-DGVM;Levis et al.,2004).Plants in IAP-DGVM are classified into 14 plant functional types(PFTs)(Table 1),of which eight are trees,three are shrubs and three are grasses.These PFTs are defined by their physical,phylogenetic and phenological parameters,and are assigned bioclimatic limits to determine their establishment and survival.At present,crops are not simulated in IAP-DGVM.More details about the model can be found in Zeng et al.(2014).

    In recent years,IAP-DGVM has undergone several major developments.These include the follow ing:(1)A shrub sub-model was established(Zeng et al.,2008;Zeng,2010).With this sub-model,IAP-DGVM can reproduce the global distribution of temperate and boreal shrubs realistically,and distinguish shrubs from grasses effectively(Zeng et al.,2008;Zeng,2010).(2)A process-based fire parameterization ofintermediate complexity was introduced(Li et al.,2012),which comprises fire occurrence, fire spread and fire impact.This fire parameterization significantly improves simulations of global fire,including burned area and fire carbon emissions(Li et al.,2013).The fire parameterization is now also adopted in the Community Earth System Model at the National Center for Atmospheric Research(NCAR)to investigate the influences of fire on carbon balance in terrestrial ecosystems,and on global land energy and water budgets(Li et al.,2014;Li et al.,2017;Li et al.,2017).(3)A new establishment parameterization scheme was developed(Song et al.,2016).This scheme significantly improves IAP-DGVM’s simulation of vegetation density by introducing soil water as an impact factor.These improvements,together with other optimized modifi cations,contribute to a better performance of IAP-DGVM in reproducing thepresent-day vegetation distribution and carbon cycle.

    Table 1. Plant functional types(PFTs)in IAP-DGVM and their corresponding abbreviations in this paper.

    The land surface model used in this study is CoLM.Starting from the code of the NCAR’s Land Surface Model(Bonan,1996),Dai et al.(2003)developed the first version of CoLM,combining the codes of the Biosphere–Atmosphere Transfer Scheme(Dickinson et al.,1993)and the IAP’s land model(Dai and Zeng,1997).Since then,CoLM has been continually improved at Beijing Normal University in many aspects and has been adopted as the land component of the Beijing Normal University Earth System Model(Ji et al.,2014).

    2.2. Experimental design

    Two global simulations were conducted within the framework of CAS-ESM.One,which coupled CoLM and CoLMDGVM,is the control(CTL)experiment,while the other,which coupled CoLM and IAP-DGVM,we refer to as IAP.Both simulations were spun up from bare ground for in excess of 1200 years.Then,a further 660 years were run to approach an equilibrium state through cycling the 33-year(1972–2004)atmospheric forcing data of Qian et al.(2006),with a T85 resolution(128×256 grid cells).The relative humidity and associated fire parameters needed for the IAP simulation were derived from Li et al.(2012)and fi xed in 2004.The atmospheric CO2was set to 365 ppmv over all simulated years.

    2.3. Observational data

    This paper focuses on the vegetation distribution,carbon cycle and leaf area index(LAI).The observed vegetation distribution and LAI came from CLM 4 surface data,which were themselves derived from Moderate Resolution Imaging Spectroradiometer(MODIS)measurements(Law rence and Chase,2007).The benchmarks for gross primary productivity(GPP),net primary production(NPP)and fire carbon emissions were from Beer et al.(2010),MODIS(Zhao and Running,2010)and version 4 of the Global Fire emissions Database(GFEDv4;Randerson et al.,2015),respectively.The average of the last five years(2000–2004)of the IAP simulation was compared with that of CTL and these benchmarks.Furthermore,to reduce the impacts of crops,the vegetation cover in each grid cell was weighted by a factor of(100%-FCcrop),where FCcrop is the fraction of crop coverage in the CLM 4 surface dataset(Zeng et al.,2014).

    3. Results and discussion

    3.1. Vegetation distribution

    Fig.1. Zonal average fractional coverage(units:%)of(a)trees,(b)shrubs,(c)grasses and(d)bare groundin CTL(blue),IAP(red)and CLM 4 surface data(OBS;black).

    In general,IAP simulated more realistic distributions for the four aggregated vegetation types(trees,shrubs,grasses and bare soil)than CTL.Over most latitudes,trees simulated by IAP were in better agreement with the observation,relative to CTL,which produced more trees(Fig.1a).In CTL,trees were overestimated over northern high latitudes,southeastern South America and A frica,with magnitudes of 50%(Fig.2a).Moreover,there was a band over central Eurasia where CTL’s trees were underestimated by more than 50%.In contrast,IAP simulated fewer trees over the tropics and more over central Eurasia than CTL,which resulted in a reduction in IAP’s biases(Fig.2b).Further investigation indicated that the new establishment scheme contributed most to the reduced biases of tropical trees,especially over transition zones(Song et al.,2016).

    In terms of shrubs,both CTL and IAP underestimated them over arctic regions,such as northeastern Canada and the northern coastline of Eurasia(Figs.1b,2c and d).Further sensitivity experiments and analysis suggested that these underestimated shrubs were lim ited by the m inimum threshold of grow ing degree days over 5°C(GDD5)set by the model,which is 350.Shrubs were unable to establish because the annual GDD5 was smaller than 350.Over northern high latitudes,shrubs could not grow and were severely underestimated in the CTL simulation,while IAP simulated a more sim ilar shrub pattern than CTL with that observed(Fig.1b).This improvement can be attributed to the boreal shrub submodel of IAP-DGVM,which can distinguish shrubs from grasses effectively(Zeng et al.,2008;Zeng,2010).However,IAP simulated more shrubs than CTL over northern middle latitudes,such as western North America and central Eurasia,which further increases the biases of CAS-ESM.Moreover,both CTL and IAP underestimated the shrub coverage in the Southern Hem isphere,such as in Australia.

    Fig.2. Differences in fractional coverage(units:%)of(a)trees,(c)shrubs,(e)grasses and(g)bare ground between CTL and observations.(b,d,f,h)as in(a,c,e,g),respectively,but between IAP and observations.

    IAP’s grasses also agreed better with the observation than those of CTL(Fig.1c).The grasses in CTL were largely overestimated over middle and high latitudes,such as northeastern Canada,central North America,middle and eastern Russia and central Eurasia,where the biases exceeded 50%.However,IAP reduced these deficiencies to around 10%over these regions(Fig.2f).Over the tropics,both CTL and IAP underestimated grasses,although IAP’s biases were a little smaller than those of CTL.The main underestimation was in southeastern South America and most regions of Africa.

    Figure 1d shows that the bare soil simulated by CTL was underestimated over middle and high latitudes,and overestimated over the tropics,while IAP simulated more bare soil over almost all latitudes except southern middle latitudes.Over high latitudes,such as northeastern Canada,the underestimated bare soil of CTL mainly resulted from its overestimated grasses(Fig.2g).IAP significantly reduced the fractional coverage of grasses over northeastern Canada(Fig.2f);however,other vegetation,such as boreal shrubs,did not grow in this region(Fig.2d).Consequently,the bare soil of IAP was overestimated over northern high latitudes(Fig.2h).In the tropics,both CTL and IAP simulated more bare soil than observations(Fig.1d).These overestimations were mainly in northeastern Africa,southern Africa and Australia for CTL,and most regions of Africa for IAP(Fig.2h).The underestimated grasses were mostly responsible for these biases of tropical bare soil.

    To investigate the contribution of each PFT to the improvements,Fig.3 shows the global average fractional coverage for each PFT in the two simulations and the observation.The total coverage of IAP’s trees(26.51%)was more consistent with that of the observation(28.51%)than that of CTL(35.60%).The smaller fractional tree coverage in IAP was mainly contributed to by a reduction in“broadleaf evergreen tropical tree”(BET;3.17%),“broadleaf deciduous tropical tree”(BDT;3.27%),“broadleaf deciduous temperate tree”(BDM;2.10%)and“broadleaf deciduous boreal tree”(BDB;2.10%).However,IAP’s BDT was less than half that of the observation,which was the major contributor to the underestimation of IAP’s total tree coverage.Further investigation indicated that the new establishment parameterization of IAP-DGVM resulted in this underestimated BDT,mainly over tropical semi-arid regions(Fig.S1 in electronic supplementary material).IAP also simulated 1.32%more“needleleaf evergreen boreal tree”(NEB)than CTL,which contributed to the band of increased tree coverage over central Eurasia apparent in Fig.2b.In terms of shrubs,the increased“broadleaf deciduous boreal shrub”(BDBsh)coverage(5.44%)in IAP was the main contributor to its better agreement with the observation than CTL.For grasses,although the total fraction in CTL was closer to that observed than IAP’s grass fraction,the“C3 arctic grass”(C3Ar)coverage in CTL was 6.33%larger than the observation,which corresponds to the severely overestimated grasses over high latitudes shown in Fig.1c.The composition of IAP’s grasses was generally in good agreement with the observation.However,the“C4 grass”(C4)coverage in IAP was 4.64%less than the observation,resulting in the underestimation of total grasses in IAP.The bare soil of CTL was close to that observed,but IAP simulated 7.67%more bare soil than the observation.The overestimated bare soil in IAP mainly resulted from its underestimated shrubs in arctic regions and underestimated grasses in the tropics(Fig.2).

    Fig.3. Global weighted average fractional coverage(%)of each PFT for CTL(blue),IAP(red)and observation(black).The abbreviations correspond to the information in Table 1.

    Fig.4. Global distribution of the dominant vegetation type obtained from(a)CTL,(b)IAP and(c)observation.The abbreviations correspond to the information in Table 1.

    Figure 4 shows the global distribution of the dominant vegetation type,which is the PFT with the highest fractional coverage.Clearly,the dominant vegetation simulated by IAP was more consistent with that obtained from the observation than that of CTL.In CTL,C3Ar dominated over most northern high-latitude regions.However,IAP showed the northern high latitudes to be dominated by bare soil,BDBsh and NEB,which was a sim ilar situation to that observed.In the tropics,CTL overestimated regions dominated by BET,as compared to observation.For example,BET was the dom inant vegetation in CTL over southeastern South America and West Africa,which,according to observation,are actually dominated by C4.With respect to CTL,IAP simulated fewer regions dominated by BET,which agreed well with observation.Nonetheless,IAP’s C4 was also not the dom inant vegetation in southeastern South America and West A frica because of the underestimated C4(Fig.3).

    3.2. Carbon fluxes

    Compared to CTL,IAPsimulated an overallmore reasonable distribution of key carbon fluxes.Over most latitudes,the GPP in IAP was sim ilar to that of Beer et al.(2010),while CTL’s GPP was overestimated in middle and high latitudes(Fig.5a).The overestimated GPP in CTL stretched over the whole of central and eastern North America,Europe,South Asia and southeastern South America,while underestimated GPP dominated over the Amazon and Africa(Fig.6a).Relative to CTL,IAP’s GPP was closer to observation,especially over the Amazon(Fig.6b).Both CTL and IAP underestimated GPP over most regions of A frica.

    Fig.5. Zonal average(a)GPP,(b)NPP and(c) fire carbon emissions(FireC)in CTL(blue),IAP(red)and the benchmarks(OBS;black).A ll units are gC m?2 yr?1.

    In terms of NPP,IAP showed a better agreement with MODIS than CTL over northern high latitudes(Fig.5b).Relative to CTL,IAP simulated lower NPP over northeastern Canada and central and eastern Russia,where CTL simulated~200 gC m?2yr?1more NPP than the observation(Fig.6c).Both CTL and IAP underestimated NPP over middle latitudes,mainly over arid and semi-arid regions.In the tropics,the simulated NPP in CTL was higher than the observation,mainly because of the overestimated NPP over the Amazon,West A frica and the Maritime Continent.On the contrary,IAP’s NPP was consistent with observations over the Amazon,Africa and the Maritime Continent.Both CTL and IAP underestimated the NPP over most parts of Africa.

    Fire carbon emissions were overestimated by IAP over middle latitudes,especially in central and western North America,northeastern China,southern South America and Australia(Fig.6f).However,the fire carbon emissions simulated by IAP were much more consistent with observation in the tropics,where CTL showed a severe underestimation(Fig.5c).Fire carbon emissions in CTL were 10 gC m?2yr?1higher than the observation in eastern North America,Europe,the Amazon,Southeast Asia and the Maritime Continent,and 50 gC m?2yr?1lower in central South America and most regions of Africa(Fig.6e).IAP reduced these errors to different degrees by increasing or decreasing fire carbon emissions in these regions,respectively.Broadly,IAP captured a better spatial distribution of fire carbon emissions than CTL.

    Overall,within the framework of CAS-ESM,IAP simulated carbon fluxes closer to observations than CTL,as summarized in Fig.7.The GPP simulated by IAP was 150.5 PgC yr?1,which is closer to the 123±8 PgC yr?1reported by Beer et al.(2010)than CTL’s GPP,and comparable to the range from 101 to 150 PgC yr?1published elsewhere(Farquhar et al.,1993;Ciais et al.,1997).Meanwhile,IAP also overestimated autotrophic respiration,which was almost the same as its counterpart in CTL.Consequently,IAP’s NPP,59.1 PgC yr?1,compared better with the expected value of 60 PgC yr?1(Castillo et al.,2012)than the result of CTL(75.31 PgC yr?1).Moreover,heterotrophic respiration in IAP was 56.20 PgC yr?1,which is closer to the 55.4 PgC yr?1from IPCC(2013)than that of CTL(74.08 PgC yr?1).Thus,the net ecosystem production in IAP was more reasonable than that in CTL,in comparison to the baseline from IPCC(2013).The fire carbon emissions in IAP were slightly higher than those of GFEDv4,because of the overestimated fire carbon emissions in the midlatitudes.Consequently,the net biome production(NBP)of IAP-DGVM was?0.2 PgC yr?1,which is outside the range of 2.63±1.22 PgC yr?1reported by other process-based terrestrial ecosystem models driven by rising CO2and by changes in climate(IPCC,2013).This nearzero value of NBP could be acceptable,however,because the results were based on the equilibrium state,which was cyclically forced by atmospheric datasets and a constant CO2value(Castillo et al.,2012).

    3.3. LAI

    Generally,both CTL and IAP overestimated LAI,although the bias in IAP was smaller than that in CTL(Fig.8).The simulated annual mean LAI in CTL and IAP was 1.0 m2m?2more than the observation over most of the northern middle and high latitudes,such as central and eastern North America,Europe,central Eurasia and southeastern China.In the tropics,CTL’s bias in LAI exceeded 5.0 m2m?2,while IAP’s was~3.0 m2m?2.In terms of seasonal variability,both CTL and IAP were consistent with observations,the largest being during June–August and the smallest during December–February.However,the simulated LAI ampli-tudes in CTL and IAP were around tw ice as large as those observed for each month.Although the LAI values in IAP were closer to observation compared to those of CTL,the improvements were not remarkable.Therefore,it is necessary to further investigate the causes of these systematically overestimated LAI values.

    Fig.6. Differences between CTL and the benchmarks(CTL minus benchmarks)in(a)GPP,(c)NPP and(e) fire carbon emissions(FireC).(b,d,f)As in(a,c,e),respectively,but between IAP and the benchmarks.All units are gC m?2 yr?1.

    4. Conclusions

    This work evaluated the performance of IAP-DGVM within the framework of CAS-ESM through comparison with that of CoLM-DGVM,as well as observations and benchmarks.The results sufficiently demonstrated that IAPDGVM can simulate a realistic vegetation distribution,including spatial patterns,total areas and compositions,as well as reasonable carbon fluxes,such as GPP,NPP and fire carbon emissions.

    The total tree coverage of IAP-DGVM was found to be in good agreement with observations,because of the reduced tropical trees and increased NEB relative to CoLM-DGVM.The shrub coverage in IAP-DGVM also showed a sim ilar distribution to that observed,which resulted from the significantly increased fractional coverage of BDBsh,with replacement of C3Ar.Meanwhile,the reduced C3Ar was the major contributor to the better consistency between the grass coverage in IAP-DGVM and that observed than between CoLMDGVM and that observed.Consequently,the global distribution of the dom inant vegetation type simulated by IAPDGVM was similar to that observed,especially over northern high latitudes.Moreover,the biases of IAP-DGVM in terms of GPP and NPP were smaller because of improvements in the GPP over middle and high latitudes,as well as in the tropics.The tropical patterns of fire carbon emissions in IAP-DGVM were much more consistent than CTL with observations.These better performances of IAP-DGVM in sim-ulating the global vegetation distribution and carbon fluxes provide a foundation to use CAS-ESM to study vegetation–climate interactions.

    Fig.7. Global means of carbon fluxes in CTL and IAP,as well as that of the benchmarks.Units:PgC yr?1.

    Fig.8. Differences in annual mean LAI between(a)CTL and observations,and(b)IAP and observations.(c)Globally averaged LAI in CTL(blue),IAP(red)and observations(black)for each month.A ll units are m2 m?2.

    Despite the above positive results,deficiencies in IAPDGVM were also found.Speci fi cally,BDT was severely underestimated because of the new establishment parameterization;BDBsh could not grow in northeastern Canada because the simulated GDD5 did not exceed the threshold;and the C4 coverage was much smaller than the observation.Therefore,improved parameterization is necessary to simulate a more realistic distribution of vegetation types.Furthermore,both simulations overestimated GPP and autotrophic respiration,which is likely associated to the parameterization in CoLM being insufficiently sensitive to the DGVM.Additionally,the fire emissions simulated by the model in the high and middle latitudes were high.An understanding of the underlying mechanisms of these biases is needed to further improve the model,which will be the subject of future research.

    The primary purpose of this study was to evaluate the performance of IAP-DGVM within the framework of CASESM,which makes the evaluations more model specific.However,this work also exerts the following influences on other studies.This study reported the improvements and deficiencies of IAP-DGVM in CAS-ESM,which is valuable information for the application of CAS-ESM,as well as a sample for other model communities in terms of a comparative benchmark.The selection of global observations of carbon fluxes in this study was limited by the spatial and temporal scale of existing datasets,which is a pivotal message for observational scientists to observe carbon fluxes with large spatial scale and continuous time.Overall,we hope that the advantages and disadvantages of the simulations reported in this paper will prove valuable for scientists seeking to investigate climate change.

    Acknowledgements.This work was supported by the National Major Research High Performance Computing Program of China(Grant No.2016YFB02008)and the National Natural Science Foundation of China(Grant Number 41705070).Fang LI and Xiang SONG are supported by the National Natural Science Foundation of China(Grant Numbers 41475099 and 41305096).

    Electronic supplementary material:Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-017-7154-7.

    Baker,I.T.,L.Prihodko,A.S.Denning,M.Goulden,S.Miller,and H.R.da Rocha,2008:Seasonal drought stress in the Amazon:Reconciling models and observations.J.Geophys.Res.,113,G00B01,https://doi.org/10.1029/2007JG000644.

    Beer,C.,and Coauthors,2010:Terrestrial gross carbon dioxide uptake:Global distribution and covariation with climate.Science,329,834–838,https://doi.org/10.1126/science.1184984.

    Betts,R.A.,P.M.Cox,M.Collins,P.P.Harris,C.Hunting ford,and C.D.Jones,2004:The role of ecosystem-atmosphere interactions in simulated Amazonian precipitation decrease and forest dieback under global climate warming.Theor.Appl.Climatol.,78,157–175,https://doi.org/10.1007/s00704-004-0050-y.

    Blok,D.,M.M.P.D.Heijmans,G.Schaepman-Strub,A.V.Kononov,T.C.Maximov,and F.Berendse,2010:Shrub expansion may reduce summer permafrost thaw in Siberian tundra.Global Change Biology,16,1296–1305,https://doi.org/10.1111/j.1365-2486.2009.02110.x.

    Bonan,G.B.,1996:A land surface model(LSM Version 1.0)for ecological,hydrological,and atmospheric studies:Technical description and user’s guide.NCAR Tech.Note NCAR/TN-417+STR,https://doi.org/10.5065/D6DF6P5X.

    Bonfils,C.J.W.,T.J.Phillips,D.M.Lawrence,P.Cameron-Sm ith,W.J.Riley,and Z.M.Subin,2012:On the influenc of shrub height and expansion on northern high latitude climate.Environmental Research Letters,7,015503,https://doi.org/10.1088/1748-9326/7/1/015503.

    Castillo,C.K.G.,S.Levis,and P.Thornton,2012:Evaluation of the new CNDV option of the community land model:Effects of dynamic vegetation and interactive nitrogen on CLM 4 means and variability.J.Climate,25,3702–3714,https://doi.org/10.1175/JCLI-D-11-00372.1.

    Charney,J.,P.H.Stone,and W.J.Quirk,1975:Drought in the Sahara:A biogeophysical feedback mechanism.Science,187,434–435,https://doi.org/10.1126/science.187.4175.434.

    Ciais,P.,and Coauthors,1997:A three-dimensional synthesis study of δ18O in atmospheric CO2:1.Surface fluxes.J.Geophys.Res.,102,5857–5872,https://doi.org/10.1029/96JD 02360.

    Cramer,W.,and Coauthors,2001:Global response of terrestrial ecosystem structure and function to CO2and climate change:Results from six dynamic global vegetation models.Global Change Biology,7,357–373,https://doi.org/10.1046/j.1365-2486.2001.00383.x.

    Dai,Y.J.,and Q.C.Zeng,1997:A land surface model(IAP94)for climate studies Part I:Formulation and validation in off-line experiments.Adv.Atmos.Sci.,14,433–460,https://doi.org/10.1007/s00376-997-0063-4.

    Dai,Y.J.,and Coauthors,2003:The common land model.Bull.Amer.Meteorol.Soc.,84,1013–1023,https://doi.org/10.1175/BAMS-84-8-1013.

    Dickinson,R.E.,A.Henderson-Sellers,and P.J.Kennedy,1993:Biosphere-Atmosphere Transfer Scheme(BATS)Version 1e as Coupled to the NCAR Community Climate Model.NCAR Tech.Note NCAR/TN-387+STR,72 pp,https://doi.org/10.5065/D67W 6959.

    Erfanian,A.,G.L.Wang,M.Yu,and R.Anyah,2016:Multimodel ensemble simulations of present and future climates over West Africa:Impacts of vegetation dynamics.Journal of Advances in Modeling Earth Systems,8,1411–1431,https://doi.org/10.1002/2016MS000660.

    Falloon,P.D.,R.Dankers,R.A.Betts,C.D.Jones,B.B.B.Booth,and F.H.Lambert,2012:Role of vegetation change in future climate under the A1B scenario and a climate stabilisation scenario,using the HadCM 3C Earth system model.Biogeosciences,9,4739–4756,https://doi.org/10.5194/bg-9-4739-2012.

    Farquhar,G.D.,J.Lloyd,J.A.Taylor,L.B.Flanagan,J.P.Syvertsen,K.T.Hubick,S.C.Wong,and J.R.Ehleringer,1993:Vegetation effects on the isotope composition of oxygen in atmospheric CO2.Nature,363,439–443,https://doi.org/10.1038/363439a0.

    Fraser,R.H.,I.Olthof,M.Carri`ere,A.Deschamps,and D.Pouliot,2011:Detecting long-term changes to vegetation in northern Canada using the Landsat satellite image archive.Environmental Research Letters,6,045502.

    Friedlingstein,P.,and Coauthors,2006:Climate-carbon cycle feedback analysis:Results from the C4M IP model intercomparison.J.Climate,19,3337–3353,https://doi.org/10.1175/JCLI3800.1.

    Guillevic,P.,R.D.Koster,M.J.Suarez,L.Bounoua,G.J.Col-latz,S.O.Los,and S.P.P.Mahanama,2002:influenc of the interannual variability of vegetation on the surface energy balance—A global sensitivity study.Journal of Hydrometeorology,3,617–629,https://doi.org/10.1175/1525-7541(2002)003<0617:IOTIVO>2.0.CO;2.

    Holling,C.S.,1973:Resilience and stability of ecological systems.Annual Review of Ecology and Systematics,4,1–23,https://doi.org/10.1146/annurev.es.04.110173.000245.

    IPCC,2013:Climate Change 2013:The Physical Science Basis.Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change.T.F.Stockeretal.,Eds.,Cambridge University Press,1535 pp.

    Ji,D.,and Coauthors,2014:Description and basic evaluation of Beijing Normal University Earth System Model(BNU-ESM)version 1.Geoscienti fi c Model Development,7,2039–2064,https://doi.org/10.5194/gmd-7-2039-2014.

    Law rence,D.,and K.Vandecar,2015:Effects of tropical deforestation on climate and agriculture.Nat.Clim.Change,5,27–36,https://doi.org/10.1038/nclimate2430.

    Law rence,P.J.,and T.N.Chase,2007:Representing a new MODIS consistent land surface in the Community Land Model(CLM 3.0).J.Geophys.Res.,112,G01023,https://doi.org/10.1029/2006JG000168.

    Levis,S.,G.B.,Bonan,M.Vertenstein,and K.Oleson,2004:The Community Land Model’s dynam ic global vegetation model(CLM-DGVM):Technical description and user’s guide.NCAR Tech.Note TN-459+IA,50 pp,https://doi.org/10.5065/D6P26W 36.

    Li,F.,and D.M.Law rence,2017:Role of fire in the global land water budget during the twentieth century due to changing ecosystems.J.Climate,30,1893–1908,https://doi.org/10.1175/JCLI-D-16-0460.1.

    Li,F.,X.D.Zeng,and S.Levis,2012:A process-based fire parameterization of intermediate complexity in a Dynam ic Global Vegetation Model.Biogeosciences,9,2761–2780,https://doi.org/10.5194/bg-9-2761-2012.

    Li,F.,S.Levis,and D.S.Ward,2013:Quantifying the role of fire in the Earth system—Part1:Improved global fire modeling in the Community Earth System Model(CESM 1).Biogeosciences,10,2293–2314,https://doi.org/10.5194/bg-10-2293-2013.

    Li,F.,B.Bond-Lamberty,and S.Levis,2014:Quantifying the role of fire in the Earth system—Part 2:Impact on the net carbon balance of global terrestrial ecosystems for the 20th century.Biogeosciences,11,1345–1360,https://doi.org/10.5194/bg-11-1345-2014.

    Li,F.,D.M.Law rence,and B.Bond-Lamberty,2017:Impact of fire on global land surface air temperature and energy budget for the 20th century due to changes within ecosystems.Environmental Research Letters,12,044014,https://doi.org/10.1088/1748-9326/aa6685.

    Li,Q.,and Y.K.Xue,2010:Simulated impacts of land cover change on summer climate in the Tibetan Plateau.Environmental Research Letters,5,015102,https://doi.org/10.1088/1748-9326/5/1/015102.

    Myers-Smith,I.H.,and Coauthors,2011:Shrub expansion in tundra ecosystems:dynamics,impacts and research priorities.Environmental Research Letters,6,045509,https://doi.org/10.1088/1748-9326/6/4/045509.

    Peng,C.H.,2000:From static biogeographical model to dynam ic global vegetation model:A global perspective on modelling vegetation dynamics.Ecological Modelling,135(1),33–54,https://doi.org/10.1016/S0304-3800(00)00348-3.

    Pielke,R.A.,R.Avissar,M.Raupach,A.J.Dolman,X.B.Zeng,and A.S.Denning,1998:Interactions between the atmosphere and terrestrial ecosystems:influenc on weather and climate.Global Change Biology,4,461–475,https://doi.org/10.1046/j.1365-2486.1998.t01-1-00176.x.

    Qian,T.T.,A.G.Dai,K.E.Trenberth,and K.W.Oleson,2006:Simulation of global land surface conditions from 1948 to 2004.Part I:Forcing data and evaluations.Journal of Hydrometeorology,7,953–975,https://doi.org/10.1175/JHM 540.1.

    Quillet,A.,C.H.Peng,and M.Garneau,2010:Toward dynamic global vegetation models for simulating vegetation-climate interactions and feedbacks:Recent developments,limitations,and future challenges.Environmental Reviews,18,333–353,https://doi.org/10.1139/A10-016.

    Randerson,J.T.,G.R.Van Der Werf,L.Giglio,G.J.Collatz,and P.S.Kasibhatla.2015:Global Fire Emissions Database,Version 4,(GFEDv4).ORNL DAAC,Oak Ridge,Tennessee,USA,https://dx.doi.org/10.3334/ORNLDAAC/1293.

    Restrepo-Coupe,N.,and Coauthors,2017:Do dynamic global vegetation models capture the seasonality of carbon fluxes in the Amazon basin?A data-model inter comparison.Global Change Biology,23,191–208,https://doi.org/10.1111/gcb.13442.

    Scheffer,M.,and Coauthors,2009:Early-warning signals for critical transitions.Nature,461,53–59,https://doi.org/10.1038/nature08227.

    Seddon,A.W.R.,M.Macias-Fauria,P.R.Long,D.Benz,and K.J.Willis,2016:Sensitivity of global terrestrial ecosystems to climate variability.Nature,531,229–232,https://doi.org/10.1038/nature16986.

    Sitch,S.,and Coauthors,2003:Evaluation of ecosystem dynamics,plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model.Global Change Biology,9,161–185,https://doi.org/10.1046/j.1365-2486.2003.00569.x.Song,X.,X.D.Zeng,J.W.Zhu,and P.Shao,2016:Development of an establishment scheme for a DGVM.Adv.Atmos.Sci.,33,829–840,https://doi.org/10.1007/s00376-016-5284-y.

    Sturm,M.,J.Holmgren,J.P.McFadden,G.E.Liston,F.S.Chapin III,and C.H.Racine,2001:Snow-shrub interactions in Arctic tundra:A hypothesis with climatic implications.J.Climate,14,336–344,https://doi.org/10.1175/1520-0442(2001)014<0336:SSIIAT>2.0.CO;2.

    Xue,Y.K.,F.De Sales,R.Vasic,C.R.Mechoso,A.Arakawa,and S.Prince,2010:Global and seasonal assessment of interactions between climate and vegetation biophysical processes:AGCM study with differentl and-vegetation representations.J.Climate,23,1411–1433,https://doi.org/10.1175/2009JCLI3054.1.

    Yu,M.,G.L.Wang,and J.S.Pal,2016:Effects of vegetation feedback on future climate change over West A frica.Climate Dyn.,46,3669–3688,https://doi.org/10.1007/s00382-015-2795-7.

    Zeng,X.D.,2010:Evaluating the dependence of vegetation on climate in an improved dynamic global vegetation model.Adv.Atmos.Sci.,27,977–991,https://doi.org/10.1007/s00376-009-9186-0.

    Zeng,X.D.,F.Li,and X.Song,2014:Development of the IAP dynam ic global vegetation model.Adv.Atmos.Sci.,31,505–514,https://doi.org/10.1007/s00376-013-3155-3.

    Zeng,X.D.,X.B.Zeng,and M.Barlage,2008:Grow ing tem-perate shrubs over arid and sem iarid regions in the Community Land Model-Dynam ic Global Vegetation Model.Global Biogeochemical Cycles,22,GB3003,https://doi.org/10.1029/2007GB003014.

    Zhao,M.S.,and S.W.Running,2010:Drought-induced reductionin global terrestrial net primary production from 2000 through 2009.Science,329,940–943,https://doi.org/10.1126/science.1192666.

    Zhu,J.W.,and X.D.Zeng,2015:Comprehensive study on the influence of evapotranspiration and albedo on surface temperature related to changes in the leaf area index.Adv.Atmos.Sci.,32(7),935–942,https://doi.org/10.1007/s00376-014-4045-z.

    Zhu,J.W.,and X.D.Zeng,2016:influences of the interannual variability of vegetation LAI on surface temperature.Atmos.Oceanic Sci.Lett.,9(4),292–297,https://dx.doi.org/10.1080/16742834.2016.1189800.

    Zhu,J.W.,and X.D.Zeng,2017:influences of the seasonal grow th of vegetation on surface energy budgets over middle to high latitudes.International Journal of Climatology,37,4251–4260,https://doi.org/10.1002/joc.5068.

    成人18禁在线播放| 性色av乱码一区二区三区2| 久久欧美精品欧美久久欧美| 美女扒开内裤让男人捅视频| 欧洲精品卡2卡3卡4卡5卡区| 久久国产精品人妻蜜桃| 久久国产精品人妻蜜桃| 久久影院123| 久久伊人香网站| 中文字幕精品免费在线观看视频| 国内久久婷婷六月综合欲色啪| 久久精品aⅴ一区二区三区四区| 十八禁网站免费在线| 亚洲九九香蕉| 亚洲欧美一区二区三区黑人| 国产精华一区二区三区| 国产极品粉嫩免费观看在线| 操出白浆在线播放| 国产亚洲欧美精品永久| 精品久久蜜臀av无| 在线av久久热| 少妇的丰满在线观看| 久久青草综合色| 成人永久免费在线观看视频| 久久久久国产精品人妻aⅴ院| 777久久人妻少妇嫩草av网站| 亚洲 欧美 日韩 在线 免费| 亚洲黑人精品在线| 国产aⅴ精品一区二区三区波| 国产欧美日韩综合在线一区二区| 别揉我奶头~嗯~啊~动态视频| 国产精品秋霞免费鲁丝片| 成人三级做爰电影| 日韩视频一区二区在线观看| www.999成人在线观看| 亚洲av第一区精品v没综合| 国产精华一区二区三区| 女警被强在线播放| 国产成人av教育| 多毛熟女@视频| 久久热在线av| 99久久99久久久精品蜜桃| av免费在线观看网站| 午夜影院日韩av| 男女床上黄色一级片免费看| 免费观看人在逋| 午夜老司机福利片| 在线国产一区二区在线| 人妻丰满熟妇av一区二区三区| 日本撒尿小便嘘嘘汇集6| 精品午夜福利视频在线观看一区| 女人精品久久久久毛片| 国产激情久久老熟女| 亚洲中文日韩欧美视频| 99在线人妻在线中文字幕| 中亚洲国语对白在线视频| 欧美一级毛片孕妇| 真人做人爱边吃奶动态| 黄色 视频免费看| 精品人妻在线不人妻| 亚洲av日韩精品久久久久久密| 1024视频免费在线观看| 成人亚洲精品av一区二区 | 成人18禁高潮啪啪吃奶动态图| 啦啦啦在线免费观看视频4| 亚洲精品av麻豆狂野| 国产亚洲精品第一综合不卡| 精品人妻在线不人妻| 女性生殖器流出的白浆| 亚洲一区二区三区欧美精品| 国产成人影院久久av| 成人国语在线视频| 国产av一区在线观看免费| 91精品国产国语对白视频| 黄色a级毛片大全视频| 美女国产高潮福利片在线看| 麻豆久久精品国产亚洲av | 亚洲人成电影观看| 女人精品久久久久毛片| 男人的好看免费观看在线视频 | 久久国产精品人妻蜜桃| 久久精品亚洲精品国产色婷小说| 淫妇啪啪啪对白视频| 女同久久另类99精品国产91| 美女大奶头视频| 亚洲精品国产精品久久久不卡| 国产精品久久久久久人妻精品电影| 这个男人来自地球电影免费观看| 曰老女人黄片| 女生性感内裤真人,穿戴方法视频| 一边摸一边抽搐一进一小说| 香蕉国产在线看| 精品久久久久久电影网| www.精华液| 亚洲 欧美一区二区三区| 亚洲精品av麻豆狂野| 日日摸夜夜添夜夜添小说| 午夜福利影视在线免费观看| 一区二区三区激情视频| 亚洲欧美日韩无卡精品| 午夜免费观看网址| 51午夜福利影视在线观看| 国产高清激情床上av| 长腿黑丝高跟| 日本a在线网址| 后天国语完整版免费观看| 在线av久久热| 久久国产亚洲av麻豆专区| 亚洲av片天天在线观看| 国产深夜福利视频在线观看| 色综合欧美亚洲国产小说| netflix在线观看网站| 久久久国产欧美日韩av| 交换朋友夫妻互换小说| 国产人伦9x9x在线观看| 成人国产一区最新在线观看| 69av精品久久久久久| 天堂影院成人在线观看| 亚洲全国av大片| 成人亚洲精品一区在线观看| 大型av网站在线播放| 最近最新中文字幕大全免费视频| 久久99一区二区三区| 超色免费av| 欧美日韩一级在线毛片| 国产激情久久老熟女| 午夜91福利影院| 一级黄色大片毛片| 久久久国产成人精品二区 | 国产欧美日韩综合在线一区二区| 精品午夜福利视频在线观看一区| 9色porny在线观看| 激情在线观看视频在线高清| 亚洲av成人不卡在线观看播放网| 一区福利在线观看| 黄色视频,在线免费观看| 亚洲av第一区精品v没综合| 亚洲va日本ⅴa欧美va伊人久久| 精品久久久久久久久久免费视频 | 久久青草综合色| 色老头精品视频在线观看| 精品电影一区二区在线| 亚洲男人天堂网一区| 久久久久久久精品吃奶| 十分钟在线观看高清视频www| 人人澡人人妻人| 曰老女人黄片| 乱人伦中国视频| 亚洲成人国产一区在线观看| 最近最新免费中文字幕在线| 女警被强在线播放| 欧美国产精品va在线观看不卡| 国产野战对白在线观看| 18美女黄网站色大片免费观看| 精品一区二区三区av网在线观看| 国产三级在线视频| 18禁裸乳无遮挡免费网站照片 | a级毛片在线看网站| 久久影院123| 每晚都被弄得嗷嗷叫到高潮| 精品福利永久在线观看| 麻豆一二三区av精品| 亚洲自拍偷在线| a级毛片在线看网站| 国产在线精品亚洲第一网站| 日韩成人在线观看一区二区三区| www.www免费av| 国产99久久九九免费精品| 最好的美女福利视频网| 嫁个100分男人电影在线观看| 国产激情欧美一区二区| 999久久久精品免费观看国产| 免费高清视频大片| 久久久国产一区二区| 亚洲 欧美 日韩 在线 免费| 狠狠狠狠99中文字幕| 国产aⅴ精品一区二区三区波| 每晚都被弄得嗷嗷叫到高潮| 日韩人妻精品一区2区三区| 久久青草综合色| 欧美日韩乱码在线| 国产激情久久老熟女| 亚洲国产欧美日韩在线播放| 少妇裸体淫交视频免费看高清 | 国产欧美日韩一区二区三区在线| 国产成人精品无人区| 脱女人内裤的视频| 欧美av亚洲av综合av国产av| 91成年电影在线观看| 免费女性裸体啪啪无遮挡网站| 1024视频免费在线观看| 精品欧美一区二区三区在线| 中文字幕av电影在线播放| 中文亚洲av片在线观看爽| 午夜日韩欧美国产| 韩国精品一区二区三区| 亚洲av美国av| 成人手机av| 国产精品亚洲av一区麻豆| 一级毛片高清免费大全| 欧美亚洲日本最大视频资源| 日韩欧美一区二区三区在线观看| 国产免费现黄频在线看| 精品久久蜜臀av无| av电影中文网址| 天天影视国产精品| 国产av又大| 又黄又粗又硬又大视频| 黄片小视频在线播放| 日本撒尿小便嘘嘘汇集6| 在线观看免费视频日本深夜| 免费女性裸体啪啪无遮挡网站| 久久久久精品国产欧美久久久| 成人av一区二区三区在线看| 亚洲一区二区三区欧美精品| 久久 成人 亚洲| 亚洲第一av免费看| 中文字幕人妻熟女乱码| 97人妻天天添夜夜摸| 亚洲精品国产精品久久久不卡| 国产伦人伦偷精品视频| 两性夫妻黄色片| 欧美性长视频在线观看| av有码第一页| 欧洲精品卡2卡3卡4卡5卡区| 精品电影一区二区在线| 色精品久久人妻99蜜桃| 国产免费男女视频| 琪琪午夜伦伦电影理论片6080| 精品欧美一区二区三区在线| 亚洲专区中文字幕在线| www.www免费av| 国产伦一二天堂av在线观看| 在线视频色国产色| 久久久久九九精品影院| 精品一区二区三区四区五区乱码| xxx96com| 亚洲熟女毛片儿| 免费在线观看完整版高清| 在线观看午夜福利视频| 丰满迷人的少妇在线观看| 变态另类成人亚洲欧美熟女 | 国产又爽黄色视频| 日韩av在线大香蕉| 久久久久精品国产欧美久久久| 麻豆成人av在线观看| 欧美黄色片欧美黄色片| 国产欧美日韩一区二区精品| av国产精品久久久久影院| 老司机靠b影院| 欧美黄色片欧美黄色片| www.www免费av| xxxhd国产人妻xxx| 一级a爱片免费观看的视频| 精品熟女少妇八av免费久了| 成人特级黄色片久久久久久久| 99国产综合亚洲精品| 在线播放国产精品三级| 国产不卡一卡二| 女警被强在线播放| 亚洲国产欧美日韩在线播放| 18美女黄网站色大片免费观看| 色播在线永久视频| 久久香蕉国产精品| 日韩大尺度精品在线看网址 | 国内毛片毛片毛片毛片毛片| 亚洲色图综合在线观看| 国产精品国产av在线观看| 一进一出抽搐动态| 麻豆久久精品国产亚洲av | 久久香蕉精品热| 看免费av毛片| 国产精品久久久人人做人人爽| 99热国产这里只有精品6| 日韩精品中文字幕看吧| 精品免费久久久久久久清纯| 91在线观看av| 国产精品爽爽va在线观看网站 | 亚洲五月天丁香| 美国免费a级毛片| 成年人免费黄色播放视频| 国产欧美日韩综合在线一区二区| 性少妇av在线| 激情视频va一区二区三区| 成人免费观看视频高清| 操出白浆在线播放| 黄色a级毛片大全视频| 岛国在线观看网站| 88av欧美| 老司机福利观看| 国产av在哪里看| 久久精品aⅴ一区二区三区四区| 亚洲狠狠婷婷综合久久图片| 久久久国产成人精品二区 | 亚洲情色 制服丝袜| 亚洲人成网站在线播放欧美日韩| www.999成人在线观看| 日本wwww免费看| 欧美不卡视频在线免费观看 | 男女之事视频高清在线观看| 91九色精品人成在线观看| 9热在线视频观看99| 日韩免费高清中文字幕av| 天天躁夜夜躁狠狠躁躁| 亚洲七黄色美女视频| 久久久国产成人精品二区 | 男女下面进入的视频免费午夜 | 一二三四在线观看免费中文在| 午夜福利在线免费观看网站| 成人精品一区二区免费| 欧美日本中文国产一区发布| 级片在线观看| 久久天堂一区二区三区四区| 中文字幕最新亚洲高清| 在线观看免费视频网站a站| 老司机午夜十八禁免费视频| 免费高清在线观看日韩| av网站免费在线观看视频| 婷婷六月久久综合丁香| 免费搜索国产男女视频| 脱女人内裤的视频| 在线观看免费视频网站a站| 欧美精品一区二区免费开放| 夜夜躁狠狠躁天天躁| 亚洲全国av大片| 国产伦人伦偷精品视频| 亚洲精品国产一区二区精华液| 欧美不卡视频在线免费观看 | 国产精品影院久久| 成人亚洲精品一区在线观看| 日本vs欧美在线观看视频| 久久久久国产精品人妻aⅴ院| 国产精品国产av在线观看| 麻豆国产av国片精品| 美国免费a级毛片| 18禁国产床啪视频网站| e午夜精品久久久久久久| 自线自在国产av| 香蕉国产在线看| 在线观看免费视频日本深夜| 男人的好看免费观看在线视频 | www.www免费av| 亚洲精品一卡2卡三卡4卡5卡| 日本免费a在线| 咕卡用的链子| 国内毛片毛片毛片毛片毛片| 久久精品亚洲熟妇少妇任你| 大型av网站在线播放| 欧美国产精品va在线观看不卡| 天天添夜夜摸| 亚洲专区中文字幕在线| 久久亚洲真实| 水蜜桃什么品种好| 久久精品91无色码中文字幕| 亚洲人成网站在线播放欧美日韩| 国产成人精品无人区| 一本综合久久免费| 天堂中文最新版在线下载| 亚洲av片天天在线观看| 国产有黄有色有爽视频| 欧美国产精品va在线观看不卡| av在线天堂中文字幕 | 日韩欧美一区二区三区在线观看| 日韩欧美在线二视频| 国产不卡一卡二| 人妻久久中文字幕网| 日韩 欧美 亚洲 中文字幕| 色综合欧美亚洲国产小说| 国产真人三级小视频在线观看| 夜夜夜夜夜久久久久| 免费看a级黄色片| 国产区一区二久久| 美女大奶头视频| 精品欧美一区二区三区在线| 日本 av在线| 亚洲国产中文字幕在线视频| 久久天堂一区二区三区四区| 最新美女视频免费是黄的| 黄片小视频在线播放| 国产精品一区二区免费欧美| 国产一区二区在线av高清观看| 人人妻人人添人人爽欧美一区卜| 99热只有精品国产| 黄色怎么调成土黄色| 亚洲av美国av| 中文字幕av电影在线播放| 免费高清视频大片| 亚洲中文日韩欧美视频| 99国产综合亚洲精品| 亚洲国产欧美网| 免费高清在线观看日韩| 搡老熟女国产l中国老女人| 国产精品自产拍在线观看55亚洲| 女人爽到高潮嗷嗷叫在线视频| 久久久久久久午夜电影 | 叶爱在线成人免费视频播放| 日韩国内少妇激情av| 亚洲aⅴ乱码一区二区在线播放 | 久久久国产成人免费| 亚洲情色 制服丝袜| 91字幕亚洲| 精品福利永久在线观看| 精品国内亚洲2022精品成人| 一区福利在线观看| bbb黄色大片| 日韩中文字幕欧美一区二区| 俄罗斯特黄特色一大片| 久久人人爽av亚洲精品天堂| 人人妻,人人澡人人爽秒播| 一区二区三区精品91| 国产亚洲精品一区二区www| 亚洲男人天堂网一区| 老司机福利观看| 99国产极品粉嫩在线观看| 国产精华一区二区三区| 亚洲一区中文字幕在线| 久久久久久大精品| 国产1区2区3区精品| 亚洲 国产 在线| 大香蕉久久成人网| 国产乱人伦免费视频| 麻豆av在线久日| 精品国产亚洲在线| 国产成人av教育| 国产亚洲精品第一综合不卡| 亚洲男人的天堂狠狠| 日韩欧美一区二区三区在线观看| 久久人妻av系列| 成年人免费黄色播放视频| 日日爽夜夜爽网站| 淫秽高清视频在线观看| 757午夜福利合集在线观看| 夜夜爽天天搞| 丰满饥渴人妻一区二区三| 欧美乱妇无乱码| 亚洲全国av大片| 最新在线观看一区二区三区| 亚洲视频免费观看视频| 人人妻人人添人人爽欧美一区卜| 国产一区二区在线av高清观看| 久久 成人 亚洲| 视频区欧美日本亚洲| 欧美乱色亚洲激情| 中亚洲国语对白在线视频| 大陆偷拍与自拍| 麻豆一二三区av精品| 精品国产乱码久久久久久男人| 欧美 亚洲 国产 日韩一| 满18在线观看网站| 黄色成人免费大全| 亚洲成a人片在线一区二区| 精品久久久久久电影网| 精品国产一区二区三区四区第35| 88av欧美| www.精华液| 成人三级黄色视频| 女警被强在线播放| 国产一区二区三区视频了| 久久亚洲真实| 国产精品乱码一区二三区的特点 | 亚洲第一青青草原| 男女高潮啪啪啪动态图| av中文乱码字幕在线| 久99久视频精品免费| 男女下面进入的视频免费午夜 | 亚洲 国产 在线| 少妇的丰满在线观看| 另类亚洲欧美激情| 欧美久久黑人一区二区| 久久久水蜜桃国产精品网| 精品国产国语对白av| 少妇裸体淫交视频免费看高清 | 夜夜爽天天搞| 免费在线观看影片大全网站| 黄色片一级片一级黄色片| 国产成人欧美| 高清毛片免费观看视频网站 | 精品人妻1区二区| 成年人黄色毛片网站| 亚洲精品在线美女| 亚洲精品国产一区二区精华液| 国产精品野战在线观看 | 欧美一级毛片孕妇| 女人被躁到高潮嗷嗷叫费观| 电影成人av| 手机成人av网站| 久久精品国产亚洲av高清一级| 在线观看免费高清a一片| 亚洲av美国av| 久久人妻av系列| 国产91精品成人一区二区三区| 狂野欧美激情性xxxx| 欧美精品一区二区免费开放| 成人三级做爰电影| 成人三级黄色视频| 性欧美人与动物交配| 岛国视频午夜一区免费看| 国产极品粉嫩免费观看在线| 亚洲 欧美 日韩 在线 免费| 亚洲av熟女| 欧美 亚洲 国产 日韩一| 日韩三级视频一区二区三区| xxx96com| 欧美激情 高清一区二区三区| 成熟少妇高潮喷水视频| 亚洲自拍偷在线| 最新在线观看一区二区三区| 97碰自拍视频| 日本wwww免费看| 久久欧美精品欧美久久欧美| 三级毛片av免费| 久久久久亚洲av毛片大全| x7x7x7水蜜桃| 男女午夜视频在线观看| 亚洲免费av在线视频| 国产亚洲精品第一综合不卡| 两人在一起打扑克的视频| av中文乱码字幕在线| 免费不卡黄色视频| 12—13女人毛片做爰片一| 日韩大码丰满熟妇| 精品少妇一区二区三区视频日本电影| 免费看十八禁软件| 婷婷精品国产亚洲av在线| 在线播放国产精品三级| 免费高清在线观看日韩| 水蜜桃什么品种好| 午夜老司机福利片| 黄色女人牲交| 亚洲午夜理论影院| 日韩大码丰满熟妇| 国产野战对白在线观看| 国产精品一区二区免费欧美| 欧美激情高清一区二区三区| 久久精品91蜜桃| 久久久国产成人精品二区 | 人人澡人人妻人| 国产精品自产拍在线观看55亚洲| 波多野结衣av一区二区av| 最近最新中文字幕大全电影3 | 国产欧美日韩一区二区三区在线| 成年女人毛片免费观看观看9| 久久精品国产清高在天天线| 在线观看舔阴道视频| 超色免费av| av天堂在线播放| 国产精品一区二区在线不卡| 精品久久久久久久毛片微露脸| 成年版毛片免费区| 国产成人av教育| 丰满的人妻完整版| 美国免费a级毛片| 成年人黄色毛片网站| 搡老岳熟女国产| 亚洲 欧美一区二区三区| 日韩av在线大香蕉| 黑人巨大精品欧美一区二区蜜桃| 日本免费a在线| 成人影院久久| 一区二区日韩欧美中文字幕| 搡老乐熟女国产| 男女做爰动态图高潮gif福利片 | 欧美人与性动交α欧美精品济南到| 在线播放国产精品三级| 国产免费av片在线观看野外av| 少妇的丰满在线观看| 村上凉子中文字幕在线| 欧美成人免费av一区二区三区| 一a级毛片在线观看| 丰满人妻熟妇乱又伦精品不卡| 在线永久观看黄色视频| 91老司机精品| 国产精品日韩av在线免费观看 | 咕卡用的链子| 国产成人精品无人区| 午夜激情av网站| 这个男人来自地球电影免费观看| 亚洲中文字幕日韩| 最近最新中文字幕大全电影3 | 99国产极品粉嫩在线观看| 99久久综合精品五月天人人| 国产一区二区在线av高清观看| 亚洲国产精品合色在线| 女警被强在线播放| 亚洲精品一二三| 不卡av一区二区三区| 99热国产这里只有精品6| 国产主播在线观看一区二区| 精品卡一卡二卡四卡免费| 一边摸一边抽搐一进一小说| 波多野结衣av一区二区av| av国产精品久久久久影院| 欧美在线一区亚洲| 国产精品影院久久| 黄色a级毛片大全视频| 久久亚洲真实| 国产野战对白在线观看| 国产99久久九九免费精品| 1024视频免费在线观看| 欧美av亚洲av综合av国产av| 啦啦啦在线免费观看视频4| 1024视频免费在线观看| 亚洲国产精品999在线| 国产1区2区3区精品| 高清毛片免费观看视频网站 | 国产极品粉嫩免费观看在线| 久久精品国产清高在天天线| 在线十欧美十亚洲十日本专区| 一级作爱视频免费观看| 久久天堂一区二区三区四区| 久久精品亚洲熟妇少妇任你| 三级毛片av免费| 国产精品久久电影中文字幕| 亚洲av电影在线进入| 妹子高潮喷水视频|