Modeling Canopy Carbon and Water Fluxes Using a Multilayered Model over a Temperate Meadow in Inner Mongolia

To incorporate canopy vertical structure in a process-based model over a temperate meadow, a multilayered model estimated canopy carbon flux (Fc) and water flux (LE) was applied by comparing with eddy covariance measurements in Inner Mongolia, China. Simulations of diurnal, seasonal CO2 and H2O fluxes and model sensitivity to parameters and variables were analyzed. The results showed that the model underestimated Fc and LE by about 0.6% and 5.0%, respectively. It was able to simulate the diurnal and seasonal variation of Fc and LE and performed well during the day and in the growing season, but poorly at night and early in the growing season. Fc was more sensitive to the leaf nitrogen content distribution coefficient and maximum catalytic activity of Rubisco, whereas LE showed greater sensitivity to the stomatal conductance parameter a1, empirical coefficient of stomatal response to saturated vapor pressure difference Vpds0, and minimum stomatal conductance of CO2gsc0. The response of Fc to environmental factors was ranked as air CO2 concentration (Ca) > air temperature (Ta) > photosynthetically active radiation (PAR) > soil water content (θsm) > vapor pressure deficit (VPD) > wind speed (u0). The response of LE to environmental factors was ranked as Ta > VPD > θsm> PAR> Ca> u0. The response of LE to vegetation characteristic parameters was greater than that of Fc.


Introduction
Exchange of carbon and water fluxes between the atmosphere and vegetation depend on complex and non-linear interplay among physiological, ecological, biochemical, and edaphic factors and meteorological conditions (Huang et al. 2016;He et al. 2017;Fell et al. 2018). Direct observation of many variables associated with energy transfer and material circulation is difficult on a large scale. The most common models for estimation of carbon and water fluxes are statistical models (Shen et al. 2011) and comprehensive processbased coupled carbon and water models Wang and Leuning 1998). The advantages of a statistical model are that it is simple in form and easy to derive. But it is more empirical and difficult to extrapolate to other regions. The physical interpretation is unclear and a statistical model lacks predictive power. A process-based carbon and water flux coupling model has a clear mechanism. The combination of plant physiological characteristics, canopy structure, soil characteristics, and climatic variables has a certain theoretical basis and clear physical meaning. Such a model not only simulates carbon and water fluxes at the single-leaf scale (Yu et al. 1998), but can be extended to canopy (Sellers et al. 1992) and ecosystem scales Wang and Leuning 1998), and even to regional scales (McMurtrie et al. 1992) and some global climate models (Sellers et al. 1996).
Process-based coupled carbon and water flux models are based on soil-plant-vegetation-atmosphere transfer interactions and include big-leaf, two-leaf, and multilayered models. A multilayered model considers parameters measured at the leaf scale and the spatial distribution of vegetation physiological parameters as well as climatic variables, and integrates the fluxes of each layer to derive the total fluxes. Thus, a multilayered model is considered to be reasonable and has been employed in many studies Shi et al. 2010;Qu and Zhao 2016).
Mongolian grasslands are sensitive to climate change and human activities. The Mongolian steppe is located between Siberian coniferous forest and the Asian Gobi Desert. Significant changes in temperature, precipitation, and other climatic variables, and increased grassland degradation and desertification are predicted for the future (IPCC 2001). Meadow, the most productive type of grassland supporting livestock production, is located in the eastern part of Inner Mongolia (Liao and Jia 1996). Meadow steppe constitutes 11% of the 7.88 × 10 5 km 2 of native grassland in Inner Mongolia (Han et al. 2008). Approximately 80% of the plant biomass in meadow steppe is composed of Phragmites australis, Leymus chinensis, Glycine soja, Plantago asiatica, and Ranunculus microphyllus. Extensive variability in photosynthetic parameters of these species lead to differences in photosynthesis productivity (Wullschleger 1993;Chen et al. 2014), thus influencing the carbon and water exchange of grassland ecosystems. In addition, each species occupies a different height within the canopy and hence shows a different competitive ability for solar radiation (Xue et al. 2017;Zhao et al. 2018). Consequently, canopy vertical structure is more difficult to incorporate in a carbon and water fluxes model for meadow (Chang et al. 2018). The majority of previous studies focused on vegetation composed of a single species (Leuning 1995) or vegetation types other than grassland (Kucharik et al. 2006;Shi et al. 2010). To the best of our knowledge, no previous study has modeled canopy carbon and water fluxes for a mixed-species meadow.
In this paper, a multilayered model estimated canopy carbon and water fluxes was applied over a temperate meadow. It was tested by comparing simulated CO 2 and H 2 O fluxes with eddy covariance measurements across half-hour to three-year time scales firstly. And then, the sensitivity of the model to physiological, environmental, and vegetation characteristic parameters were analyzed.

