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

Abstract

To incorporate canopy vertical structure in a process-based model over a temperate meadow, a multilayered model estimated canopy carbon flux (F_c) and water flux (LE) was applied by comparing with eddy covariance measurements in Inner Mongolia, China. Simulations of diurnal, seasonal CO₂ and H₂O fluxes and model sensitivity to parameters and variables were analyzed. The results showed that the model underestimated F_c and LE by about 0.6% and 5.0%, respectively. It was able to simulate the diurnal and seasonal variation of F_c and LE and performed well during the day and in the growing season, but poorly at night and early in the growing season. F_c 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 a₁, empirical coefficient of stomatal response to saturated vapor pressure difference V_pds0, and minimum stomatal conductance of CO₂g_sc0. The response of F_c to environmental factors was ranked as air CO₂ concentration (C_a) > air temperature (T_a) > photosynthetically active radiation (PAR) > soil water content (θ_sm) > vapor pressure deficit (VPD) > wind speed (u₀). The response of LE to environmental factors was ranked as T_a > VPD > θ_sm> PAR> C_a> u₀. The response of LE to vegetation characteristic parameters was greater than that of F_c.

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 process-based coupled carbon and water models (Leuning et al. 1995; 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 (Leuning et al. 1995; 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 (Leuning et al. 1995; 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⁵ km² 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₂ and H₂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.

Materials and Methods

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 (Leuning et al. 1995; 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 (Leuning et al. 1995):

$$Q_{sl} (\xi ) = Q_{lb} + Q_{sh} (\xi )$$

(1)

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:

$$Q_{sh} (\xi ) = Q_{ld}^{{\prime }} (\xi ) + Q_{lbs} (\xi )$$

(2)

where $Q_{\text{ld}}^{{\prime }}$ 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₂ (g_tc) and H₂O (g_tw) from the substomatal cavities to the air outside the leaf boundary layer are calculated as:

$$g_{tc} = \frac{1}{{{1 \mathord{\left/ {\vphantom {1 {g_{sc} + {1 \mathord{\left/ {\vphantom {1 {g_{bc} }}} \right. \kern-0pt} {g_{bc} }}}}} \right. \kern-0pt} {g_{sc} + {1 \mathord{\left/ {\vphantom {1 {g_{bc} }}} \right. \kern-0pt} {g_{bc} }}}}}}$$

(3)

$$g_{tw} = \frac{1}{{{1 \mathord{\left/ {\vphantom {1 {g_{sw} }}} \right. \kern-0pt} {g_{sw} }} + {1 \mathord{\left/ {\vphantom {1 {g_{bw} }}} \right. \kern-0pt} {g_{bw} }}}}$$

(4)

The stomatal conductance of CO₂ (g_sc) is calculated as (Leuning 1995):

$$g_{sc} = g_{sc0} + \frac{{a_{1} \, \cdot \,f_{sm} \, \cdot \,A_{n} }}{{(C_{s} - \varGamma^{*} )\, \cdot \,(1\, + \,Vpd_{s} /Vpd_{s0} )}}$$

(5)

where A_n is the CO₂ assimilation rate, g_sc0 is the minimum stomatal conductance of CO₂ when A_n = 0 at the light compensation point, C_s is the CO₂ concentration, Γ^* is the CO₂ compensation point, and the parameter a₁ is associated with the intercellular CO₂ concentration (C_i) at saturating irradiance as determined by the relationship 1/a₁ = 1 − C_i/C_s (Leuning et al. 1995). 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₂ (g_bc) is proportional to the boundary layer conductance of water vapor (g_bw) as:

$$g_{bc} = {{g_{bw} } \mathord{\left/ {\vphantom {{g_{bw} } {\left( {1.37\, \cdot \,m} \right)}}} \right. \kern-0pt} {\left( {1.37\, \cdot \,m} \right)}}$$

(6)

where the constant 1.37 is the ratio of diffusivity of CO₂ 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:

$$g_{bw} = {{D_{w} } \mathord{\left/ {\vphantom {{D_{w} } {B_{t} }}} \right. \kern-0pt} {B_{t} }}$$

(7)

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₂O (g_sw) is calculated as (Leuning 1995):

$$g_{sw} = 1.56\, \cdot \,g_{sc} \, \cdot \,m$$

(8)

where the constant 1.56 is the ratio of molecular diffusivity for water vapor and CO₂ in air (Leuning et al. 1995).

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 (Leuning et al. 1995), whereas the sensible heat is calculated by the energy-balance model:

$$LE = \frac{{s\, \cdot \,Y\, \cdot \,R_{n}^{*} + \rho_{a} \, \cdot \,c_{p} \, \cdot \,Vpd_{a} \, \cdot \,g_{bh} }}{{s\, \cdot \,Y + \gamma \, \cdot \,{{g_{bh} } \mathord{\left/ {\vphantom {{g_{bh} } {g_{tw} }}} \right. \kern-0pt} {g_{tw} }}}}$$

(9)

$$H = Y\, \cdot \,\left( {R_{n}^{*} \, - \,LE} \right)$$

(10)

where $R_{n}^{*}$ is the net isothermal radiation absorbed by the leaf, Y is a coefficient associated with radiation conductance (Leuning et al. 1995), 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

Farquhar et al. (1980) presented a biochemical model of photosynthetic carbon assimilation in C₃ plants:

$$A_{n} = \hbox{min} \left\{ {A_{v} ,A_{j} } \right\} - R_{d}$$

(11)

where min{} denotes ‘minimum of’, A_v is ribulose bisphosphate carboxylase-oxygenase (Rubisco) activity, A_j is the ribulose bisphosphate (RuP₂) regeneration rate allowed by electron transport, and R_d is ‘dark respiration’ (defined as the CO₂ 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):

$$R_{d} = 0.0089 \cdot V_{c\hbox{max} }$$

(12)

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 (Leuning et al. 1995):

$$u(\xi ) = u_{0} \, \cdot \,\exp \left( { - k_{u} \, \cdot \,\xi } \right)$$

(13)

$$V_{c\hbox{max} } = V_{cmtop} \, \cdot \,\left( { - k_{n} \, \cdot \,\xi /\xi_{t} } \right)$$

(14)

where u₀ 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:

$$A_{c} (t) = \sum\limits_{i = 1}^{n} {[A_{sl} } f_{sl} + A_{sh} f_{sh} ]X$$

(15)

$$LE_{c} (t) = \sum\limits_{i = 1}^{n} {\left[ {LE_{sl} f_{sl} + LE_{sh} f_{sh} } \right]} X$$

(16)

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_slf_sl + A_shf_sh] is the assimilation rate of one layer, [LE_slf_sl + LE_shf_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.

Fig. 1

Map of site’s location

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.

Observed Data Collection

Microclimate and CO₂/H₂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.

Table 1 Introduction of eddy covariance, routine meteorological and solar radiation measurements

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₂ concentration source (LiCor 6400-01 CO₂ 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:

$${\text{MAE}} = \frac{1}{{n}}\sum\limits_{{{{i}} = 1}}^{{n}} {\left| {X_{{\bmod ,{{i}}}} - X_{\text{obs},i} } \right|}$$

(17)

$${\text{RMSE}} = \sqrt {\frac{1}{{n}}\sum\limits_{{{{i}} = 1}}^{{n}} {\left( {X_{{\bmod ,{{i}}}} - X_{\text{obs},i} } \right)^{2} } }$$

(18)

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 (2010–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₀, photosynthetically active radiation PAR, wind speed u₀, atmospheric pressure P, atmosphere The CO₂ 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₂ and H₂O fluxes (A_c and LE_c) used for validation were eddy covariance measurements (A_ecoand LE_eco) minus the soil CO₂ (R_soil) and H₂O fluxes (LE_soil), respectively.

$$A_{eco} = A_{c} + R_{soil}$$

(19)

$$LE_{eco} = LE_{c} + LE_{soil}$$

(20)

R_soil is estimated with consideration of the soil temperature at 0.05 m soil depth (T_s) using an exponential equation:

$$R_{soil} = a3\, \cdot \,\exp \left( {b3\, \cdot \,T_{s} } \right)$$

(21)

According to Chen et al. (2011), the value of a3 is 0.678 and that of b3 is 0.106.

Results and Discussion

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.

Fig. 2

Comparison of modeled (Y-axis) against measured (X-axis) 30 min carbon and water flux a carbon flux (F_c), and b water flux (LE) during May to August, 2010–2011. The solid lines were the best fit linear regressions (y = ax + b), and the dashed lines represented 1:1 relation

Table 2 Linear regressions between measured and modeled 30 min carbon and water fluxes and model error estimation

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⁻² s⁻¹), 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₂ and H₂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₂ release to absorption after sunrise, and reverted from CO₂ absorption to release before sunset. The maximum net absorption of CO₂ 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.

Fig. 3

Comparison of diurnal variation of modeled (line) and measured (dot) canopy carbon and water flux during growing season (May–August, 2011) aF_c, bLE. 3 days’ data in each month were randomly selected

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.

Fig. 4

Seasonal variation of modeled and measured canopy daily carbon and water fluxes for May to August, 2010–2011

Sensitivity of Multilayered Model to Parameters and Variables

We analyzed sensitivity of multilayered model to the physiological parameters (the parameter a₁, the empirical coefficient of stomatal response to saturated vapor pressure difference Vpd_s0, the minimum stomatal conductance of CO₂ 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₀, soil water content θ_sm, air temperature T_a, atmospheric CO₂ 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₁, 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₀ 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.

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

Table 3 Amplitude of canopy cumulative carbon and water fluxes after changes of parameters

CO₂ 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₀, whereas the rank order of sensitivity of LE to environmental factors was T_a > VPD > θ_sm > PAR > C_a > u₀. 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.

Discussion

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₂ and H₂O fluxes affected by values of key input parameters and variables of the model (Liu et al. 2015). The parameter a₁ is associated with the intercellular CO₂ concentration at light saturation. Leuning et al. (1995) determined in the multilayered model that a₁ ranged in value from 4 to 9 depending on the leaf nitrogen content. Wang and Leuning (1998) reported the a₁ 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 CO₂ m⁻² s⁻¹. 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₂ 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₂ 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⁻² s⁻¹ 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⁻² s⁻¹, k_n as 2.5, and V_cmax as 23.4–99.6 μmol m⁻² s⁻¹ 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₁, V_pds0 and g_sc0).