Leaf 13 C data constrain the uncertainty of the carbon dynamics of temperate forest ecosystems

Stable carbon isotope discrimination occurred in plant biophysical and biogeochemical processes can help understand plant physiology and soil biogeochemistry with respect to carbon cycling. Here, we incorporated the photosynthetic carbon isotope discrimination into a process-based land surface model (iTEM) to test if stable carbon isotope composition (δC) can impose additional constraint on model parameters. Sequential data assimilation was implemented at six eddy covariance flux tower sites using carbon flux observations including gross primary productivity (GPP) and net ecosystem exchange (NEE) with and without considering foliar δC (δCf) measurement constraints, respectively. Our model-data fusion showed that δCf can provide useful constraint on photosynthetic (Vcmax25, the maximum rate of carboxylation at 25°C) and stomatal (g1, the slope of stomatal function) parameters as well as the posterior carbon fluxes. When including δCf measurement, g1 spatially varies among the six sites and is significantly correlated with annual precipitation. We incorporated the statistical relationship between g1 and annual precipitation into iTEM, which is then used to quantify the regional carbon dynamic in temperate forest ecosystems of the Northern Hemisphere. Compared with the simulation only conditioned on carbon flux observations, regional carbon flux estimations performed slightly better against the FLUXCOM products and the uncertainties of modeled carbon fluxes were reduced by 27%. Our study demonstrated that δCf data constrains carbon flux uncertainties across space.