Model Description
The model used in this paper was developed by Leuning et al. (1995), tested on a wheat crop canopy by Leuning et al. (1998), and revised for a mixed forest by Shi et al. (2010). The multilayered model at the canopy scale included four components: radiation absorption, coupling of a single-leaf stomatal conductance-energy balance-photosynthetic model, spatial distribution of physiological parameters, and integration of space and time.
In this paper, the canopy radiation absorption, coupling of a single-leaf stomatal conductance-energy balance-photosynthetic model, and spatial distribution of physiological parameters were similar to those of the previous studies Shi et al. 2010), whereas the integration of space and time was revised to be compatible with meadow steppe vegetation.

Canopy Radiation Absorption
Radiation absorbed by sunlit leaves (Q sl ) is expressed as  where ξ is the cumulative leaf area index measured from the top of the canopy downwards, Q lb is the diffuse radiation absorbed by all leaves at canopy depth ξ, and Q sh is the radiation absorbed by shaded leaves. Q sh is estimated by the formula: where Q ′ ld and Q lbs are the absorbed components of incoming diffuse and scattered beam radiation, respectively, and are associated with the extinction coefficient for radiation and the leaf area index (Goudriaan and van Laar 1994).

Coupling of Single-Leaf Stomatal Conductance-Energy Balance-Photosynthetic Model
Stomatal and Boundary Layer Conductance Model Total conductance of CO 2 (g tc ) and H 2 O (g tw ) from the substomatal cavities to the air outside the leaf boundary layer are calculated as: The stomatal conductance of CO 2 (g sc ) is calculated as (Leuning 1995): where A n is the CO 2 assimilation rate, g sc0 is the minimum stomatal conductance of CO 2 when A n = 0 at the light compensation point, C s is the CO 2 concentration, Γ * is the CO 2 compensation point, and the parameter a 1 is associated with the intercellular CO 2 concentration (C i ) at saturating irradiance as determined by the relationship 1/a 1 = 1 − C i / C s . The function f sm reflects the soil moisture limit (θ sm ) on stomatal conductance (Stewart and Verma 1992), Vpd s is the vapor pressure saturation deficit, and Vpd s0 is the empirical coefficient of stomatal response to saturated vapor pressure difference.
The boundary layer conductance of CO 2 (g bc ) is proportional to the boundary layer conductance of water vapor (g bw ) as: where the constant 1.37 is the ratio of diffusivity of CO 2 and water vapor in the boundary layer (Von Caemmerer and Farquhar 1981), and m is the molar volume of gas that is required to convert g bc from velocity to molar units.
The boundary layer conductance of water vapor (g bw ) is associated with its thickness and the molar diffusivity of water vapor (D w ) in air: where B t is the boundary layer thickness (Aphalo and Jarvis 1993). D w is calculated as described by Amthor (1994) from the air temperature (T a ).
The stomatal conductance of H 2 O (g sw ) is calculated as (Leuning 1995): where the constant 1.56 is the ratio of molecular diffusivity for water vapor and CO 2 in air .

Leaf Energy-Balance Model
Net radiation absorbed by the leaves is divided into sensible and latent heat. The Penman-Monteith equation is used to calculate latent heat , whereas the sensible heat is calculated by the energybalance model: where R * n is the net isothermal radiation absorbed by the leaf, Y is a coefficient associated with radiation conductance , LE and H are the latent and sensible heat exchanges between a leaf and its surroundings, respectively, ρ a is the air density, Vpd a is the vapor pressure saturation deficit in air, and c p is the specific heat capacity of dry air.

Leaf Photosynthesis Model
where min{} denotes 'minimum of', A v is ribulose bisphosphate carboxylase-oxygenase (Rubisco) activity, A j is the ribulose bisphosphate (RuP 2 ) regeneration rate allowed by electron transport, and R d is 'dark respiration' (defined as the CO 2 evolved other than through the photorespiratory pathway). R d is proportional to the maximum catalytic activity of Rubisco (V cmax ), following (Farquhar et al. 1980): Given that A n , LE and g sc are not independent, the values of C i and g sc must be determined in an iterative fashion. The model is operated in accordance with the method of Shi et al. (2010).

