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

    Estimating the carbon fluxes of forests with an individual-based forest model

    2017-07-08 02:27:31EdnadigAndreasHuthFriedrichBohnCorinnaRebmannandMatthiasCuntz
    Forest Ecosystems 2017年2期

    Edna R?dig,Andreas Huth,Friedrich Bohn,Corinna Rebmannand Matthias Cuntz,4

    Estimating the carbon fluxes of forests with an individual-based forest model

    Edna R?dig1*,Andreas Huth1,2,3,Friedrich Bohn1,Corinna Rebmann1and Matthias Cuntz1,4

    Background:Capturing the response of forest ecosystems to inter-annual climate variability is a great challenge. In this study,we tested the capability of an individual-based forest gap model to display carbon fluxes at yearly and daily time scales.The forest model was applied to a spruce forest to simulate the gross primary production (GPP),respiration and net ecosystem exchange(NEE).We analyzed how the variability in climate affected simulated carbon fluxes at the scale of the forest model.

    Results:Six years were simulated at a daily time scale and compared to the observed eddy covariance(EC)data.In general,the seasonal cycle of the individual carbon fluxes was correctly described by the forest model.However,the estimated GPP differed from the observed data on the days of extreme climatic conditions.Two new parameterizations were developed:one resulting from a numerical calibration,and the other resulting from a filtering method.We suggest new parameter values and even a new function for the temperature limitation of photosynthesis.

    Conclusions:The forest model reproduced the observed carbon fluxes of a forest ecosystem quite well.Of the three parameterizations,the calibrated model version performed best.However,the filtering approach showed that calibrated parameter values do not necessarily correctly display the individual functional relations.The concept of simulating forest dynamics at the individual base is a valuable tool for simulating the NEE,GPP and respiration of forest ecosystems.

    Forest model,Temperate forest,Carbon fluxes,Eddy covariance,FORMIND

    Background

    Inter-annual climate variations can strongly influence the productivity of forest ecosystems.The heat wave of 2003,for example,caused a reduction of approximately 30%to the gross primary production(GPP)over Europe (Ciais et al.2005).This extreme event was followed by several studies to understand ecosystem responses and their underlying mechanisms(e.g.,Zaitchik et al.2006; Granier et al.2007).Models predict that such extreme events will become more frequent and intense in the future(Meehl and Tebaldi 2004).

    However,correctly capturing the responses of an ecosystem to varying climatic conditions with vegetation models is a major challenge(Keenan et al.2012).In this study,we test the potential of a forest gap model that considers forest structure at the local scale to estimate daily carbon fluxes and their response to climate variations in a spruce forest in Germany.

    Forest models have a long tradition in ecological research(Botkin et al.1972;Shugart 1984;Pacala et al. 1993;K?hler and Huth 2004).They have been successfully applied to analyze forest succession,tree species composition and biomass(e.g.,Fischer et al.2016). Capturing the competition between individuals enables these types of models to reproduce reality better than vegetation models which operate at a larger scale (Smith et al.2001).Traditionally,studies using forest models have focused on forest structure and dynamics, but they often neglected carbon exchange with the atmosphere(Bugmann 2001).

    In this study,we used an individual-based forest gap model(FORMIND)that simulates the growth of individual trees by calculating its photosynthesis and respiration (K?hler and Huth 1998).In addition,a soil carbon module is included(derived from Sato et al.2007).These model characteristics established a base to capture carbon fluxesat the ecosystem level,in addition to biomass and treesize distributions.

    Up to now,the forest model FORMIND has mainly been used to reproduce the average forest conditions in long-term studies at yearly time scales(Gutiérrez et al. 2009;Fischer et al.2014;Bohn et al.2014).The influence of short-term climate variability on individual model processes (GPP and respiration)has not yet been tested.Therefore, we here compared daily simulation output of the forest model with observed eddy covariance data of two spruce forests and analyzed the following two aspects.First,we wanted to test whether the local forest model is generally capable of displaying daily carbon fluxes.Second,we analyzed whether the model processes correctly respond to variable climate inputs.In this context,we evaluated the response of three different model parameterizations. We explored how we can use eddy covariance data to improve the simulation of carbon fluxes with an individualbased forest model.The simulation time also covered the heat wave of 2003 which enabled to include an extreme event in the analyses.

    The following questions will guide us through this study:(1)How well does an individual-based forest gap model simulate the daily and yearly carbon fluxes of a temperate forest ecosystem?(2)How can EC data be used to improve the concept of limiting factors in forest gap models?

    Methods

    The study sites

    This study focused on a forest site located at Wetzstein Mountain,part of the Thuringian Forest in central-east Germany where measured carbon fluxes and inventory data are available(Wetzstein flux tower,Rebmann et al. 2010).Observed carbon fluxes were derived with the eddy covariance(EC)method,a technique that observes the local carbon flux dynamics of the vegetation and monitors inter-annual changes(Aubinet et al.1999). The Wetzstein forest is dominated by even-aged Norway spruce(Picea abies)stands on clay loam.In addition,we analyzed another Norway spruce stand at Tharandt,a study site in the Ore Mountains in Germany where EC-data were available(Tharandt Anchor station,Grünwald and Bernhofer 2007).The stand characteristics of both sites are summarized in Table 1.

    The forest model FORMIND

    FORMIND(K?hler and Huth 2004;Fischer et al.2016) is an individual-based forest gap model in which growth is calculated for each tree individually.The approach uses patches to describe the vertical and horizontal forest structures.The main processes of the model include establishment,growth,mortality and competition.Important driving factors are daily means of incoming light(photosynthetic photon flux density,PPFD),temperature and precipitation(based on a model verison for temperate forests as in Bohn et al.2014).In this study,the model was applied to an even-aged spruce forest(1 ha).Establishment and mortality were deactivated for the short simulation time of 6 years.A full model description can be found in Fischer et al.(2016)and at www.formind.org.

    Table 1 Site characteristics for Wetzstein(Rebmann et al. 2010;Martina Mundt,personal communication)and Tharandt (Grünwald and Bernhofer 2007)

    The model runs with daily variable observed climate inputs of PPFD,day length,temperature and precipitation measured onsite.PPFD and day length serve as the driving forces for forest productivity.The sum of the GPP over all trees thus equals the GPP of the ecosystem. The ecosystem respiration is the sum of the respiration of all trees plus that of the soil and deadwood.The NEE is calculated as the difference between the ecosystem GPP and the ecosystem respiration(Fischer et al.2014). A positive NEE corresponds to increasing carbon stocks.

    Gross primary production

    Photosynthesis is calculated at the leaf level using a light-response function and is then integrated over the entire canopy(Thornley and Johnson 1990).The GPP of an individual tree under optimal climatic conditions (Huth and Ditzer 2000)equals

    in μmol(CO2)?m?2?s?1,where pmax(μmol(CO2)?m?2?s?1) is the maximum photosynthetic rate of the tree species (here,spruce),α is the initial slope of the light-responsecurve(μmol(CO2)?μmol(photons)–1),k is the light extinction factor,and m is the transmission coefficient of the leaves.Iindis the fraction of the PPFD at daily time step t that reaches the top of the individual tree.Ac(m2) is the crown area,and ψ(s)the photosynthetically active period of the time scale.Under non-optimal climatic conditions,GPPptreeis limited by the available soil water (SW)and temperature(T)(Bohn et al.2014):

    where φSWis the water reduction factor([0,1]),and φTthe temperature reduction factor([0,1]).

    The temperature reduction factor φTis derived from the LPJ-model(Sitch et al.2003)which includes two ramp functions(Gutiérrez et al.2012):

    where T(°C)is the daily mean air temperature at time step t.TCO2;h,TCO2;l,Tcoldand Thot(°C)are speciesspecific parameters representing the higher and lower temperature limits for CO2assimilation and the monthly mean air temperatures of the warmest and coldest months when production can still occur.

    In this study,we also tested a new temperature reduction curve φT*.It is distributed around the optimal temperature for photosynthesis Topt(°C)and the width Tsig(°C)(June et al.2004,reduction of the electron transport rate with n=2)since Eq.3 could not be properly fitted to the observed data.We suggest fitting this bell-shaped curve because it only relies on two parameters instead of four parameters(Eq.3):

    We use a water reduction factor,φSW,as proposed by Granier et al.(1999):

    where SWpwpis the permanent wilting point,SWmsw=SWpwp+0.4(SWfc?SWpwp)is the minimum soil water content for maximum photosynthesis,and SWfcis the field capacity.Available soil water is calculated from the daily precipitation,interception by leaves,above-and below-ground water runoff,and transpiration of trees (Fischer et al.2014).

    Respiration

    The respiration of a tree is the sum of its maintenance respiration,Rm,and its growth respiration,Rg,a constant fraction of(GPP-Rm).The maintenance respiration is calculated as follows:

    where Rbis a base respiration,a fraction of standing biomass of the tree(Bohn et al.2014,detailed description in supplementary information A3).κ(T)describes the influence of the daily mean air temperature T on respiration(Prentice et al.1993):

    with constants Trefand Q10(Bohn et al.2014).

    Field data and data filtering

    We compared the simulation results of the forest model with the EC data of the Wetzstein site(Table 1).For the Wetzstein site,the EC data were pre-processed as described in Rebmann et al.(2010).The net ecosystem exchange(NEE)was gap filled since the data are compared at daily time scales and partitioned into GPP and respiration.We use an algorithm that extrapolates day-time ecosystem respiration from night-time respiration considering temperature sensitivities(Reichstein et al.2005). We filtered the EC data to identify days that are affected by specific limitations.Optimal temperature or soil water conditions were defined for days when the daily mean GPP was maximal(98th percentile for the years 2003 to 2008).We assumed that on those days the GPP is not affected by any limitation.The filtered range of the optimal temperature(daily daytime mean)conditions was identified at 7.3°C 16.0%.Optimal light conditions were defined for days when values rise above the monthly 80th percentile.We define night-time as time periods when PPFD<20 μmol (CO2)?m?2?s?1.When we use normalized GPP values in our analyses,we normalize GPP values yearly by its annual 98th percentile.

    The model setup

    The forest model was run with daily time steps for three different parameterizations:literature-based(M1),numerically calibrated(M2)and filter-based(M3)parameterizations (Table 2).The literature-based parameterization(M1)is based on Bohn et al.(2014)for a spruce forest where the parameter values are derived from inventory data and the literature.The soil parameter values were adapted to theclay loam soil type as in Maidment(1993).The calibrated parameterization(M2)is based on parameters derived from a numerical calibration against the NEE,GPP and respiration data(Lehmann and Huth 2015,see Additional file 1 for details).The filter-based parameterization(M3)arose from filtering the EC data(same data as used for calibration of M2)for optimal climatic conditions(see Field data and data filtering)to isolate individual processes.Model functions were directly fitted through filtered data to derive new parameter values and a new temperature reduction curve (Eq.4).

    Table 2 Model parameter values for the literature-based(M1)with references(ref.),calibrated(M2)and filter-based(M3)model version

    All model setups were initialized according to the inventory data for Wetzstein and Tharandt(Table 1).Trees were spread equally over the 25 patches of the 1 ha model area (at Wetzstein,410 stems?ha?1with a mean stem diameter of 0.33 m).The deadwood pool was filled with 4.14 tC?ha?1(Wetzstein inventory,personalcommunication from Martina Mund,University of Goettingen).The fastdecomposing soil stock was initialized with 2.0 tC?ha?1, and the slow-decomposing soil stock with 1.5 tC?ha?1(means in the climax stage of long-term simulations). The simulation period at the Wetzstein site was from 2003 to 2008 and at Tharandt from 1999 to 2004.All model simulations were deterministic since none of the model setups included recruitment or stochastic mortality.

    Results

    Simulation of daily carbon fluxes at Wetzstein forest

    The measured and modeled(parameterization M1)GPP and respiration at Wetzstein forest for the dry year 2003 are shown in Fig.1.The forest model performed well for daily GPP.It reproduced the seasonal cycle,and the daily fluctuations showed similar magnitudes to those observed.Differences were observed during times of low soil water availability,very low and high temperatures and bright days.In late summer(July-September),the simulated respiration diverges from the observational data.In general,simulated respiration shows stronger fluctuations than the observed respiration.

    Simulation of carbon fluxes at Wetzstein forest for three different parameterizations

    We compared the simulations of all three model versions with the observed fluxes at Wetzstein at the daily time scale (Table 3,Fig.2).The correlation is best for the numerically calibrated parameterization(M2),closely followed by the filter-based parameterization(M3).Simulated respiration correlates with the observed values with an R2of 0.44–0.54; NEE with an R2of 0.63–0.66.Simulated and observed GPP values match best with an R2of 0.73–0.82.High GPP values above 30 tC?ha?1?a?1can only be reached with the filter-based parameterization(M3,Fig.2).

    In a second step,we calculated annual GPP,NEE and respiration for the Wetzstein forest(Fig.3a).The simulated annual NEE values fit the observed annual values best for the calibrated parameterization(M2, deviation from observed data by 3.5% for average values throughout the simulation years).The simulated annual GPP fits quite well and only diverges from the observed data by 0.8%for the calibrated (M2)and 2.6%for the filter-based(M3)version;simulated respiration by 0.3%(M2)and 0.9%(M3).The literature-based model version clearly underestimates all ecosystem fluxes.

    Fig.1 Time series of the daily observed and simulated data for 2003 at Wetzstein.Observed(grey line)and modeled(M1,black line)(a)GPP and (b)ecosystem respiration(R)of the ecosystem with the literature-based parameterization.Time series of the observed(c)daytime temperature (T),(d)soil water content(SW),and(e)daytime PPFD

    Simulation of carbon fluxes at another spruce forest

    We applied the forest model to another spruce forest (Tharandt)for the years 1999–2004(Fig.3b).The model parameters were the same as at Wetzstein.Only the climate input and the initialization of the forest model(forest state based on inventory data)were adapted to the site.

    Both,the filter-based(M3)and calibrated(M2)model version,performed well for the daily GPP(R2=0.61 (M2),R2=0.61(M3),RMSE=8.83 tC?ha?1?a?1(M2)and RMSE=9.57 tC?ha?1?a?1(M3),Table 3).The simulated respiration was even better than at Wetzstein,at the daily scale(R2=0.62);the NEE performed worse at Tharandt than at Wetzstein(R2=0.38).

    The annual simulated GPP is best for the calibrated (M2)and filter-based(M3)parameterizations(2 and 7% deviations).In April 2001,the tree density was reduced by 30%(at Tharandt and in the forest model),which is reflected in the lower GPP values after 2001.In the dry year 2003,the simulated annual GPP is 0.47 tC?ha?1lower than in the previous year for the calibrated modelversion(M2)and 1.37 tC?ha?1lower for the filter-based (M3).The observed GPP in 2003 is 1.87 tC?ha?1lower than in 2002.The annual carbon budget(NEE in Fig.3bmatches the observed budget well(M2:mean deviation of 20%from observations;M3:deviation of 16%). Note that NEE is the balance between GPP and respiration and the bias is consequently more sensitive.Respiration is partly overestimated for the model versions M2 and M3(deviations of 14%and 19%).

    Table 3 R2and root mean square error(RMSE,in tC?ha?1?a?1)at the calibration and validation sites

    Fig.2 Simulated vs.observed ecosystem carbon fluxes(NEE,GPP and respiration)at Wetzstein for the three parameterizations for 2192 simulated days

    Limiting factors for productivity

    As an example,we analyze here the GPP limitation due to temperature(Fig.4).Daily GPP values were normalized and filtered for optimal soil water conditions and sunny days to distinguish them from days with other limitations.Filtered,normalized data thus represent the reduction factor due to temperature.

    On colder days,the literature-based model(dashed line)shows a reduction in GPP which seemed stronger than the observed.For example,at the freezing point, the observed GPP was reduced to 20%of its optimum, whereas the model reduced the GPP to nearly 0.For higher temperatures,the observed data show a reduction of photosynthesis starting at 20°C,whereas the literature-based parameterization barely showed a reduction in photosynthesis until 30°C was reached.For the calibrated model version,the reduction curve showed a steep slope at 0°C,attained its highest values between 3°C and 20°C and reached zero at 30°C.The solid black line shows the best fit of the bell-shaped curve(Eq.4)through filtered EC-data.

    We further analyzed the relationships between the PPFD and productivity,temperature and respiration,as well as the soil water and productivity(see Additional file 1 and Table 1 for the derived filter-based parameter values).

    Discussion

    Simulating carbon fluxes at daily and yearly time scales

    We analyzed the overall model performance for three parameterizations at the daily scale(Table 3,Fig.2)and the annual scale(Fig.3).At the Wetzstein site,the forest model performed best for the numerical calibration.This performance is not surprising since the calibration method aims for the least error(Lehmann and Huth 2015).In any case,it is satisfying that we found a parameter combination that reproduces the observations with such good performance at the daily and annual time scales.This shows that the simplifications assumed in this forest model are sufficient to reproduce the complex interactions of climate and ecosystem fluxes even with the literature-based model version(M1).By comparison,LPJ-Guess simulated daily GPP with a similar performance(RMSE of 6.50–8.94 tC?ha?1?a?1and R2of 0.62–0.72)for a pine stand in the Netherlands in a similar study(the range results from different model setups, Vermeulen et al.2015).

    Fig.3 Annual ecosystem carbon fluxes(NEE,GPP and respiration) for the three parameterizations:literature-based(M1,dashed), calibrated(M2,dotted)and filter-based(M3,solid);and the observed data(circles)for(a)the Wetzstein site and(b)the Tharandt validation site

    Fig.4 Temperature limitation of photosynthesis at Wetzstein.Parameter values are listed in Table 2.The three parameterizations are the literaturebased(dashed)and the calibrated(dotted)parameterizations with the original function Eq.3),and the filter-based parameterization with the new formulation(Eq.4).The observed,daily mean temperature values of the Wetzstein site were filtered for optimal soil water and light conditions and normalized(see Methods:Observational data and data filtering)

    Simulated respiration shows a stronger bias than GPP. The deviation from the observed data can be explained in several ways.On the one hand,respiration is strongly coupled to GPP in the forest model.Consequently, strong GPP fluctuations induce strong fluctuations in respiration which is also seen at Tharandt.On the other hand,observed respiration shows a rather smooth curve (Fig.1b),which is a consequence of the flux-partitioning method that derives respiration from the measured NEE (Reichstein et al.2005).Respiration data thus arise from a modeling approach.The forest model additionally uses a simplified soil module(as in SEIB-DGVM,Sato et al. 2007).Some of the rhizosphere processes are neglected, such as the release of organic compounds by roots (Nguyen 2009).Rhizomicrobial respiration might have effects on short-term CO2efflux but should have no effect on the long-term carbon stock in the soil(Kuzyakov 2006).Also,note that different forest histories,such as different forest management at Wetzstein and Tharandt, might thus lead to different soil pools and respiration rates.The fact that the soil carbon pool of the forest model was initialized with a soil pool in a steady state might lead to an overestimation of ecosystem respiration at Tharandt.

    The simulation results of Tharandt after the dry year 2003 show that the filter-based version(M3)reproduced the drought and heat event of 2003 better than the calibrated version(M2).Hence,parameter values of individual processes of the filter-based version(M3)might be more appropriate than the calibrated ones(M2)although the overall performance of M2 is better than the performance of M3(Table 3).This might indicate that the numerical calibration compensates individual processes of extreme events to achieve an overall best performance throughout all simulation years while the filter-based version aims for the best parameterization of individual processes.

    The overall performance of all parameterizations let us conclude that carbon fluxes of forest ecosystems can also be modeled with individual-based models.The characteristic of simulating each individual tree has the advantage to investigate in plant-population dynamics, forest structure and their interaction with carbon dynamics in future studies.

    Analyzing limiting factors for productivity with the help of EC data

    To test whether the parameter values(M1,M2)correctly describe the limiting factors of photosynthesis and respiration in the forest model,we filtered the EC data. The filtering singled out time steps that are mainly influenced by a single constraint such as temperature or water stress.This approach resulted in an additional parameterization of the forest model for which the functional relations between the variable climate and model processes were directly fitted through the filtered data(M3).As an example,we here discuss the relationship between temperature and productivity.The other relationships are discussed in detail in the Additional file 1.

    The relationship between the measured GPP and temperature(Fig.4)showed that the literature-based parameterization,which was validated for long-term simulations(Bohn et al.2014),underestimated GPP for low and overestimated it for high temperatures.The function was taken from the LPJ model(Sitch et al. 2003).In LPJ,the function was applied as a“bulk”temperature-response function for different plant functional types in temperate forests to reproduce the current vegetation distribution(Stephan Sitch,personal communication).In FORMIND,this function seems to be sufficient to estimate the current mean carbon budgets from long-term simulations(Bohn et al.2014),but the model version is not able to display inter-annual changes correctly.A similar study that tested the temperature-response curve in LPJ-Guess in a pine stand concluded that the curve needs to be shifted to lower temperatures(Vermeulen et al.2015).This finding also agrees with the fact that trees can photosynthesize down to a temperature of about?5°C(root zone temperature),while growth seems to stop at approximately 0°C(K?rner and Paulsen 2004).The numerical calibration experiment also shows that shifting the curve towards lower temperatures results in a better fit to the observed data.The EC data indicate photosynthetic activity down to about?10°C.This is plausible considering that the air temperature is generally lower than temperature in the root zone during winter.In addition, the trees at the Wetzstein site are adapted to low temperatures and cold winter conditions due to the elevation.Note,that positive ecosystem GPP does not necessarily implicate tree growth since respiration compensates for GPP at low temperatures and low PPFD in the forest model.

    GPP was not limited due to high temperatures with the literature-based parameterization(M1).High temperatures might go along with a vapor pressure deficit limitation(K?rner 1994;Lasslop et al.2010),which is not considered in our model approach.This assumption is supported by the results of the numerical calibration (the limiting factor is 0 at 31.6°C).This reduction in temperature might compensate for the missing vapor pressure deficit limitation.

    Since the original reduction curve did not fit the filtered data,we introduced a new functional relationship for the temperature reduction factor(Eq.4,methods). The new curve originates from a normal distribution that describes the temperature dependence of the photosynthetic electron transport rate(June et al.2004).A completely normal distribution(n=2)led to a very small plateau in the temperature reduction function,which means a small range of optimal temperature conditions. n=4 led to a much wider range of optimal conditions. An advantage of the introduced function is that it uses only two instead of four parameters.We therefore suggest using the less complex bell-shaped curve for future studies.

    The fact that the original temperature curve could not be fitted through filtered data properly supports the assumption that a pure calibration against ecosystem fluxes does not necessarily result in optimal parameter values of the individual model processes.We can conclude that eddy covariance data and the filtering-approach can give important insights into the correct parameterization of model processes(the limiting factors).

    Sources of uncertainty

    This study comes with a variety of uncertainties from various sources that must be considered.The first source of uncertainty comes from the EC data.NEE is measured at a half-hourly scale.However,it comes with data gaps that are filled to compare the observed data with the simulated data at a daily time scale.The gap-filled NEE is based on a modeling procedure(Reichstein et al.2005).At the Tharandt site,for example,the uncertainty of the gap filling methods is up to 10%(Grünwald and Bernhofer 2007).In addition, GPP and respiration are not directly measured,but are partitioned from the NEE,which is based on another modeling procedure(Reichstein et al.2005).When we analyzed the drought event of 2003 at Tharandt,we found that GPP was 1.87 tC?ha?1less than in the previous year.A multi-site study on the event in 2003 reported a reduction of 2.08 tC?ha?1at Tharandt(Ciais et al.2005).These deviations demonstrate the uncertainties implied by gap filling and the partitioning of EC data,especially at the annual time scale. The second potential source of uncertainties comes from the filtering method and the concept of limiting factors.The forest model,and thus also the filtering method,consider only temperature and water as limiting factors.The vapor pressure deficit and its influence on productivity(as in BIOME-BGC,Kimball et al.1997; 3PG,Landsberg and Waring 1997),for example,is not considered.However,we still assume that the filtered data are reasonable for singling out different constraints in forest productivity.

    Conclusion

    The model version only based on literature values(M1, Bohn et al.2014)is capable of reproducing the seasonal cycle and daily fluctuations of carbon fluxes.However, this model version underestimates carbon fluxes at both spruce stands on the annual time scale.The calibrated model version(M2)derived from a numerical calibration (Lehmann and Huth 2015)against the observed NEE, GPP and respiration performs best at the daily andannual time scales.Deviations of the individual processes from the observed data seem to compensate each other, so that,in sum,they reproduce the observed net fluxes well.The third parameterization resulted from a fit through filtered data(M3).We identified a new functional relationship between temperature and GPP.Its mean performance at both sites differs only slightly from the calibrated parameterization,but it shows a closer match to observations for the extreme event at Tharandt in 2003. This shows that we should not blindly trust in a numerical calibration,although its overall performance is best.

    The presented filter method improved carbon flux estimates for both spruce stands by improving the model processes.The consideration of the individual limiting factors for productivity(Fig.4,Additional file 1:App.1) is essential to correctly reveal the impact of inter-annual climate variations on carbon fluxes.Therefore,we favor the filter-based model version for future studies.We can conclude that an individual-based forest model is a valuable tool that allows analyses of daily and yearly carbon fluxes in addition to the traditional analyses of forest successions and biomass.

    Additional file

    Additional file 1:Simulation results for the three parameterizations. (DOCX 826 kb)

    Abbreviations

    EC:Eddy covariance;GPP:Gross primary productions;NEE:Net ecosystem exchange;PPFD:Photosynthetic photon flux density;R:Ecosystem respiration;RMSE:Root mean square error;SW:Soil water

    Acknowledgements

    Measurements were supported by the CarboEurope-IP(European Commission, Directorate-General Research,Sixth Framework Programme,Priority 1.1.6.3:Global Change and Ecosystem(Contract No.GOCECT-2003-505572))of the Max Planck Institute for Biogeochemistry in Jena.We thank Martina Mund of the Georg-August-University of G?ttingen for the inventory information at the Wetzstein site.EC data at the Tharandt site were kindly provided by the Department of Meteorology at TU Dresden.We thank Rico Fischer,Sebastian Paulick and Franziska Taubert for their support and discussions.

    Funding

    ER,FB and AH were supported by the Helmholtz-Alliance Remote Sensing and Earth System Dynamics.ER was supported by the Helmholtz Impulse and Networking Fund through the Helmholtz Interdisciplinary Graduate School for Environmental Research(HIGRADE).

    Availability of data and materials

    Simulation results are presented as additional supporting files.The forest model FORMIND is freely available on www.formind.org.

    Authors’contributions

    ER,AH,MC worked together in designing this manuscript.FB and ER worked together on the parameterization of the forest model.CR processed the eddy covariance data.All authors were involved in the writing process.ER coordinated their contributions and led the writing.All authors read and approved the final manuscript.

    Competing interests

    The authors declare that they have no competing interests.

    Consent for publication

    Not applicable.

    Ethics approval and consent to participate

    Not applicable.

    Author details

    1UFZ-Helmholtz Centre for Environmental Research,Permoserstr.15,04318 Leipzig,Germany.2University of Osnabrück,Barbarastra?e 12,49076 Osnabrück,Germany.3German Centre for Integrative Biodiversity Research (iDiv)Halle-Jena-Leipzig,Deutscher Platz 5e,04103 Leipzig,Germany.4INRA-Université de Lorraine,UMR1137 Ecologie et Ecophysiologie Forestières,54280 Champenoux,France.

    Aubinet M,Grelle A,Ibrom A,Rannik ü,Moncrieff J,Foken T,Kowalski AS,Martin PH,Berbigier P,Bernhofer C,Clement R,Elbers J,Granier A,Grünwald T, Morgenstern K,Pilegaard K,Rebmann C,Snijders W,Valentini RM,Vesala T (1999)Estimates of the Annual Net Carbon and Water Exchange of Forests: The EUROFLUX Methodology.Adv.Ecol.Res.pp 113–175.https://doi.org/10. 1016/S0065-2504(08)60018-5

    Bohn FJ,Frank K,Huth A(2014)Of climate and its resulting tree growth:Simulating the productivity of temperate forests.Ecol Modell 278:9–17.doi:10.1016/j. ecolmodel.2014.01.021

    Botkin DB,Janak JF,Wallis JR(1972)Some ecological consequences of a computer model of forest growth.J Ecol 60:849.doi:10.2307/2258570

    Bugmann H(2001)A review of forest Gap models.Clim Change 51:259–305

    Ciais P,Reichstein M,Viovy N,Granier A,Ogée J,Allard V,Aubinet M,Buchmann N,Bernhofer C,Carrara A,Chevallier F,De Noblet N,Friend AD,Friedlingstein P,Grünwald T,Heinesch B,Keronen P,Knohl A,Krinner G,Loustau D,Manca G,Matteucci G,Miglietta F,Ourcival JM,Papale D,Pilegaard K,Rambal S, Seufert G,Soussana JF,Sanz MJ,Schulze ED,Vesala T,Valentini R(2005) Europe-wide reduction in primary productivity caused by the heat and drought in 2003.Nature 437:529–33.doi:10.1038/nature03972

    Fischer R,Armstrong A,Shugart HH,Huth A(2014)Simulating the impacts of reduced rainfall on carbon stocks and net ecosystem exchange in a tropical forest.Environ Model Softw 52:200–206.doi:10.1016/j.envsoft.2013.10.026

    Fischer R,Bohn F,Dantas de Paula M,Dislich C,Groeneveld J,Gutiérrez AG, Kazmierczak M,Knapp N,Lehmann S,Paulick S,Pütz S,R?dig E,Taubert F, K?hler P,Huth A(2016)Lessons learned from applying a forest gap model to understand ecosystem and carbon dynamics of complex tropical forests.Ecol Modell 326:124–133.doi:10.1016/j.ecolmodel.2015.11.018

    Granier A,Bréda N,Biron P,Villette S(1999)A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands.Ecol Modell 116:269–283.doi:10.1016/S0304-3800(98)00205-1

    Granier A,Reichstein M,Bréda N,Janssens IA,Falge E,Ciais P,Grünwald T, Aubinet M,Berbigier P,Bernhofer C,Buchmann N,Facini O,Grassi G, Heinesch B,Ilvesniemi H,Keronen P,Knohl A,K?stner B,Lagergren F, Lindroth A,Longdoz B,Loustau D,Mateus J,Montagnani L,Nys C,Moors E, Papale D,Peiffer M,Pilegaard K,Pita G,Pumpanen J,Rambal S,Rebmann C, Rodrigues A,Seufert G,Tenhunen J,Vesala T,Wang Q(2007)Evidence for soil water control on carbon and water dynamics in European forests during the extremely dry year:2003.Agric For Meteorol 123:123–145.doi:10.1016/j. agrformet.2006.12.004

    Grünwald T,Bernhofer C(2007)A decade of carbon,water and energy flux measurements of an old spruce forest at the Anchor Station Tharandt.Tellus B 59:387–396.doi:10.1111/j.1600-0889.2007.00259.x

    Gutiérrez AG,Armesto JJ,Aravena J-C,Carmona M,Carrasco NV,Christie DA, Pe?a M-P,Pérez C,Huth A(2009)Structural and environmental characterization of old-growth temperate rainforests of northern Chiloé Island,Chile:Regional and global relevance.For Ecol Manage 258:376–388. doi:10.1016/j.foreco.2009.03.011

    Gutiérrez AG,Armesto JJ,Díaz MF,Huth A(2012)Sensitivity of North Patagonian temperate rainforests to changes in rainfall regimes:a process-based,dynamic forest model.Biogeosciences Discuss 9:6293–6333.doi:10.5194/bgd-9-6293-2012

    Huth A,Ditzer T(2000)Simulation of the growth of a lowland Dipterocarp rain forest with FORMIX3.Ecol Modell 134:1–25.doi:10.1016/S0304-3800(00)00328-8

    June T,Evans JR,Farquhar GD(2004)A simple new equation for the reversible temperature dependence of photosynthetic electron transport:a study on soybean leaf.Funct Plant Biol 31:275–283.doi:10.1071/FP03250

    Keenan TF,Baker I,Barr A,Ciais P,Davis K,Dietze M,Dragoni D,Gough CM,Grant R, Hollinger D,Hufkens K,Poulter B,Mccaughey H,Raczka B,Ryu Y,Schaefer K, Tian H,Verbeeck H,Zhao M,Richardson AD(2012)Terrestrial biosphere model performance for inter-annual variability of land-atmosphere CO 2 exchange. Glob Chang Biol 18:1971–1987.doi:10.1111/j.1365-2486.2012.02678.x

    Kimball JS,White MA,Running SW(1997)BIOME-BGC simulations of stand hydrologic processes for BOREAS.J Geophys Res 102:29043

    K?hler P,Huth A(2004)Simulating growth dynamics in a South-East Asian rainforest threatened by recruitment shortage and tree harvesting.Clim Change 67:95–117.doi:10.1007/s10584-004-0713-9

    K?hler P,Huth A(1998)The effects of tree species grouping in tropical rainforest modelling:Simulations with the individual-based model Formind.Ecol Modell 109:301–321.doi:10.1016/S0304-3800(98)00066-0

    K?rner C(1994)Leaf Diffusive Conductances in the Major Vegetation Types of the Globe.In:Schulze E-D,Caldwell MM(eds)Ecophysiology of Photosynthesis.Springer Berlin Heidelberg,Berlin,Heidelberg

    K?rner C,Paulsen J(2004)A world-wide study of high altitude treeline temperatures.J Biogeogr 31:713–732.doi:10.1111/j.1365-2699.2003.01043.x

    Kuzyakov Y(2006)Sources of CO2 efflux from soil and review of partitioning methods.Soil Biol Biochem 38:425–448.doi:10.1016/j.soilbio.2005.08.020

    Landsberg JJ,Waring RH(1997)A generalised model of forest productivity using simplified concepts ofradiation-use efficiency,carbon balance and partltlonlng. For Ecol Manage 95:209–228.doi:10.1016/S0378-1127(97)00026-1

    Lasslop G,Reichstein M,Papale D,Richardson AD,Arneth A,Barr A,Stoy P, Wohlfahrt G(2010)Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach:critical issues and global evaluation.Glob Chang Biol 16:187–208.doi:10.1111/j.1365-2486.2009.02041.x

    Lehmann S,Huth A(2015)Fast calibration of a dynamic vegetation model with minimum observation data.Ecol Modell 301:98–105.doi:10.1016/j.ecolmodel. 2015.01.013

    Maidment D(1993)Handbook of hydrology.McGrawHill Inc,New York

    Meehl GA,Tebaldi C(2004)More intense,more frequent,and longer lasting heat waves in the 21st century.Science 305:994–997.doi:10.1126/science.1098704

    Nguyen C(2009)Rhizodeposition of Organic C by Plant:Mechanisms and Controls. In:Sustainable Agriculture.Springer Netherlands,Dordrecht,pp 97–123

    Pacala SW,Canham CD,Silander Jr JA(1993)Forest models defined by field measurements:I.The design of a northeastern forest simulator.Can J For Res 23:1980–1988.doi:10.1139/x93-249

    Piao S,Luyssaert S,Ciais P,Janssens IA,Chen A,Cao C,Fang J,Friedlingstein P, Luo Y,Wang S(2010)Forest annual carbon cost:a global-scale analysis of autotrophic respiration.Ecology 91:652–661

    Prentice IC,Sykes MT,Cramer W(1993)A simulation model for the transient effects of climate change on forest landscapes.Ecol Modell 65:51–70.doi:10. 1016/0304-3800(93)90126-D

    Rebmann C,Zeri M,Lasslop G,Mund M,Kolle O,Schulze E-D,Feigenwinter C (2010)Treatment and assessment of the CO2-exchange at a complex forest site in Thuringia,Germany.Agric For Meteorol 150:684–691.doi:10.1016/j. agrformet.2009.11.001

    Reichstein M,Falge E,Baldocchi D,Papale D,Aubinet M,Berbigier P,Bernhofer C, Buchmann N,Gilmanov T,Granier A,Grunwald T,Havrankova K,Ilvesniemi H, Janous D,Knohl A,Laurila T,Lohila A,Loustau D,Matteucci G,Meyers T, Miglietta F,Ourcival J-M,Pumpanen J,Rambal S,Rotenberg E,Sanz M, Tenhunen J,Seufert G,Vaccari F,Vesala T,Yakir D,Valentini R(2005)On the separation of net ecosystem exchange into assimilation and ecosystem respiration:review and improved algorithm.Glob Chang Biol 11:1424–1439. doi:10.1111/j.1365-2486.2005.001002.x

    Sato H,Itoh A,Kohyama T(2007)SEIB–DGVM:A new Dynamic Global Vegetation Model using a spatially explicit individual-based approach.Ecol Modell 200: 279–307.doi:10.1016/j.ecolmodel.2006.09.006

    Shugart HH(1984)A Theory of Forest Dynamics.The Blackburn Press,New Jersey

    Sitch S,Smith B,Prentice IC,Arneth A,Bondeau A,Cramer W,Kaplan JO,Levis S, Lucht W,Sykes MT,Thonicke K,Venevsky S(2003)Evaluation of ecosystem dynamics,plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model.Glob Chang Biol 9:161–185.doi:10.1046/j.1365-2486.2003.00569.x

    Smith B,Prentice IC,Sykes MT(2001)Representation of vegetation dynamics in the modelling of terrestrial ecosystems:comparing two contrasting approaches within European climate space.Glob Ecol Biogeogr 10:621–637. doi:10.1046/j.1466-822X.2001.t01-1-00256

    Sonntag M(1998)Klimaveraenderung und Waldwachstum:TREEDYN3-Simulationen mit einer Analyse modellstruktureller Unsicherheiten.PhD thesis.Universitaet Gesamthochschule Kassel.

    Thornley JHM,Johnson IR(1990)Plant and crop modelling:a mathematical approach to plant and crop physiology.Oxford University Press.https://doi. org/10.1016/0308-521X(91)90064-H

    Vermeulen MH,Kruijt B,Kabat P(2015)Modelling short term variability in carbon and water exchange in a Dutch pine forest.Earth Syst Dyn 14:5575.doi:10. 5194/esd-6-485-2015

    Zaitchik BF,Macalady AK,Bonneau LR,Smith RB(2006)Europe’s 2003 heat wave: A satellite view of impacts and land-Atmosphere feedbacks.Int J Climatol 26:743–769.doi:10.1002/joc.1280

    1 November 2016 Accepted:10 April 2017

    *Correspondence:edna.roedig@ufz.de

    1UFZ-Helmholtz Centre for Environmental Research,Permoserstr.15,04318 Leipzig,Germany

    Full list of author information is available at the end of the article

    ?The Author(s).2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/),which permits unrestricted use,distribution,and reproduction in any medium,provided you give appropriate credit to the original author(s)and the source,provide a link to the Creative Commons license,and indicate if changes were made.

    美女视频免费永久观看网站| 黄片播放在线免费| 亚洲美女黄色视频免费看| 亚洲色图 男人天堂 中文字幕 | 青春草亚洲视频在线观看| 黄片播放在线免费| 国产精品欧美亚洲77777| 黑人欧美特级aaaaaa片| 日韩欧美精品免费久久| 老司机影院成人| 在线观看www视频免费| 国产爽快片一区二区三区| 精品久久久精品久久久| √禁漫天堂资源中文www| 搡女人真爽免费视频火全软件| 欧美日韩成人在线一区二区| 亚洲美女视频黄频| 国产 精品1| 26uuu在线亚洲综合色| 丝袜人妻中文字幕| 日韩人妻精品一区2区三区| 久久久国产精品麻豆| 人妻一区二区av| 国产色婷婷99| 久久人人97超碰香蕉20202| 黄色 视频免费看| 成人二区视频| 99热全是精品| 啦啦啦啦在线视频资源| 老司机影院成人| 内地一区二区视频在线| 99热这里只有是精品在线观看| 日韩一区二区三区影片| 老司机亚洲免费影院| av女优亚洲男人天堂| 日本欧美国产在线视频| 一二三四中文在线观看免费高清| 国产综合精华液| av在线播放精品| 男人爽女人下面视频在线观看| 三级国产精品片| 99久久人妻综合| 国产色婷婷99| 日韩电影二区| 精品人妻熟女毛片av久久网站| 性高湖久久久久久久久免费观看| 国产熟女欧美一区二区| 久久久亚洲精品成人影院| 亚洲精品,欧美精品| 日本与韩国留学比较| 精品国产乱码久久久久久小说| 又大又黄又爽视频免费| 香蕉国产在线看| 99香蕉大伊视频| 日韩中文字幕视频在线看片| 免费大片18禁| 国内精品宾馆在线| 国产精品嫩草影院av在线观看| av线在线观看网站| 中文字幕制服av| 黄色 视频免费看| 又黄又爽又刺激的免费视频.| 久久精品国产亚洲av天美| 国产精品国产三级国产专区5o| 黄片播放在线免费| 观看美女的网站| 国产精品久久久久久av不卡| 人妻人人澡人人爽人人| 9热在线视频观看99| 久久99热6这里只有精品| 18禁观看日本| 巨乳人妻的诱惑在线观看| 亚洲精品一二三| 人体艺术视频欧美日本| 午夜日本视频在线| 少妇人妻久久综合中文| 爱豆传媒免费全集在线观看| 日韩欧美精品免费久久| 精品视频人人做人人爽| 亚洲精品久久久久久婷婷小说| 国产亚洲最大av| 亚洲欧美日韩另类电影网站| 97超碰精品成人国产| av在线观看视频网站免费| 精品第一国产精品| 蜜桃国产av成人99| 90打野战视频偷拍视频| 午夜激情久久久久久久| 久久久精品区二区三区| 99re6热这里在线精品视频| 男人操女人黄网站| 久久精品国产综合久久久 | 精品一区二区三区四区五区乱码 | 不卡视频在线观看欧美| 欧美日韩精品成人综合77777| 日日摸夜夜添夜夜爱| 女人被躁到高潮嗷嗷叫费观| 久久久久久久久久人人人人人人| 18禁裸乳无遮挡动漫免费视频| 一本大道久久a久久精品| 久久久久久久亚洲中文字幕| 亚洲精品视频女| 在线观看国产h片| 国产精品一二三区在线看| 99国产精品免费福利视频| 哪个播放器可以免费观看大片| 亚洲精品av麻豆狂野| 丝瓜视频免费看黄片| 亚洲成人av在线免费| 女人精品久久久久毛片| 免费黄色在线免费观看| 久久女婷五月综合色啪小说| 免费不卡的大黄色大毛片视频在线观看| 在线 av 中文字幕| 国产午夜精品一二区理论片| 国产黄频视频在线观看| videossex国产| 午夜福利视频在线观看免费| 精品国产一区二区三区久久久樱花| av在线观看视频网站免费| 秋霞伦理黄片| 亚洲欧洲精品一区二区精品久久久 | 亚洲图色成人| 午夜福利视频精品| 黑人欧美特级aaaaaa片| 午夜91福利影院| 这个男人来自地球电影免费观看 | 在线 av 中文字幕| 亚洲成人一二三区av| 亚洲精品美女久久av网站| 免费观看a级毛片全部| 国产日韩一区二区三区精品不卡| 熟女人妻精品中文字幕| 免费人成在线观看视频色| 久久99一区二区三区| 日韩中文字幕视频在线看片| 亚洲国产色片| 国产一区有黄有色的免费视频| 久久毛片免费看一区二区三区| 欧美激情极品国产一区二区三区 | 9191精品国产免费久久| 日韩制服丝袜自拍偷拍| av有码第一页| 日本av免费视频播放| 午夜久久久在线观看| 99九九在线精品视频| 久久精品国产鲁丝片午夜精品| 九九在线视频观看精品| 国产极品粉嫩免费观看在线| 国产成人av激情在线播放| 韩国高清视频一区二区三区| 午夜福利乱码中文字幕| 成人毛片a级毛片在线播放| 色婷婷av一区二区三区视频| av在线播放精品| 一二三四在线观看免费中文在 | 狠狠婷婷综合久久久久久88av| 大香蕉久久成人网| 十分钟在线观看高清视频www| 高清不卡的av网站| 日韩 亚洲 欧美在线| 一个人免费看片子| 亚洲成人av在线免费| 一级,二级,三级黄色视频| 国产成人免费观看mmmm| 男女高潮啪啪啪动态图| 国产av国产精品国产| 一级黄片播放器| 精品视频人人做人人爽| 国产精品久久久久久久电影| 欧美日韩视频精品一区| 青春草亚洲视频在线观看| 一二三四在线观看免费中文在 | 蜜桃在线观看..| 成人毛片a级毛片在线播放| 五月伊人婷婷丁香| 亚洲欧美清纯卡通| 国产亚洲精品久久久com| 99久国产av精品国产电影| 免费黄频网站在线观看国产| 欧美最新免费一区二区三区| 五月伊人婷婷丁香| av在线播放精品| 欧美3d第一页| 免费看不卡的av| 亚洲国产日韩一区二区| 国产男女超爽视频在线观看| 久久99精品国语久久久| 亚洲国产av影院在线观看| 久久ye,这里只有精品| 中文乱码字字幕精品一区二区三区| 97在线人人人人妻| 岛国毛片在线播放| 伊人亚洲综合成人网| 少妇人妻久久综合中文| 久久久久久久精品精品| 精品一区二区三区视频在线| 丝瓜视频免费看黄片| 国产精品久久久av美女十八| av有码第一页| 国产免费福利视频在线观看| 日韩精品有码人妻一区| 久久午夜综合久久蜜桃| 久久久久国产精品人妻一区二区| 久久久久久久亚洲中文字幕| 纵有疾风起免费观看全集完整版| 一级a做视频免费观看| 国产黄色视频一区二区在线观看| 亚洲欧美成人精品一区二区| 国产精品偷伦视频观看了| 青青草视频在线视频观看| 国产成人欧美| 大香蕉久久网| 女人精品久久久久毛片| 又黄又粗又硬又大视频| av在线app专区| 欧美日韩综合久久久久久| 日韩三级伦理在线观看| 99久久精品国产国产毛片| 午夜精品国产一区二区电影| 大香蕉久久网| 午夜av观看不卡| 一二三四中文在线观看免费高清| 乱码一卡2卡4卡精品| 午夜福利在线观看免费完整高清在| 亚洲综合色惰| 在现免费观看毛片| 久久久久久久国产电影| 人人妻人人澡人人爽人人夜夜| 久久精品国产自在天天线| 日韩中文字幕视频在线看片| 我要看黄色一级片免费的| 精品视频人人做人人爽| 欧美精品一区二区大全| 精品午夜福利在线看| 国产精品久久久久久av不卡| 高清在线视频一区二区三区| 亚洲经典国产精华液单| 亚洲一码二码三码区别大吗| 亚洲天堂av无毛| 中文字幕人妻熟女乱码| 一级,二级,三级黄色视频| 国产极品粉嫩免费观看在线| 欧美国产精品一级二级三级| 免费观看无遮挡的男女| 99热这里只有是精品在线观看| av在线老鸭窝| 国产熟女欧美一区二区| www.色视频.com| 美国免费a级毛片| 如何舔出高潮| 80岁老熟妇乱子伦牲交| 狠狠婷婷综合久久久久久88av| 欧美bdsm另类| 中国国产av一级| 老司机亚洲免费影院| 国产xxxxx性猛交| 国产精品久久久久久久电影| 久久久久久伊人网av| 久久99热这里只频精品6学生| 久久精品国产综合久久久 | 国产午夜精品一二区理论片| 成人亚洲精品一区在线观看| 观看av在线不卡| 免费大片18禁| 欧美成人午夜精品| 一区二区日韩欧美中文字幕 | 草草在线视频免费看| av天堂久久9| 国产免费视频播放在线视频| 欧美日韩av久久| 国产欧美日韩一区二区三区在线| 午夜影院在线不卡| 亚洲精品美女久久久久99蜜臀 | 精品国产露脸久久av麻豆| 永久免费av网站大全| 婷婷色综合大香蕉| 亚洲一级一片aⅴ在线观看| 中文字幕最新亚洲高清| 久久午夜福利片| 久久人人爽人人片av| videos熟女内射| 欧美+日韩+精品| 久久99一区二区三区| 少妇被粗大猛烈的视频| 又粗又硬又长又爽又黄的视频| 国产av国产精品国产| 亚洲欧美中文字幕日韩二区| 在线亚洲精品国产二区图片欧美| 亚洲三级黄色毛片| 亚洲精品美女久久av网站| 国产女主播在线喷水免费视频网站| 99香蕉大伊视频| 永久免费av网站大全| 如日韩欧美国产精品一区二区三区| 亚洲精品乱久久久久久| 国产综合精华液| 亚洲欧美中文字幕日韩二区| 高清毛片免费看| 国产午夜精品一二区理论片| 精品国产露脸久久av麻豆| 国产精品无大码| 啦啦啦在线观看免费高清www| 日韩成人av中文字幕在线观看| av不卡在线播放| 18禁观看日本| 免费在线观看完整版高清| 亚洲国产精品一区二区三区在线| 观看av在线不卡| 91久久精品国产一区二区三区| 夜夜骑夜夜射夜夜干| 亚洲国产欧美日韩在线播放| 久久婷婷青草| 国产有黄有色有爽视频| 午夜av观看不卡| 欧美国产精品一级二级三级| 性色avwww在线观看| 晚上一个人看的免费电影| 另类亚洲欧美激情| 天堂中文最新版在线下载| 国产国语露脸激情在线看| 久久ye,这里只有精品| 免费日韩欧美在线观看| 男人舔女人的私密视频| 宅男免费午夜| 国产成人精品无人区| 99热网站在线观看| 久久青草综合色| 亚洲精品成人av观看孕妇| 一级爰片在线观看| 黑丝袜美女国产一区| 五月伊人婷婷丁香| 久久久久人妻精品一区果冻| 色婷婷久久久亚洲欧美| 高清在线视频一区二区三区| 岛国毛片在线播放| 欧美精品一区二区免费开放| 男女边吃奶边做爰视频| 欧美人与性动交α欧美软件 | 国产极品天堂在线| 国产男人的电影天堂91| 亚洲婷婷狠狠爱综合网| 久久久久久人人人人人| 久久av网站| 亚洲一区二区三区欧美精品| 国产男女内射视频| 在现免费观看毛片| 丁香六月天网| 久久久久久久久久人人人人人人| 国产熟女午夜一区二区三区| 欧美性感艳星| 一级黄片播放器| 亚洲人成网站在线观看播放| 纵有疾风起免费观看全集完整版| 日韩大片免费观看网站| 咕卡用的链子| 亚洲图色成人| 大香蕉久久网| 欧美日韩精品成人综合77777| 在线观看国产h片| 99热全是精品| 另类亚洲欧美激情| 女性生殖器流出的白浆| 亚洲av欧美aⅴ国产| 国产国语露脸激情在线看| 18禁裸乳无遮挡动漫免费视频| 性色avwww在线观看| 97在线人人人人妻| 成人亚洲欧美一区二区av| 日本wwww免费看| 丝袜美足系列| 久久精品国产a三级三级三级| 女人久久www免费人成看片| 最黄视频免费看| 一区二区三区四区激情视频| 亚洲国产欧美日韩在线播放| 亚洲精品美女久久av网站| 亚洲天堂av无毛| xxxhd国产人妻xxx| 婷婷成人精品国产| 日本黄色日本黄色录像| 日本午夜av视频| 久久青草综合色| 国产欧美日韩一区二区三区在线| www.av在线官网国产| 99香蕉大伊视频| 日日摸夜夜添夜夜爱| 国产淫语在线视频| 香蕉国产在线看| 午夜老司机福利剧场| 在线观看免费视频网站a站| 大香蕉久久成人网| 亚洲高清免费不卡视频| 夜夜爽夜夜爽视频| 91aial.com中文字幕在线观看| 两个人看的免费小视频| 狠狠精品人妻久久久久久综合| 色网站视频免费| 国产男女内射视频| 性色av一级| 波野结衣二区三区在线| 高清av免费在线| 国产精品一二三区在线看| 伦理电影大哥的女人| 成人18禁高潮啪啪吃奶动态图| 乱码一卡2卡4卡精品| 街头女战士在线观看网站| 中文字幕精品免费在线观看视频 | 多毛熟女@视频| 水蜜桃什么品种好| 秋霞伦理黄片| 日韩人妻精品一区2区三区| 日本猛色少妇xxxxx猛交久久| 成年人午夜在线观看视频| 夜夜爽夜夜爽视频| 国产免费福利视频在线观看| 久久久久视频综合| 大香蕉97超碰在线| av黄色大香蕉| 精品熟女少妇av免费看| 国产av国产精品国产| 18禁国产床啪视频网站| 日韩制服骚丝袜av| 蜜桃在线观看..| 久久99精品国语久久久| videossex国产| 少妇被粗大的猛进出69影院 | 亚洲精品av麻豆狂野| 涩涩av久久男人的天堂| 91午夜精品亚洲一区二区三区| 男男h啪啪无遮挡| 成人午夜精彩视频在线观看| 五月伊人婷婷丁香| 91在线精品国自产拍蜜月| 国产爽快片一区二区三区| 久久精品人人爽人人爽视色| 人人妻人人澡人人看| 好男人视频免费观看在线| 涩涩av久久男人的天堂| 精品熟女少妇av免费看| 亚洲综合色网址| av女优亚洲男人天堂| 考比视频在线观看| 在线观看免费日韩欧美大片| 久久97久久精品| 日韩在线高清观看一区二区三区| 亚洲av.av天堂| 亚洲精品国产av蜜桃| 一区二区三区精品91| 少妇猛男粗大的猛烈进出视频| 久久人人97超碰香蕉20202| 日本免费在线观看一区| 亚洲精品乱久久久久久| 国产精品久久久久久精品电影小说| 欧美变态另类bdsm刘玥| 成人手机av| 久久97久久精品| 高清在线视频一区二区三区| 成年美女黄网站色视频大全免费| 午夜激情久久久久久久| 亚洲精品日本国产第一区| 日韩成人伦理影院| 极品人妻少妇av视频| 天天躁夜夜躁狠狠躁躁| 日产精品乱码卡一卡2卡三| 中文字幕人妻丝袜制服| 国产精品.久久久| 精品第一国产精品| 午夜免费男女啪啪视频观看| 啦啦啦在线观看免费高清www| 91成人精品电影| 亚洲国产精品专区欧美| 中文字幕制服av| 波多野结衣一区麻豆| 久久精品久久久久久久性| 亚洲精品久久久久久婷婷小说| 97人妻天天添夜夜摸| 天天操日日干夜夜撸| 国产老妇伦熟女老妇高清| 久久人妻熟女aⅴ| 午夜福利影视在线免费观看| 99视频精品全部免费 在线| 亚洲国产最新在线播放| 黑丝袜美女国产一区| 美女福利国产在线| www.av在线官网国产| 精品国产乱码久久久久久小说| 日本爱情动作片www.在线观看| 久久毛片免费看一区二区三区| 亚洲一级一片aⅴ在线观看| 国产日韩欧美亚洲二区| 美女主播在线视频| 成人二区视频| 在线观看免费日韩欧美大片| 亚洲五月色婷婷综合| 日韩av在线免费看完整版不卡| 亚洲欧美中文字幕日韩二区| 成人毛片60女人毛片免费| 色婷婷久久久亚洲欧美| 亚洲欧美一区二区三区黑人 | 日韩中文字幕视频在线看片| 一个人免费看片子| av又黄又爽大尺度在线免费看| 精品亚洲成a人片在线观看| av免费观看日本| 热re99久久精品国产66热6| av播播在线观看一区| 亚洲成人一二三区av| 中文欧美无线码| 亚洲精品中文字幕在线视频| 国产一区二区在线观看日韩| 纯流量卡能插随身wifi吗| 少妇人妻 视频| 精品视频人人做人人爽| xxxhd国产人妻xxx| 国产福利在线免费观看视频| 成人二区视频| 777米奇影视久久| 蜜桃在线观看..| 午夜福利视频精品| 男女边摸边吃奶| 18禁观看日本| 一本色道久久久久久精品综合| 午夜视频国产福利| 9191精品国产免费久久| 汤姆久久久久久久影院中文字幕| 国产日韩欧美亚洲二区| 欧美性感艳星| 日本-黄色视频高清免费观看| 欧美3d第一页| 777米奇影视久久| 日本黄色日本黄色录像| 曰老女人黄片| 日韩不卡一区二区三区视频在线| 一区二区日韩欧美中文字幕 | 国产伦理片在线播放av一区| 夜夜骑夜夜射夜夜干| 国产精品成人在线| 国产日韩欧美在线精品| 啦啦啦中文免费视频观看日本| 国产免费一级a男人的天堂| 晚上一个人看的免费电影| 国产国拍精品亚洲av在线观看| www.av在线官网国产| 国产精品一区www在线观看| 亚洲欧洲国产日韩| 成人午夜精彩视频在线观看| 久久午夜福利片| 超色免费av| 日本91视频免费播放| 国产成人av激情在线播放| 日本黄大片高清| 女的被弄到高潮叫床怎么办| 国产精品三级大全| 人人妻人人澡人人爽人人夜夜| 99热6这里只有精品| 一级毛片 在线播放| 欧美人与善性xxx| 在线 av 中文字幕| 人人妻人人澡人人看| 蜜桃在线观看..| 婷婷色综合www| 男人操女人黄网站| 午夜免费观看性视频| 国产免费视频播放在线视频| 亚洲av国产av综合av卡| 久久久国产欧美日韩av| 看免费av毛片| 欧美精品一区二区大全| 乱人伦中国视频| 两个人看的免费小视频| 最近最新中文字幕大全免费视频 | 国产男女内射视频| 久久女婷五月综合色啪小说| 青春草国产在线视频| 精品久久蜜臀av无| 在线天堂中文资源库| 一二三四在线观看免费中文在 | 狠狠精品人妻久久久久久综合| www.熟女人妻精品国产 | 一边摸一边做爽爽视频免费| 少妇的丰满在线观看| 国产在线免费精品| 99热全是精品| 大码成人一级视频| 久久久久久久久久久免费av| a级毛色黄片| 少妇的丰满在线观看| 大香蕉久久网| 日韩一区二区三区影片| 国内精品宾馆在线| 观看av在线不卡| 国产精品国产av在线观看| 精品久久久精品久久久| av免费观看日本| 亚洲人成网站在线观看播放| 欧美人与性动交α欧美软件 | 国产又爽黄色视频| 中文字幕av电影在线播放| 丝袜人妻中文字幕| 久久国产精品男人的天堂亚洲 | 母亲3免费完整高清在线观看 | 18禁国产床啪视频网站| 欧美国产精品va在线观看不卡| 亚洲欧美成人综合另类久久久| 国产又色又爽无遮挡免| 午夜福利影视在线免费观看| 国产一区二区在线观看日韩| 亚洲精华国产精华液的使用体验| 在现免费观看毛片| 国产麻豆69| 国产探花极品一区二区| 51国产日韩欧美| 国产精品人妻久久久影院| 90打野战视频偷拍视频|