INTRODUCTION
Quantification of ecosystem carbon fluxes and their responses to environmental changes would improve our understanding of carbon-climate feedbacks. Land surface models (LSM) incorporated with key biophysical and biochemical processes, including photosynthesis, respiration, and evapotranspiration, are important tools to predicting carbon exchanges between terrestrial biosphere and the atmosphere (IPCC 2014). To date, the model-data fusion approach encompassing model parameter optimization and uncertainty quantification using eddy flux tower observation (e.g., carbon and water fluxes) has been widely used for quantifying ecosystem responses to changing climate (Williams et al. 2009).
As an important tracer in the biogeochemical process, the stable carbon isotope composition (designated as δ 13 C) provides novel insights on plant physiological, soil biological, physical, and chemical processes (Brüggemann et al. 2011) and is potentially useful to constrain the processbased ecosystem models. The carbon isotope composition in atmospheric CO 2 (designed as δ 13 C atm ) has been used to provide an additional constraint on carbon budget quantifications (Tans et al. 1993, Ciais et al. 1995 and to examine the relative contribution of C3 and C4 plants to total primary productivity (Sage et al. 1999, Griffis et al. 2010. At the local scale, the δ 13 C measurement has been adopted in ecological research as integrators of how organisms interact with and respond to the abiotic and biotic environment changes. For example, the δ 13 C variation of plant tissue can improve our understanding on the response of plants to changes in precipitation, vapor pressure deficit (Farquhar and Richards 1984, Bowling et al. 2002, Diefendorf et al. 2010) and the intrinsic water use efficiency (Saurer et al. 2004, Seibt et al. 2008). In addition, the high-frequency in situ carbon isotopic flux measurement provides a novel tool to partitioning the net ecosystem exchange (NEE) into photosynthesis and respiration (Bowling et al. 2001).
Abundant in situ foliar δ 13 C (designated as δ 13 C f ) measurements have been published and the relationship between the δ 13 C f variation and environmental conditions (Marshall et al. 2007), edaphic factors (Arens et al. 2000) and plant attributes (Kaplan et al. 2002) has been well documented. A theoretical treatment of carbon isotope discrimination (the offset between δ 13 C atm and δ 13 C f , designated as Δ) in C3 and C4 vegetation photosynthesis (Farquhar and Richards 1984) has been incorporated into modeling to estimate plant carbon discrimination at different spatial and temporal scales (Fung et al. 1997, Suits et al. 2005, Scholze et al. 2008. As for the 13 C fractionation process during photosynthesis, the photosynthetic fractionation in C3 plant is modeled as: where a (2.9‰), b (4.4‰), c (1.8‰), and d (28.2‰) are constants representing fractionation due to diffusion at the boundary layer, stomatal cavity, CO 2 entering solution and Rubisco CO 2 fixation, respectively. C a , C s , C i , and C c are the CO 2 concentration in the air, at the boundary layer, stomata, and chloroplast, respectively. In literature, a simplified two-stage (diffusion and carboxylation) model is often applied to estimate photosynthetic carbon isotope fractionation: where A (4.4‰) and B (28‰) are constants representing fractionation due to diffusion and carboxylation, respectively. According to Farquhar photosynthesis formula (Farquhar et al. 1980), the carbon assimilation rate is a function of the intercellular CO 2 concentration (C c or C i ), and therefore, the in situ observation of carbon fluxes combined with the measurements of isotopic signal of plant tissues, which integrate all the biophysical and biogeochemical processes, may have good constraints on the parameters related to biophysical and biochemical processes at the canopy scale. Previous studies have attempted to combine in situ δ 13 C f measurement for calibrating LSMs (Aranibar et al. 2006;Raczka et al. 2016;Duarte et al. 2017 This study is to assess if δ 13 C f measurement can constrain the key parameters related to the carbon fluxes at multiple forest sites for a process-based LSM (iTEM; Chen and Zhuang 2014, Liu et al. 2014. By implementing a rigorous model-data fusion approach using observation of both carbon flux and δ 13 C f observations, we aimed at addressing two questions: (1) Can δ 13 C f measurement consistently improve model parameterization across multiple sites? (2) What are the regional implications to carbon fluxes when δ 13 C f is included as an additional constraint to model parameters?

Land surface model (iTEM) and photosynthetic fractionation process
In the integrated Terrestrial ecosystem model (iTEM), the canopy is modeled with a one-layer, two-big-leaf approach (Dai et al. 2004), which diagnoses energy budget, leaf temperature, evapotranspiration, and photosynthesis separately for v www.esajournals.org sunlit and shaded leaves. The boundary layer turbulent processes are modeled based on the Monin-Obukhov similarity theory. The hydrological processes include the interception, through fall of precipitation, snow accumulation, sublimation and melt, surface runoff, surface evapotranspiration, water infiltration, and redistribution in soil and subsurface drainage. These algorithms allow the model to simulate the response of land surface processes to changing direct and diffusive radiation regimes, such as surface energy balance, thermal dynamics, leaf and canopy conductance, and surface evapotranspiration. In addition, the iTEM kept the approach of TEM in modeling ecosystem carbon and nitrogen dynamics and their interactions . Two carbon pools and four nitrogen pools are used to represent carbon and nitrogen pools in a terrestrial ecosystem, including vegetation carbon, soil organic carbon, vegetation labile nitrogen, vegetation structural nitrogen, soil organic nitrogen, and soil inorganic nitrogen. More details of the iTEM are documented in Chen (2013).
As for the 13 C discrimination process during photosynthesis, the 13 C fractionation during photosynthesis is modeled as: where a (2.9‰), b (4.4‰), c (1.8‰), and d (28.2‰) were constants representing fractionation due to diffusion at the boundary layer, stomatal cavity, CO 2 entering solution and Rubisco CO 2 fixation, respectively. g b is the boundary layer conductance, and C a , C s , C i , and C c represented the CO 2 concentration in the air, at the boundary layer, stomata, and chloroplast, respectively. We followed Evans and von Caemmerer (1996), assuming mesophyll conductance (g m ) is proportional (4800) to the Rubisco content (V m ). δ 13 C atm was assumed as constant (−7.8‰) (White et al. 2015). It was not clear how the postphotosynthetic fractionation and leaf carbon allocation would affect δ 13 C f , here we assumed that δ 13 C f equals the isotopic composition of assimilated carbon (δ 13 C A ) (Aranibar et al. 2006). More detailed steps about the incorporation of δ 13 C f into iTEM can be found in the Appendix S1.

Data
We selected six temperate forest sites from the European Flux network (Table 1, http://www. europe-fluxdata.eu/home). The in situ observed vapor pressure, air temperature, incoming radiation, wind speed, and air pressure were used to drive iTEM. The daily time step net ecosystem exchange (NEE) measurement using the eddy covariance method were from the Level 4 dataset. The NEE data were previously frictionvelocity-filtered, gap-filled, and partitioned into the component fluxes of gross primary productivity (GPP) and ecosystem respiration (Reco) using an online tool (http://www.bgc-jena.mpg.de/ MDIwork/eddyproc/, Reichstein et al. (2005)).
Here we used the gap-filled NEE values (and derived GPP) based on the marginal distribution sampling method for model-data fusion experiment.
The foliar 13 C compositions (δ 13 C f ) at the six sites were obtained from the Work Package 5 (WP5, entitled "Isotopic Studies") of v www.esajournals.org CARBOEUROFLUX project (http://www. weizmann.ac.il/EPS/wp5/results.html). One of the WP5 goals is to monitor and elucidate the ecosystem-scale processes by analyzing the carbon isotope signal in organic tissues as indicators for ecophysiological response to environmental change. The foliar samples were generally taken from the dominant species at the sites with 2-3 fully expanded and sunlight leaves at the top canopy (https://www.weizmann.ac.il/EPS/wp5/ sample.html). Compared with the continuous measurement of carbon fluxes (NEE) at these sites, the δ 13 C f measurements were conducted every one or two months during growing season. The observed δ 13 C f values of sunlit leaves reflected the integral of photosynthetic discrimination, expressed by the Eq. 3 above. For the regional simulations, model runs were conducted at a 3-hourly time step for the period 2003-2010 and at a spatial resolution of 1°×1°for the temperate forest region in the Northern Hemisphere. The meteorological data of air temperature, wind speed, radiation, CO 2 , precipitation, water vapor concentration and surface air pressure, the initial conditions, soil properties, and the vegetation distribution were from Chen and Zhuang (2014).

Data assimilation and experiment setup
Our model-data fusion framework employed the ensemble Kalman Filter (EnKF) (Evensen 1994) with an augmented state formulation that allows model parameters to be estimated (Evensen 2009). EnKF incorporates an ensemble of forecasts to estimate background error covariance, which is a Monte Carlo approximation to the traditional Kalman filter. Then, an analysis step improves the forecast by calculating an updated ensemble of model vectors based on available observations: K is the Kalman gain. x f t and x a t are both augmented model state vectors that comprise model states (S) and parameters (p) for each ensemble (i) before and after the analysis step, respectively. y t is the observation vector, and H is the operator converting the model states to observation space. P f t is the background error covariance matrix, and R is the observation error covariance matrix. Five vegetation related parameters in iTEM were selected for estimation at our selected sites (Table 2). These parameters were found to be sensitive in Chen and Zhuang (2014). The prior parameter distributions were assumed to be uniformly distributed with the range of values reported in White et al. (2000). The EnKF implementation then used two types of data (y) for the analysis step, gross primary productivity (GPP) and net ecosystem exchange (NEE) (E0 experiment). To assess if δ 13 C f can further constrain iTEM parameters, another EnKF experiment was implemented using both GPP, NEE and available δ 13 C f observation at these sites (E1 experiment). Based on preliminary EnKF results, an ensemble size of 100 was chosen as a balance of performance and computational feasibility. The data assimilation period started from 2001 to 2002 for all six sites and the conditioned parameters at the final time step for the two EnKF implementations were used for posterior site-level simulations and regional extrapolations. Since EnKF had to forecast the model output for 100 times at each time step, we explored the OPENMPI tool to write the above EnKF algorithm as parallel mode to reduce

RESULTS AND DISCUSSIONS
Data-conditioned parameters and posterior sitelevel simulations Results in Fig. 1 show the data-conditioned parameters at each site with the two model-data fusion experiments. Both E0 and E1 had similar estimates in the nitrogen related parameters (N fall , N max and N up ), but with different V max25 and g1 values. In addition, the uncertainties of V max25 and g1 were narrowed to smaller ranges at all sites when δ 13 C f observations were assimilated in iTEM. For example, the V max25 of DeTha site was estimated to range from 60 to 77 μmol CO 2 Ám −2 Ás −1 when only GPP and NEE were assimilated into iTEM ( Fig. 1; E0 experiment), while its uncertainty can be reduced with additional δ 13 C f observation constraint (E1 experiment), ranging from 63 to 71 μmol CO 2 Ám −2 Ás −1 . However, the uncertainties of nitrogen parameters (N fall , N max , and N up ) did not change significantly, this may be because: (1) The three parameters were relatively well constrained in E0; and/or (2) their variability had negligible effects on δ 13 C f simulations.
Figs. 2, 3, respectively, show the posterior simulated GPP and NEE using the data-conditioned parameters from E0 and E1 experiment. Overall, both experiments can capture the in situ carbon Fig. 1. Comparison of data-conditioned parameters from the two model-data fusion experiments (E0 and E1) at six forest sites. E0 used GPP and NEE observations as constraint, while E1 used GPP, NEE and δ 13 C f observations as constraint. More details on the parameter description in the figure refer to Table 2. v www.esajournals.org fluxes variation reasonably well, with the adjusted r 2 ranging from 0.4 to 0.74 (Appendix S1: Fig. S1). However, if only GPP and NEE were assimilated into iTEM (E0 experiment), most δ 13 C f predictions have poor performance against observations during growing season (Fig. 2).
When δ 13 C f observation was additionally included into model, the posterior GPP and NEE simulations improved slightly at most sites (Appendix S1: Fig. S1). In addition, the uncertainties of carbon flux predictions were significantly reduced (Appendix S1: Fig. S2, we only showed GPP statistical test here, NEE results were similar). For example, using the parameters from E0 experiment, the daily averaged GPP and NEE values at DE-Tha site were 7.32 AE 0.9 g C/d and −3.32 AE 0.6 g C/d. In comparison, the uncertainties of daily fluxes were reduced to 0.63 and 0.42 g C/d using the conditioned parameters from E1 experiment.
The poorly constrained V max25 and g1 (larger uncertainties in E0) without the δ 13 C f observations raises the issue of parameter equifinality. Zhuang (2008, 2009) demonstrated that the terrestrial ecosystem model (TEM) calibrated by using eddy flux tower data was able to reproduce the observed carbon fluxes with very different sets of parameters. Here we also showed that Fig. 2. Comparison between observed (blue) and simulated (red) daily GPP, NEE and δ 13 C f during growing season (May-August) at six forest sites using data-conditioned parameters from E0 experiment. E0 used GPP and NEE observations as constraint.
v www.esajournals.org both model-data fusion experiments had similar results in the posterior carbon fluxes (GPP and NEE) with different V max25 and g1 values, and however, E0 experiment generated greater uncertainties and large discrepancy existed between the δ 13 C f predictions and observations. This indicated that GPP and NEE data may only have weak constraint on relevant parameters in the coupled stomatal-photosynthesis modules. The leaf δ 13 C signal reflects the coupling process between carbon and water fluxes during the photosynthetic activity, enabling these fluxes to be integrated and mutually constrained within the model. These isotopic values are related to both stomatal conductance and canopy photosynthesis, which together determine the CO 2 concentration in the chloroplast (Cc), and consequently the photosynthetic discrimination (Farquhar et al. 1989). Therefore, leaf δ 13 C can shed light on physiological responses to environmental conditions such as drought, temperature, and relative humidity (Williams et al. 2010, Voelker et al. 2014, Hartl-Meier et al. 2015) and provides useful constraint on the photosynthetic and stomatal conductance-related parameters (V max25 and g1 in our study). Our model-data fusion results at multiple forest sites were consistent with other studies (Aranibar et al. 2006;Duarte et al. 2017), Fig. 3. Comparison between observed (blue) and simulated (red) daily GPP, NEE and δ 13 C f during growing season (May-August) at six forest sites using data-conditioned parameters from E1 experiment. E1 used GPP, NEE, and δ 13 C f observations as constraint.
v www.esajournals.org which also confirmed that relevant photosynthetic parameters can be further constrained by using isotopic signal of leaves. We further found that the modeled canopy conductance in E1 parameterization schemes showed better performance against observation than that in E0 experiment (Appendix S1: Table S1), indicating that the calibrated iTEM in the E1 experiment may better capture the observed response of canopy conductance to environmental change (e.g., vapor pressure deficit).
The slope of stomatal model (g1) By assimilating the δ 13 C f observation (E1 experiment), we found g1 to be spatially variable among the six sites and there was a significantly positive relationship between the g1 estimates and annual precipitation (Fig. 4), indicating that the stomatal conductance of temperate forest is more sensitive under favorable environment. Previous observations reported a large variation in g1 values among plant function types (PFTs) and even within similar or same species (Baldocchi and Meyers 1998; Lin et al. 2015;Miner et al. 2017) and this variation was generally found to increase with growth temperature and soil water availability (Bauerle et al. 2003;Lin et al. 2015). As for the water availability effect on g1, vegetation in less humid region generally has specific structural issues that can endure unfavorable periods (Hacke and Sperry 2001). Such tissue requires large investment of carbon (Wright et al. 2003), resulting in more conservative stomatal behaviors and more efficient water use (Cowan and Farquhar 1977). The reduction of parameterized g1 in Ball-Berry model may help to reduce water loss in hot and dry environment in waterlimited regions, suggesting a limitation of carbon assimilation by water availability. Based on the long-standing theory of optimal photosynthesisstomatal process, Medlyn et al. (2011) proposed a tractable stomatal model, which is functionally equivalent to the widely used empirical stomatal model (e.g., Ball-Berry model). Based on Medlyn's model, Lin et al. (2015) compiled a global database of stomatal conductance from a wide range of field studies and they found that marginal water costs of carbon uptake (similar term Fig. 4. Statistical relationship between annual precipitation and the slope of stomatal model derived from the E1 experiment among the six forest sites. v www.esajournals.org with g1 in Ball-Berry model) was lower in waterlimited regions, which is consistent with our model-data fusion results, indicating the wetter site had larger g1 estimation. No significant statistical relationship is found between g1 and other climate variables (Fig. S3), while Lin et al. (2015) demonstrated that g1 increased with temperature. It is possible that our model-data fusion experiments were based on only six sites, which may not capture the relationship across a wide range of climate gradients. In addition, the global database of stomatal conductance was derived from leaf scale gas measurement, while the iTEM model uses a big-leaf model, indicating the parameterized g1 actually represents the averaged properties of the whole canopy.
We further showed that there was a spatial variability in this constraint on model parameters. For example, using the final dataconditioned g1 parameter (Fig. 3, E1 experiment) at the DETha site for simulations at other five sites resulted in large deviations from the observed annual mean δ 13 C f (Appendix S1: Fig.  S4). This suggested that the single site-level derived parameter may not be valid for other sites and the spatially variable g1 should be considered for accurately simulating δ 13 C f . Although g1 can greatly impact on transpiration and carbon fluxes (Leuning et al. 1998;Luo et al. 2001;Barnard and Bauerle 2013;Bauerle et al. 2014), current earth system models assume it to be spatially constant (e.g., in the Community Land Model 4.5, g1 is 9.0 in deciduous forest) without considering its possible change resulting from the plant adaption to climate gradient. De Kauwe et al. (2015) first attempted to represent the spatially variable g1 parameter in a land surface model based on plant functional types (PFTs) and environmental conditions. They found that the g1 variability greatly affected annual fluxes of transpiration across evergreen needleleaf forest, tundra and C4 grass regions. Here, using the g1 parameter estimates conditioned on both carbon fluxes and isotopic signals (E1 experiment) at the six sites, we incorporated the statistical relationship between g1 and annual precipitation into iTEM and conducted the regional simulation for temperate forest in the Northern Hemisphere (here we refer as E1-R1). We compared the simulation with results using the parameters derived from E0 experiment (E0-R0).
Since no significant relationships were found between the other four parameter estimates and environmental conditions, we simply used the averaged estimates in each PFT (DBF and ENF). Thus, in both regional simulations, only g1 was precipitation dependent in E1-R1, while other parameters were spatially constant within one PFT. It should be noted that a number of studies found that g1 can vary seasonally, in response to soil water availability (Tenhunen et al. 1990;Misson et al. 2004), light (Bunce 1998; Aasamaa and Sõber 2011), and elevated CO 2 (Kellomaki and Wang 1997;Tausz-Posch et al. 2013). The temporal variations have been reported in trees (Falge et al. 1996;Misson et al. 2004;Barnard and Bauerle 2013;Gimeno et al. 2016) and crops (Leakey et al. 2006;Miner et al. 2017). However, other observations show contradictory results (Xu and Baldocchi 2003;Weber and Gates 1990;Liozon et al. 2000;Medlyn et al. 2001) and the extent to which g1 change is still debated. Therefore, the impacts of potential g1 seasonal changes are not considered here.

Regional simulations
Both simulations (E0-R0 and E1-R1) can well capture the spatial GPP and NEP variability across the temperate forest region in the North Hemisphere (Figs. 5, 6). The Gulf Coast and parts of the Southeast US were the most productive regions, with GPP around 2000 g CÁm −2 Áyr −1 and NEP greater than 500 g CÁm −2 Áyr −1 , while the North China and Tibet-Plateau generally had relatively low GPP (500 g CÁm −2 Áyr −1 ) and NEP (100 g CÁm −2 Áyr −1 ). The comparison between the two simulations showed that E0-R0 was predicted to be slightly higher than E1-R1, with the GPP difference ranging from 50 to 150 g CÁm −2 Áyr −1 and NEP difference from 50 to 100 g CÁm −2 Áyr −1 , respectively. Larger differences were observed in the regions characterized by high productivity, such as the Southeast US and East China. We further evaluated the regional simulations with a data-driven product-Fluxcom (http://www.fluxcom.org). The Fluxcom dataset includes energy and carbon flux predictions which were upscaled to the globe based on machine learning methods that integrate FLUXNET site-level observations, satellite remote sensing, and meteorological data (Jung et al. 2019). From 2003 to 2010, the spatiotemporally v www.esajournals.org averaged GPP and NEP in Fluxcom in the study region were 1047 and 534 g CÁm −2 Áyr −1 , which was comparable with the results in E0-R0 (GPP: 1240 g CÁm −2 Áyr −1 ; NEP: 495 g CÁm −2 Áyr −1 ) and E1-R1 (GPP: 1123 g CÁm −2 Áyr −1; NEP: 465 g CÁm −2 Áyr −1 ). However, the parameter estimates from the E0 experiment performed worse than those from E1 experiment in reproducing the annual GPP and NEP from Fluxcom (Appendix S1: Fig. S5). This indicated that the conditioned model parameters with additional δ 13 C f measurement better constrain regional carbon fluxes.
We further found that the uncertainties of estimated GPP and NEP were larger in the simulations without the δ 13 C f constraint (E0-R0), and the spatial differences were obvious in high productivity areas (Figs. 5,6). For example, the standard deviation of GPP was around 150 g Fig. 5. GPP estimation and its uncertainty (stand deviation-std) by E0-R0 simulation in temperate forest regions of the Northern Hemisphere. The differences (_Diff) of GPP and its uncertainty are calculated as E0-R0 minus E1-R1. E0-R0 and E1-R1 stand for the regional simulations using parameters from the E0 and E1 modeldata fusion experiment, respectively (Units: g CÁm −2 Áyr −1 ). v www.esajournals.org CÁm −2 Áyr −1 in the Southeast US in E0 and in E1, around 120 g CÁm −2 Áyr −1 . The uncertainty induced by the wide parameters space also exhibited a strong seasonal variation (Appendix S1: Fig. S6). For both E0-R0 and E1-R1, the GPP and NEP uncertainty was larger in summer season (from June to August), but small in winter (from November to January). This was consistent with the site-level calibration and validation ( Fig. 3), which showed a relatively larger uncertainty in growing season (non-growing season is not presented here). It was also noted that the GPP and NEP uncertainties had the most reduction in summer season when iTEM used additional δ 13 C f constraints. For example, the summer GPP uncertainty in the Southeast US was reduced by 25%, from 40 to 30 g CÁm −2 Ámonth −1 (Appendix S1: Fig. S3). Fig. 6. NEP estimation and its uncertainty (stand deviation-std) by E0-R0 simulation in temperate forest regions of the Northern Hemisphere. The differences (_Diff) of NEP and its uncertainty are calculated as E0-R0 minus E1-R1. E0-R0 and E1-R1 stand for the regional simulations using parameters from the E0 and E1 modeldata fusion experiment, respectively (Units: g CÁm −2 Áyr −1 ). v www.esajournals.org

Implication and limitation
By implementing a rigorous model-data fusion approach at multiple forest sites, our results confirmed that in situ δ 13 C f measurement can consistently constrain LSM parameters and carbon flux uncertainties, implying that δ 13 C f observations may help improve model parameterization at other sites with different climate or species. Furthermore, the significant relationship between g1 and climate gradient was found in the experiment that includes δ 13 C f constraints.
Our model-data fusion also demonstrated that the spatially variable g1 was essential for accurately predicting site-level δ 13 C f , challenging the current modeling paradigm that vegetation parameters share the same values within one PFT. Recent studies began to explore statistical methods (Verheijen et al. 2015;Butler et al. 2017) or probabilistic approach (Wang et al. 2009;Pappas et al. 2016) to explicitly incorporate the spatial variation into ESMs based on the global trait database (TRY) (Kattge et al. 2011). Although these studies present important advances in representing spatially variable parameters, they failed to consider the uncertainties of the sparse database and the derived correlation between vegetation parameters and climate is generally weak for spatial extrapolations (Wright et al. 2005;Ordoñez et al. 2009;Coyle et al. 2014). By comparison, our data fusion approach can effectively constrain model parameters by maximizing the use of in situ observation and remote sensing data with different sources and spatiotemporal resolutions. This allows us to efficiently build the statistical relationship between the data-conditioned parameters and biotic and abiotic factors. Future efforts should also focus on identifying the potential temporal variation in plant parameters (Williams et al. 2009; Liu and Ng 2019) based on more plant trait data.
A few limitations need to be addressed in the future. First, simply assuming that the mesophyll conductance is proportional to the Rubisco content (Evans and von Caemmerer 1996) may bias our predicted carbon fluxes. Leaf gas-exchange measurements suggest that the mesophyll conductance varies with leaf development (Warren 2008), nutrient availability (Warren 2004), radiation (Grassi andMagnani 2005, Niinemets et al. 2006), and leaf temperature (Bernacchi et al. 2002). Although previous studies empirically described the mesophyll conductance as a function of plant physiological traits and climate variability (Suits et al. 2005, Sun et al. 2014, the understanding of the response of mesophyll conductance to environmental conditions is still limited (Bodin et al. 2013). Second, the effect of soil moisture on stomatal conductance is not considered here, which may affect the predictability of plant responses under drought conditions. Soil water stress and vapor pressure deficit can occur simultaneously, and however, the underlying mechanistic or physiological knowledge of the interaction between soil water availability and stoma is often difficult to be modeled (Damour et al. 2010). In addition, plant morphology often changes under water shortage condition (e.g., the increase root to shoot ratio), and these adaptive strategies could in turn attenuate the negative effects on stomatal conductance. Whether stomatal (soil water shortage directly influences the stomatal closure) or non-stomatal (soil water shortage influences the photosynthetic capacity or mesophyll conductance) limitation dominates the carbon flux is still debated (Egea et al. 2010;Keenan et al. 2010, Kauwe et al. 2015, in particular concerning the response at different time scales. Both experiment which links stomatal function to plant traits and environmental variables and theoretical studies would help understand the physiology of observed stomatal behavior, which can be integrated into future modeling studies.

CONCLUSION
Land surface models usually consist of multiple computational modules integrating biological, hydrological, and physical processes and a group of parameters associated with these processes and controls. Carbon isotope discrimination by terrestrial ecosystems, which is affected by variable environmental conditions (e.g., air humidity, temperature, and light), have the potential to further constrain model parameters.
Here we incorporated photosynthetic carbon isotope discrimination data into a process-based ecosystem model (iTEM) using an Ensemble Kalman Filter (EnKF) to constrain model parameters at six temperate forest sites. The carbon flux observations with and without δ 13 C f measurements were used. When only carbon flux observations v www.esajournals.org were assimilated into iTEM, a wide range of parameter values were found to be able to capture in situ observation, resulting in large deviations from δ 13 C f observations and great uncertainties in the posterior carbon fluxes. In contrast, including δ 13 C f measurements can narrow the estimated ranges of maximum carboxylation rate at 25°C (V cmax25 ) and slope of stomatal conductance (g1) and reduce the simulated carbon flux uncertainties. We further found that when δ 13 C f was used as constraint, the g1 parameter estimates co-varied with annual precipitation among the six forest sites, and its spatial variation was essential for accurately predicting site-level δ 13 C f . The comparison of regional simulations using the parameter sets from the two model-data fusion experiments showed that including δ 13 C f measurement led to slightly better performance against the FLUXCOM products and less uncertainties in the carbon flux predictions. Our model-data fusion study at multiple forest sites demonstrated the potential of using δ 13 C f to constrain carbon flux uncertainties.

ACKNOWLEDGMENTS
This study is supported through funded projects to Q. Z. by the NASA Land-Use and Land-Cover Change program (NASANNX09AI26G), Department of Energy (DE-FG02-08ER64599), the NSF Division of Information and Intelligent Systems (NSF-1028291), and the NSF Carbon and Water in the Earth Program (NSF-0630319). The supercomputing resource is provided by the Rosen Center for Advanced Computing at Purdue University.