Spatial Distribution of Physiological Parameters
In parameterizing multilayered models, physiological and environmental characteristics should change with depth in the canopy, and sunlit and shaded leaves should be treated separately. The distribution of key environmental and plant factors at different canopy depths (ξ), wind speeds (u), and V cmax0 are estimated as follows  where u 0 is the wind speed above the top of the canopy, V cmtop is the value of V cmax at ξ = 0, ξ t is the total leaf area index, k u is the extinction coefficient for wind speed, and k n is the distribution coefficient for leaf nitrogen.

Integration of Space and Time
To apply the model to a temperate meadow composed of mixed species, we adjusted the method of dividing the canopy into multilayers. The total canopy thickness of the meadow is about 1.2 m based on the results of a field survey in the vicinity of the flux tower (Chen et al. 2014). The canopy is divided into 12 layers based on the vertical distribution of the dominant species.
The photosynthesis rate A c (t) and transpiration rate LE c (t) of each layer of the canopy of one species at time t are calculated as: where A sl and A sh are the assimilation rates of sunlit and shaded leaves, respectively, of each species, f sl and f sh are the fractions of sunlit and shaded leaf areas, LE sl and LE sh are the evapotranspiration rates of sunlit and shaded leaves of each species, [A sl f sl + A sh f sh ] is the assimilation rate of one layer, [LE sl f sl + LE sh f sh ] is the evapotranspiration rate of one layer, n is the total number of layers (12), X is the proportion of biomass in a layer of a certain species. Finally, A c (t) and LE c (t) at each point of the canopy of each species are respectively summed, which is the net photosynthesis rate A n and the evapotranspiration rate LE c of the grass canopy at that time point.

Studied Location
The experiment was carried out in a meadow of Horqin Grasslands in Horqinzuoyihou County, eastern Inner Mongolia, China (as shown in Fig. 1). A measurement tower 4 m in height was established at the study site (43°17′37″ N, 122°16′42″ E, 203 m altitude). The site experiences a temperate continental climate, cold and dry in the dormant season, and relatively warm and wet in the growing season. The mean annual precipitation is 378.8 mm, and the annual mean air temperature is 6.2 °C. Approximately 80% of the plant biomass at the study site was composed of Phragmites australis, Leymus chinensis, Glycine soja, Plantago asiatica, and Ranunculus microphyllus. The average height of the vegetation was 0.45 m in the mid-growing season from June to August. The grass at the study site is harvested for livestock forage in September. The field is flat and characterized by the Haplumbrepts soil type (USDA Taxonomy). Soil organic matter content is 8-10%, soil depth is 40-100 cm, and soil pH is 8.5-10.5.

Microclimate and CO 2 /H 2 O/Energy Flux Measurements
Eddy covariance and microclimatic instruments were installed on the meteorological tower at the experimental site and have been in operation since August 2007. Eddy covariance, routine meteorological, and solar radiation instruments were installed on the meteorological tower in the meadow. Details of the instrumentation are presented in Table 1. Original signals from the eddy covariance system were sampled at 10 Hz, whereas data from the other instruments were collected at 0.5 Hz. All of the above-mentioned meteorological data were stored as 30-min averages.

Photosynthetic Parameters of Dominant Species
To quantify photosynthetic parameters, such as V cmax , photosynthetic characteristics of the five dominant species were measured in 2010 and 2011 with a LI-6400 portable photosynthesis system (LiCor Inc., Lincoln, NE, USA) equipped with an artificial irradiance source (LiCor 6400-02B red/ blue LED light source) and a CO 2 concentration source (LiCor 6400-01 CO 2 mixer; Chen et al. 2014).

Measurement of Leaf Area Index
Leaf area index was measured by a plant canopy analyzer (LAI-2000, LiCor Inc., Lincoln, NE, USA). The observation site was selected as a cross section of about 200 m long in the southwest of the flux tower, and 12 points were selected before sunrise or after sunset, or on cloudy days. In the early (May) and end (September) growing seasons, the leaf area index is observed every 3-5 days, while every 5-10 days during the mid growing season (June-August).

Measurement of Soil Evaporation
Six small lysimeters (30 cm in height), constructed of PVC, were consisted of outer tubes with diameter of 12 cm and inner tubes with diameter of 10 cm. Lysimeters No. 1, 2, and 3 were placed under the grass canopy, and No. 4, 5, and 6 were placed in a grass-free plot nearby (area 1 m × 2 m, in which the grass was cut every few days to ensure the plot was grass-free). The outer tube was fixed in the soil so that the surface was level with the neighboring soil to avoid damage to the nearby soil structure during operation. When preparing the undisturbed soil, the bottomless lysimeter inner cylinder was vertically pressed into the soil, leaving 0.5 cm of the tube exposed above the soil surface. The inner cylinder containing the undisturbed soil was dug out, and the excess soil at the bottom of the cylinder was excised so that the bottom surface of the cylinder was flat. The bottom surface of the inner cylinder was sealed with a bottom cover with micro holes. The cylinder was weighed at 30-min intervals using a load cell (model 203A, Wuhan Lihe Electric Meter Factory, Wuhan, China). To ensure that the soil moisture content in the lysimeter was similar to that of the surrounding soil, the undisturbed soil in the lysimeter was replaced every 5-7 days.

Calibration Phase
Numerical solution method was used in the process: 1. Give g sc a small initial value, such as 0.1; 2. Simulate the stomatal conductance and boundary layer conductance under the condition of g sc , and obtain the conductance of each resistance segment; 3. Substituting the obtained conductance into the energy balance equations to calculate the leaf temperature (T l ), latent heat flux (LE) and sensible heat flux(H); 4. Substituting the calculated T l , and a given C i value into the photosynthesis model, calculate the net photosynthesis rate A n at this time. Calculate a new C i value according to the given g sc value and the relationship between A n and C i . If the difference between this C i value and the previously given C i value is less than a certain smaller critical value (such as 0.01), the model goes to the simulation of stomatal conductance, otherwise repeat the photosynthesis model, until it meets the requirements; 5. Substituting the modeled A n and C s into the stomatal conductance model to obtain a new g sc , if the difference between the value and the previously given g sc is less than a certain smaller critical value (such as 0.0001), the model goes to the next period, otherwise repeat Eqs. 2-5, until it meets the requirements.

Evaluation Phase
Model calculated F c and LE values were compared to the direct observed values using standard statistics and regression analysis. In this paper, we used Mean absolute error (MAE) and root mean squared error (RMSE) as aggregate indicators of model performance. They were computed as: where X mod,i and X obs,i are the modeled and observed F c and LE values, respectively for each day i (or 30 min) and n is the sample number.

Simulation Processes
Using the above mentioned multilayered model to simulate the canopy net photosynthesis and evapotranspiration rate from May to August (2010August ( -2011 in Horqin meadow grassland. The time step is 30 min. Using 30 min meteorological data from the meteorological observation tower in Horqin meadow grassland test site as the environmental input variable to calculate the photosynthesis within canopy in the corresponding time, thereby a series of daily and season changes in photosynthesis obtained. Furthermore, the model sensitivity to parameters and variables were analyzed. Detailed description is as follows: 1. Input of environmental factors Environmental factor data collected every half hour from May 1st to August 30th (2010-2011) (temperature T a , water vapor pressure e a , solar radiation S 0 , photosynthetically active radiation PAR, wind speed u 0 , atmospheric pressure P, atmosphere The CO 2 concentration in the C a ) input into the model. 2. Simulation of radiation absorption Based on the canopy stratification and the corresponding cumulative leaf area index ξ t , the solar radiation absorbed by the illuminated and shaded leaves in each layer and the photosynthetically active radiation are calculated. 3. Simulation of the foliar process Including photosynthesis model, energy balance equation, stomatal conductance model. Numerical solution method was used in the process as that in the calibration phase described. 4. Simulation of canopy scale The A sl , A sh , LE sl , and LE sh of the shaded leaves obtained by the foliar process simulation are substituted into the net photosynthesis calculation formula to obtain the A n and LE values of a certain layer at that time, and then return to step 2. The simulation of the other layer of the canopy, the final layers are summed, and the total photosynthesis rate A c and evapotranspiration rate LE c of the entire canopy at that time are obtained. 5. Daily values of carbon and water fluxes were calculated to analyze the seasonal variation. 6. Model sensitivity analyses when the physiological parameters, environmental factors and vegetation parameters changed.

Model Validation
The canopy CO 2 and H 2 O fluxes (A c and LE c ) used for validation were eddy covariance measurements (A eco and LE eco ) minus the soil CO 2 (R soil ) and H 2 O fluxes (LE soil ), respectively.
R soil is estimated with consideration of the soil temperature at 0.05 m soil depth (T s ) using an exponential equation: According to Chen et al. (2011), the value of a3 is 0.678 and that of b3 is 0.106.

Verification of the Multilayered Model
Observation data recorded from May to August, 2010 and 2011 in the Korqin meadow were selected to validate the simulation of carbon and water fluxes by the multilayered model. The simulated 30-min values for carbon flux F c and water flux LE and the observed ones showed good correspondence (Fig. 2). The linear regression line was close to the 1:1 line. The slope, intercept, correlation coefficient, sample number, mean absolute error (MAE), and root mean square error (RMSE) are shown in Table 2. The multilayered model underestimated F c and LE by about 0.6% and 5.0%, respectively.
The simulation results of the model were not uniform in different days, months, and years. For the different years, F c was overestimated by 8.6% and underestimated by 4.9%, whereas LE was underestimated by 3.1% and 6.6%, in 2010 and 2011, respectively. For different months, F c was underestimated by 18.8% and overestimated by 1.9%, and LE was underestimated by 34.3% and 1.8%, early in the growing season (May) and during the growing season (June to August), respectively. For daytime (photosynthetically active radiation PAR > 5 μmol m −2 s −1 ), F c and LE were overestimated by 2.0% and underestimated by 5.8%, respectively.

Simulation of Diurnal Variation of Carbon and Water Fluxes
Simulation of diurnal variation of carbon and water fluxes was analyzed using observation data for three consecutive days per month from May to August in 2011 (Fig. 3). The CO 2 and H 2 O fluxes showed distinct diurnal cycles, represented by bell-shaped curves, with low and steady values at night and unimodal distribution of values during the day. The carbon flux was positive at night and negative in the daytime. The ecosystem changed from CO 2 release to absorption after sunrise, and reverted from CO 2 absorption to release before sunset. The maximum net absorption of CO 2 was attained at about 11:30. The water vapor flux was small at night, increased after sunrise, and gradually decreased after sunset. The maximum water vapor flux was attained at around 12:30. Trends for the simulated and observed 30-min values of F c and LE were essentially identical. Correspondence between the simulated and observed values was superior during the day than at night.

Simulation of Seasonal Variation of Carbon and Water Fluxes
The seasonal variation (Fig. 4) showed that the simulation values for F c were consistent with the observed values. The total F c early in the growing season (May) was low, whereas in the growing season (June-August) was much higher. These results are associated with the LAI, leaf physiological activity, and environmental factors such as temperature and radiation. Compared with the observed eddy-covariance values, the simulated values for F c were higher early in the growing season (May) because leaf physiological activities are poor and the chlorophyll content is lower. The simulation of LE was not reliable, the simulated value was low early in the growing season (May) because the LAI is lower in May and soil evaporation is a large component of LE.

Sensitivity of Multilayered Model to Parameters and Variables
We analyzed sensitivity of multilayered model to the physiological parameters (the parameter a 1 , the empirical coefficient of stomatal response to saturated vapor pressure difference Vpd s0 , the minimum stomatal conductance of CO 2 at the light compensation point g sc0 , the leaf nitrogen content distribution coefficient k n , and maximum Rubisco catalytic activity V cmax ), environmental factors (PAR, vapor pressure deficit VPD, wind speed u 0 , soil water content θ sm , air temperature T a , atmospheric CO 2 concentration C a ), and vegetation parameters (leaf area index LAI). The cumulative canopy carbon and water fluxes from May to August with changes in the parameters were shown in Fig. 5. The amplitude of canopy cumulative F c and LE after changes of parameters is shown in Table 3. The cumulative canopy F c and LE declined with decrease in the value of a 1 , Vpd s0 , g sc0 , PAR and vice versa. But they increased when k n , V cmax decreased and vice versa. The cumulative F c and LE decreased concomitant with decrease in temperature. With progressive increase in temperature, the cumulative F c and LE initially increased but thereafter decreased. VPD had little effect on F c , but had a greater influence on LE. With a reduction in VPD, cumulative F c increased and LE decreased, and vice versa. The effect of u 0 on cumulative F c and LE was small; when the wind speed increased by 10%, the influence on simulated values of F c and LE was less than 0.1%. C a had a greater effect on F c and less impact on LE. The effect of θ sm on F c and LE was small.
CO 2 flux was most sensitive to k n and V cmax , whereas water vapor flux was most sensitive to other physiological parameters. The sensitivity of F c to individual environmental factors was ranked as C a > T a > PAR > θ sm > VPD > u 0 , whereas the rank order of sensitivity of LE to environmental factors was T a > VPD > θ sm > PAR > C a > u 0 . The response of LE to vegetation parameters was greater than that of F c . Overall, the multilayered model could be used to simulate the multi-species meadow canopy productivity under climate change.

Common Issues When Simulating Carbon and Water Flux
The multilayered model underestimated F c and LE by about 0.6% and 5.0%, respectively at meadow site. Similar phenomenon found by some scholars (Shi et al. 2010;Chang et al. 2018;Reed et al. 2018). Underestimations of F c at meadow site are expected to have resulted from soil respiration which was simply to the exponential relationship of soil temperature. Underestimations of LE are expected to have resulted from canopy evaporation, transpiration, and soil evaporation estimation errors. Systematic measurements on each component of evapotranspiration are essential for the further evaluation and calibration of latent heat flux estimations for vegetated surfaces (Zhang et al. 2016).
The model performed well during the day and in the growing season, but poorly at night and early in the growing season. The reason for this difference is that F c during the day mainly includes canopy photosynthesis, and vegetation and soil respiration, and canopy photosynthesis activity is much greater than respiration activity. The simulated values for the canopy were similar to the observed values. The F c at night is mainly from vegetation and soil respiration. Soil respiration is affected by soil moisture content, nutrient concentrations, and soil microbe activities (Huang et al. 2015;Drollinger et al. 2017). The current model does not consider the latter two factors. The low wind speed led to a lower air temperature at the soil surface than at the reference height, resulting in an atmospheric inversion layer, but the model does not account for this temperature difference. The reason for poor estimation early in the growing season is that physiological and ecological processes in the canopy leaves were mainly considered, but LAI and V cmax at the beginning and end of the growing season each differ substantially, and hence resulted in the deviations in the simulations. Mover, leaf chlorophyll concentration was not measured owing to the limitation of the installed instruments, therefore the impact of chlorophyll content is ignored and the F c value is small in May. Thus, the percentage error is relatively large. The reason for the dispersion on LE simulated data is that soil evaporation is a large component of LE. It is not only related to LAI, but also associated with solar radiation, temperature, the saturated vapor pressure difference, wind speed, and other environmental factors. In addition, relatively few leaves have developed in the vegetation in May, and thus the model mainly simulated canopy evaporation and neglected that of the stem. Overall, the multilayered model could be used to simulate the multi-species meadow canopy productivity.

Values of the Parameters of the Model
Simulation of CO 2 and H 2 O fluxes affected by values of key input parameters and variables of the model (Liu et al. . The parameter a 1 is associated with the intercellular CO 2 concentration at light saturation. Leuning et al. (1995) determined in the multilayered model that a 1 ranged in value from 4 to 9 depending on the leaf nitrogen content. Wang and Leuning (1998) reported the a 1 value to be 11 in a two-leaf model. The empirical coefficient Vpd s0 reflects the sensitivity of leaf stomatal conductance to the surface saturated vapor pressure difference Vpd s . Leuning et al. (1995) determined the Vpd s0 value to be 1.5-3.5 kPa. Harley et al. (1992) reported the value of g sc0 to be 0.03-0.05 mol Fig. 5 Sensitivity analysis of the multi-layered model for the parameters (take 2011 for example). X-axis was the date, and Y-axis was the canopy cumulative of carbon and water fluxes before and after changes of parameters. The black solid lines were the initial values CO 2 m −2 s −1 . The parameter k n reflects the degree of attenuation of leaf nitrogen in the canopy. The majority of leaf nitrogen is a constituent of photosynthesis-related enzymes.
The value of k n varies among studies, and between sparse and dense canopies, with reported values ranging from 0.08 to 2.8 (Schieving et al. 1992;Archontoulis et al. 2011).  Leuning et al. (1995) showed that change in k n (0.0-1.2) had little impact on the CO 2 assimilation rate under any light condition if the leaf nitrogen content was high. However, if the leaf nitrogen content was low, the total daily CO 2 assimilation increased by about 10-16% if k n increased from 0.0 to 1.2. The maximum Rubisco catalytic response rate V cmax varies from 23.4 to 99.6 μmol m −2 s −1 among species and observation seasons (Chen et al. 2014). In this paper, we evaluated a1 as 5, Vpd s0 as 1.5 kPa, g sc0 as 0.03 mol m −2 s −1 , k n as 2.5, and V cmax as 23.4-99.6 μmol m −2 s −1 separately among five species. We further analyzed the sensitivity of multilayered model to parameters and variables, the results indicated that F c was more sensitive to k n and V cmax , whereas LE showed greater sensitivity to other physiological parameters (a 1 , V pds0 and g sc0 ).