A global scale mechanistic model of photosynthetic capacity (LUNA V1.0)

. Although plant photosynthetic capacity as determined by the maximum carboxylation rate (i.e., V c , max25 ) and the maximum electron transport rate (i.e., J max25 ) at a reference temperature (generally 25 ◦ C) is known to vary considerably in space and time in response to environmental conditions, it is typically parameterized in Earth system models (ESMs) with tabulated values associated with plant functional types. In this study, we have developed a mechanistic model of leaf utilization of nitrogen for assimilation (LUNA) to predict photosynthetic capacity at the global scale under different environmental conditions. We adopt an optimality hypothesis to nitrogen allocation among light capture, electron transport, carboxylation and respiration. The LUNA model is able to reasonably capture the measured spatial and temporal patterns of photosynthetic capacity as it explains ∼ 55 % of the global variation in observed values of V c , max25 and ∼ 65 % of the variation in the observed values of J max25 . Model simulations with LUNA under current and future climate conditions demonstrate that modeled values of V c , max25 are most affected in high-latitude regions under future cli-mates. ESMs that relate the values of V c , max25 or J max25 to plant functional types only are likely to substantially overes-timate future global photosynthesis.


Introduction
Photosynthesis is one of the major components of the ecosystem carbon cycle (Canadell et al., 2007;Sellers et al., 1997) and is thus a key ingredient of Earth system models (ESMs) (Block and Mauritsen, 2013;Hurrell et al., 2013). Most ESMs are based on the photosynthesis model developed by Farquhar et al. (1980). The maximum carboxylation rate scaled to 25 • C (i.e., V c,max25 ; µmol CO 2 m −2 s −1 ) and the maximum electron transport rate scaled to 25 • C (i.e., J max25 ; µmol electron m −2 s −1 ) in the model have been generally accepted as the main proxies of photosynthetic capacity. V c,max25 and J max25 are key biochemical parameters in photosynthesis models as they control the carbon fixation process (Farquhar et al., 1980). Large spatial and temporal variations in estimates of the gross primary productivity exist across ESMs (Schaefer et al., 2012), which have been partly attributed to uncertainties in V c,max25 (Bonan et al., 2011). Accurate estimates of V c,max25 and J max25 are of paramount importance to simulate the gross primary productivity as errors in these two entities may be exacerbated when upscaling from leaf to ecosystem level (Hanson et al., 2004).
Many studies have demonstrated that it is particularly difficult to predict accurately the global scale variations in V c,max25 and J max25 (Bonan et al., 2011;Rogers, 2014). One important reason that contributes to this rather poor predictability is a lack of understanding of the processes that control the values of V c,max25 and J max25 (Maire et al., 2012;Xu et al., 2012) despite the fact that V c,max25 has been measured and studied more extensively than most other photosynthetic parameters Leuning, 1997;Wullschleger, 1993). Many empirical studies have shown that V c,max25 and J max25 (or field-based surrogates) correlate with leaf nitrogen content (Medlyn et al., 1999;Prentice et al., 2014;Reich et al., 1998;Ryan, 1995;Walker et al., 2014). Therefore, a constant relationship between the leaf nitrogen content and V c,max25 or J max25 is commonly utilized by many ecosystem models (Bonan et al., 2003;Haxeltine and Prentice, 1996;Kattge et al., 2009). However, the relationship between leaf nitrogen content, V c,max25 and J max25 varies with light intensity, temperature, nitrogen availability and the atmospheric CO 2 concentration (Friend, 1991;Reich et al., 1995;Ripullone et al., 2003). Thus, the presumed relationship between V c,max25 , J max25 and leaf nitrogen content might introduce significant simulation biases of future photosynthetic rates, which in turn may also affect predictions of the downstream carbon cycle and other climate processes that are dependent on the modeled photosynthetic rates (Bonan et al., 2011;Knorr and Kattge, 2005;Rogers, 2014).

Published by Copernicus
To better describe the relationship between photosynthetic capacity and its driving environmental conditions, we have developed a global scale mechanistic model of leaf utilization of nitrogen for assimilation (LUNA). This model takes into explicit consideration the key environmental variables including temperature, radiation, humidity, CO 2 concentration and day length to explain the complex dependencies between leaf nitrogen, V c,max25 and J max25 . Using an optimality hypothesis, the LUNA model allocates leaf nitrogen to different processes, thereby predicting the values of V c,max25 and J max25 under different environmental conditions. We estimate the parameters of LUNA by fitting the model against globally distributed observations of V c,max25 and J max25 . We then use the calibrated LUNA model to assess the impacts of future climate change on photosynthesis by estimating the summer-season net photosynthetic rate using LUNA's predicted values of V c,max25 and J max25 under historic and future climate conditions.

Overview
The LUNA model (version 1.0) is based on the nitrogen allocation model developed by Xu et al. (2012), which optimizes nitrogen allocated to light capture, electron transport, carboxylation and respiration. Xu et al. (2012) considered a few model assumptions to derive the optimized nitrogen distributions, including (i) storage nitrogen is allocated to meet requirements to support new tissue production; (ii) respira-tory nitrogen is equal to the demand implied by the sum of maintenance respiration and growth respiration; and (iii) light capture, electron transport and carboxylation are colimiting to maximize photosynthesis. The model of Xu et al. (2012) has been tested for three different sites only without global-scale calibration of its parameters. Here, we expand the work of Xu et al. (2012) by using a global data set of observations of photosynthetic capacity to derive accurate values of the model parameters and by incorporation of several refinements to support global-scale prediction of V c,max25 and J max25 . Specifically, this revised model considers additional environmental variables such as day length and humidity, and honors variations in the balance between the light-limited electron transport rate and the Rubisco-limited carboxylation rate. We use an efficient Markov chain Monte Carlo simulation approach, the Differential Evolution Adaptive Metropolis algorithm (DREAM (ZS) ) (Laloy and Vrugt, 2012;Vrugt et al., 2008Vrugt et al., , 2009, to fit the nitrogen allocation model to a large data set of observed V c,max and J max values collected across a wide range of environmental gradients (Ali et al., 2015). After model fitting, sensitivity analyses are performed to gauge the response of the model to changes in its parameter values and the key environmental drivers including temperature, photosynthetic active radiation, day length, relative humidity and atmospheric CO 2 concentration. Finally, mean summer-season V c,max25 and J max25 values and their impacts on net photosynthesis are estimated for the globe, using climate projections from the Community Climate System Model (CCSM) (Gent et al., 2011).

Model description
The structure of the LUNA model is based on Xu et al. (2012), where the leaf nitrogen is divided into four different pools including structural nitrogen, photosynthetic nitrogen, storage nitrogen and respiratory nitrogen. We assume that plants optimize their nitrogen allocation to maximize the net photosynthetic carbon gain, defined as the gross photosynthesis (A) minus the maintenance respiration for photosynthetic enzymes (R psn ), under given environmental conditions and a leaf nitrogen use strategy as determined by the parameters of the LUNA model. The model includes the following four unitless parameters: (1) J maxb0 , which specifies the baseline proportion of nitrogen allocated for electron transport rate; (2) J maxb1 , which determines response of electron transport rate to light; (3) t c,j 0 , which defines the baseline ratio of Rubisco-limited rate to light-limited rate; and (4) H , which determines the response of electron transport rate to relative humidity. The LUNA model predicts the values of V c,max25 and J max25 based on the optimal amounts of nitrogen allocated for carboxylation and electron transport. The model inputs are area-based leaf nitrogen content, leaf mass per unit leaf area and the driving environmental conditions including temperature, CO 2 , radiation, relative humidity and day length.
It is important to stress here that the outcome of the optimality concept used in LUNA is conditional on the plant's nitrogen use strategies built into the model. Thus, it is possible that the optimal values of V c,max25 and J max25 predicted by the LUNA model for future climate conditions could produce lower values of the net photosynthetic carbon gain than fixed values of V c,max25 and J max25 without the use of a nitrogen use strategy. An example of this is shown in Fig. S1 in the Supplement, where the net photosynthetic carbon gain pertaining to the optimal nitrogen allocations predicted by the LUNA model for the elevated temperature is lower than its counterpart derived from a fixed nitrogen allocation for the ambient temperature. A complete description of the LUNA model and the associated optimality hypothesis and algorithms appears in Appendix A. This optimality approach is introduced and tested by Xu et al. (2012) for only three different sites, and here we evaluate its usefulness and applicability at the global scale with improvements to account for largescale variability. Optimality approaches are important tools for land surface models because they provide testable hypotheses for specific plant functions (Dewar, 2010;Franklin et al., 2012;Schymanski et al., 2009;Thomas and Williams, 2014).

Data and temperature response functions
Details of data collection are reported in Ali et al. (2015). Specifically, we conducted an exhaustive literature search with Google Scholar to obtain publications that contained the words V c,max , J max , maximum carboxylation rate or maximum electron transport rate. To rapidly subset the most appropriate and relevant publications, we use simultaneously the wording leaf nitrogen content, leaf mass per area or specific leaf area. Individual values of V c,max , J max , area-based leaf nitrogen content (LNC a ; g N m −2 leaf) and leaf mass per unit leaf area (LMA; g dry mass m −2 leaf) are then obtained by digitizing the experimental data depicted graphically in the selected papers. We use all of the data from Ali et al. (2015) with the exception of one study that collected seasonal data on V c,max and J max during a prolonged drought (Xu and Baldocchi, 2003), as the LUNA model does not take into consideration the potential enzyme deterioration due to water stress but rather simulates only the optimal nitrogen allocation based on monthly climate conditions. This resulted in a data set of 766 observations of V c,max and 643 data points of J max ranging from the tropics to the arctic with a total of 125 species. The data include evergreen and deciduous species from arctic, boreal, temperate and tropical areas at different times of the year and from various canopy locations (Fig. S2).
To allow comparisons of V c,max and J max data collected at different temperatures, we first standardize our data to a common reference temperature (25 • C). To do this, we employ temperature response functions (TRFs) for V c,max and J max . Because of issues related to the possibility of acclimation to temperature, there is not yet scientific consensus on which TRF to use (Yamori et al., 2006). Therefore, we use two alternative temperature response functions to evaluate the potential impact of our selection of the TRFs on the outcome of our analysis. The first temperature response function (TRF1) uses the formula of Kattge and Knorr (2007), which accounts empirically for the potential of thermal acclimation to growth temperature. Following the Community Land Model version 4.5, the growth temperature is constrained between 11 and 35 • C (Oleson et al., 2013) to limit the extent of acclimation to growth temperatures found in the calibration data set. The second temperature response function (TRF2) does not consider the thermal acclimation by fixing the acclimation coefficients in TRF1 . Please refer to Appendix B for details on TRF1 and TRF2 used herein.
Because the LUNA model is based on the C 3 photosynthetic pathway, in this study, we only consider C 3 species. Typically, plant species are grouped into several simple plant functional types (PFTs) in ESMs because of computational limitations and gaps in ecological knowledge. The LUNA model does not differentiate among the PFTs of C 3 species due to a limited coverage of environmental conditions for individual PFTs. Thus, a single set of parameter estimates is used for the PFTs of all C 3 species.

Parameter estimation
The four parameters of LUNA are difficult to measure directly in the field. In this study, we estimate their values by fitting our model against observations of V c,max25 and J max25 using the DREAM (ZS) algorithm (Vrugt et al., 2008(Vrugt et al., , 2009Laloy and Vrugt, 2012). This method uses differential evolution as a genetic algorithm for population evolution with a Metropolis selection rule to decide whether candidate points should replace their parents or not. This simple Markov chain Monte Carlo (MCMC) method exhibits excellent sampling efficiencies on a wide range of model calibration problems, including multimodal and high-dimensional search problems. A detailed description of DREAM (ZS) appears in Vrugt et al. (2008Vrugt et al. ( , 2009) and Laloy and Vrugt (2012) and interested readers are referred to these publications for further details. Uniform prior parameter distributions are used to constrain the potential parameter values and the Gaussian likelihood function is used to quantify the distance between the modeled V c,max25 and J max25 values and their observed counterparts. Convergence plots of the DREAM (ZS) sampled LUNA parameters to the posterior distribution are presented in Figs. S3 and S4.

Model evaluations
In this study, two different goodness-of-fit metrics are used to quantify the performance of the LUNA model against the V c,max and J max data. These are the coefficient of determination (r 2 ) (Whitley et al., 2011) and the model efficiency (ME) (Whitley et al., 2011). The r 2 metric ranges between 0 and 1 and measures how much of the observed dispersion of V c,max or J max is explained by the model. A related metric, the model efficiency is calculated using where y i andŷ i denote the observed and LUNA-simulated values, respectively, andȳ signifies the mean of the observations. This metric measures the proportion of the variance in V c,max or J max explained by the 1 : 1 line between model predictions and observations (Mayer and Butler, 1993;Medlyn et al., 2005). The ME ranges between 0 and 1, where a ME of unity corresponds to a perfect match between the modeled and measured data and a ME of zero indicates that the model predictions are only as accurate as the mean of the measured data.

Model sensitivity analysis
To better understand how the simulated LUNA output of V c,max25 and J max25 depends on its four parameters and the key model inputs, a one-at-a-time (OAT) sensitivity analysis is performed. In the first analysis, we focus on the model parameters only, and perturb their calibrated values, one at a time, with ±15 %. The second sensitivity analysis considers separately the effect of the key environmental variables on the simulated values of V c,max25 and J max25 and perturbs the mean values of the day length (hours), daytime radiation (W m −2 ), temperature ( • C), relative humidity (unitless) and carbon dioxide (ppm) with ±15 %.

Changes in V c,max25 and J max25 under future climate projections
The global surface temperature could increase as much of 3.9 • C by the year 2100 relative to present day (Friedlingstein et al., 2014), with large variations across different regions of the globe (Raddatz et al., 2007). Given the dependence of photosynthesis on temperature, it is critical to examine how much future photosynthesis is likely to change in different regions. In this study, we investigate how the LUNA predicted values of V c,max25 and J max25 will change under future climate conditions and how they will affect future estimates of A net , the net photosynthesis rate. The impacts of future climate on V c,max25 , J max25 and A net are quantified by calculating their values for the leaf layer at the top canopy during the summer season under historic and future climate conditions. Appendix C summarizes the calculation of A net . We use model outputs from the climate carbon cycle Coupled Model Intercomparison Project Phase 5 (CMIP5) (Meehl et al., 2000) to obtain projections of future climate. Climate modelers have developed four representative concentration pathways (RCPs) for the 21st century (Taylor et al., 2013). Each of them corresponds to a different level of greenhouse gas emission. In this study, we use historic and future climate conditions simulated by the CCSM 4.0 model under scenario RCP8.5, which considers the largest greenhouse gas emissions. We do not consider herein other models and emission scenarios as the main purpose of our study is not to do a complete analysis under all CMIP5 outputs but rather to estimate the potential impact of nitrogen allocation on photosynthesis. Specifically, the 10-year climate record between 1995 and 2004 is used as a benchmark for historic conditions, whereas the climate data between 2090 and 2099 is used for future conditions. We present the LUNA's predicted V c,max25 and J max25 values for the months of the peak growing season. Data from the NOAA Earth System Research Laboratory (Riebeek, 2011) showed that the maximum amount of carbon dioxide drawn from the atmosphere occurs in August and February for the large land masses of the Northern and Southern hemispheres, respectively. As a result, June, July and August are assumed herein to best represent the summer season for the Northern Hemisphere and December, January and February are considered representative for the summer season on the Southern Hemisphere. In this study, V c,max25 and J max25 values are predicted by LUNA using the average values of climate variables for the summer season.
We conduct a third sensitivity analysis to quantify the impacts of future changes in climate variables such as temperature, CO 2 concentration, radiation and relative humidity on the simulated values of V c,max25 and J max25 . While the first two sensitivity analyses in section 2.6 assume current mean climate conditions, this sensitivity analysis investigates global patterns in sensitivity of V c,max25 and J max25 to future changes in climate variables across different biomes. Specifically, we calculate the percentage difference in the LUNA predicted values of V c,max25 and J max25 using historic and future values of each climate variable. All other climate variables are set at their default (or historic) values.

Model-data comparison of V c,max25 and J max25
The DREAM (ZS) algorithm provides us with the posterior means and standard deviations of the four LUNA parameters (Table 1). The calibrated LUNA model explains ∼ 54 % of the variance in the observed values of V c,max25 across all species (Fig. 1a) and ∼ 65 % of the variance in the observed values of J max25 (Fig. 1b) using the temperature response function TRF1. This response function considers explicitly the thermal acclimation. If TRF2 is used in LUNA, the model is able to explain ∼ 57 % of the variance in the observed values of V c,max25 (Fig. 1c) and ∼ 66 % of the variance in the observed values of J max25 (Fig. 1d) across all species. When the LUNA predictions with TRF1 are compared with the observation data with seasonal cycles, the model explains ∼ 67 Table 1. Mean values and standard deviations (parentheses) of LUNA parameters estimated by using the Differential Evolution Adaptive Metropolis Snooker updater (DREAM (ZS) ) sampling technique for temperature response functions TRF1 and TRF2. The parameters include (1) J maxb0 (unitless): baseline proportion of nitrogen allocated for electron transport rate, (2) J maxb1 (unitless): electron transport rate response to light availability, (3) t c,j 0 (unitless): baseline ratio of Rubisco limited rate to light limited rate and (4) H (unitless): electron transport rate response to relative humidity. is explained by the model (Fig. S7b), yet a similar performance is observed for herbaceous plants and trees (Fig. S7a, c). The statistics for the predictions of J max25 are very similar to those reported previously for TRF1 ( Fig. S6d-f).

Model sensitivity analysis
Sensitivity analysis shows that all four LUNA model parameters (Table 1) have a positive effect on V c,max25 (Fig. 2a, c) and J max25 (Fig. 2b, d) regardless which temperature response function is used. The parameter t c,j 0 has the strongest effect on V c,max25 (Fig. 2a, c) while J maxb0 has the strongest impact on J max25 (Fig. 2b, d). Parameter H has a much lesser control on the simulated values of both V c,max25 and J max25 ( Fig. 2a-d).
Sensitivity analysis of the LUNA model output to its main climate variables shows that radiation most strongly affects the simulated V c,max25 values, whereas an increasingly smaller impact is observed for the day length, temperature, CO 2 concentration and relative humidity (Fig. 3a, c). The LUNA predicted values of J max25 appear most sensitive to day length, followed by temperature, radiation, relative humidity and CO 2 concentration ( Fig. 3b, d). These findings are independent of the TRF being used.

Impacts of climate change on V c,max25 and J max25
Across the globe, a similar pattern is observed for TRF1 and TRF2 in the simulated values of V c,max25 and J max25 (Figs. 4 and S8). Under historical climate conditions, the higher latitudes are predicted to have relatively high values of V c,max25 and J max25 while lower latitudes are predicted to have relatively low values of V c,max25 and J max25 (Fig. 4a, c for TRF1; Fig. S8a, c for TRF2). Future climate conditions are likely to decrease significantly the V c,max25 values for most vegetated lands in large part due to a predicted rise in the temperature and CO 2 concentration (Fig. 4b for TRF1 and Fig. S8b for TRF2). A somewhat opposite trend is observed for J max25 with decreasing values at higher latitudes and increasing values at the lower latitudes ( Fig. 4d for TRF1 and Fig. S8b for TRF2).
Our results show that the LUNA-simulated V c,max25 values are most sensitive to the changes in CO 2 concentration, followed by temperature, radiation and relative humid- to changes in model parameters. Each parameter (J maxb0 , J maxb1 , t c,j 0 , and H ) is varied one at a time by ±15 % of its fitted value. The values of environmental variables are held fixed at their mean values with day length = 14 h, daytime radiation = 182 W m −2 , temperature = 14 • C, relative humidity = 0.6 (unitless) and CO 2 concentration = 393 ppm. V c,max25 and J max25 values are first obtained at changed parameter values and the percentage changes in V c,max25 and J max25 are then calculated relative to the baseline values of V c,max25 and J max25 predicted based the default parameter values. Positive values indicate that the increase in a specific model parameter leads to larger values of V c,max25 or J max25 , while negative values indicate that the increase in a specific model parameter leads to smaller values of V c,max25 or J max25 . ity ( Fig. 5a-d for TRF1 and Fig. S9a-d for TRF2). The variable J max25 appears most sensitive to the changes in temperature, and then radiation, relative humidity and CO 2 (Fig. 6a-d for TRF1 and Fig. S10a-d for TRF2). Across the globe, temperature has negative impacts on V c,max25 when using TRF1 (Fig. 5a); however, when TRF2 is used, V c,max25 is found to increase at the lower latitudes (Fig. S9a).
The simulations of LUNA demonstrate that the future summer-season mean photosynthetic rate at the top leaf layer might be substantially overestimated if acclimation of V c,max25 and J max25 under future climate conditions (i.e., using historic values of V c,max25 and J max25 ) is not explicitly considered (Fig. 7a, b). This is especially true for regions with high temperatures (Fig. S11). Future estimates of the summer-season mean photosynthesis rate are much higher under TRF1 than TRF2 (Fig. 7b). The omission of acclimation could lead to a 10.1 and 16.3 % overestimation in the to changes in environmental variables including day length (D), daytime radiation (R), temperature (T ), relative humidity (RH) and CO 2 concentration. Each environmental variable is varied one at a time by ±15 % around their mean values with day length = 14 h, daytime radiation = 182 W m −2 , temperature = 14 • C, relative humidity = 0.6 (unitless) and CO 2 concentration = 393 ppm. The model parameters (J maxb0 , J maxb1 , t c,j 0 , and H ) are held fixed at their fitted values. V c,max25 and J max25 values are first obtained at changed environmental conditions and percentage changes in V c,max25 and J max25 are calculated relative to the baseline values of V c,max25 and J max25 under the mean climate conditions in the data. Positive values indicate that the increase in a specific environmental variable leads to larger values of V c,max25 and J max25, while negative values indicate that the increase in a specific environmental variable leads to smaller values of V c,max25 and J max25 . global net photosynthetic rate at the top canopy layer for TRF1 and TRF2, respectively.

Model limitations
The LUNA model built on the assumption that nitrogen is allocated according to optimality principles explains a large part of the global-scale variability observed in V c,max25 (∼ 55 %) and J max25 (∼ 65 %), regardless which TRF is used. This approach used by LUNA also mimics accurately seasonal cycles and PFT-specific values of V c,max25 and J max25 (Figs. S5-S7), and has a much better overall predictive power than a multiple linear regression with LNC a and LMA as main predictors. Such a linear model is able to explain only, (µmol CO 2 m −2 s −1 ), b: J max25 (µmol electron m −2 s −1 )) and the difference in either V c,max25 (b) or J max25 (d) due to changed climatic conditions in the future. The difference is calculated by subtracting the photosynthetic capacity predicted by the LUNA model under the historical climate conditions from that under the future climate conditions. The historical climate is represented by the 10-year monthly averages over years 1995-2004 and the future climate is represented by the 10-year monthly averages over years 2090-2099. The model is run by using TRF1, which is a temperature response function that considered the thermal acclimation.  for both TRFs, about 22 % of the variance in the observed V c,max25 values (Fig. S12a, d) and approximately 13 % of the variance in the observed J max25 values (Fig. S12b, d). These results suggest that our model includes many of the key variables that determine the spatial and temporal variation of V c,max25 and J max25 across the globe. The remaining portion of the variance that cannot be explained by the LUNA model can be related to variability within the 125 species considered in this study. There are inherent intraspecific variations in leaf traits (Valladares et al., 2000) and in photosynthetic capacity (Moran et al., 2016). Data availability limits the number of species that can be considered in the present analysis and favors a single LUNA calibration for all species. Indeed, the data for individual species normally did not cover a sufficiently large range of environmental conditions. When more data become available for individual species in the future, we can revisit the calibration procedure and fit our model to specific PFTs pending a sufficiently large enough coverage of environmental conditions. We posit that such a model will be able to adequately capture a larger portion of the variability observed in V c,max25 and J max25 .
Other deficiencies of LUNA might be related to unexplored nutrient limitations and other plant physiological properties. For example, low phosphorus concentrations can reduce considerably the nitrogen use efficiency of tropical plants with typically modest to low nitrogen (Cernusak et al., 2010;Reich and Oleksyn, 2004), suggesting that our LUNA model can be enhanced by considering multiple dif-ferent nutrient limitations simultaneously (Goll et al., 2012;Walker et al., 2014;Wang et al., 2010). Our treatment of the photosynthetic capacity can also be improved by incorporating a species-specific mesophyll and stomatal conductance (Medlyn et al., 2011), by analyzing properties such as leaf life span (Wright et al., 2004), and by considering soil pH, nutrient availability and water availability (Maire et al., 2015). Unexplored nutrient limitations and other plant physiological properties could also play a factor in the limitation of our model. For example, the nitrogen use efficiency of tropical plants (typically modest to low nitrogen) can be diminished by low phosphorus (Cernusak et al., 2010;Reich and Oleksyn, 2004), suggesting that our model could be improved by considering multiple nutrient limitations (Goll et al., 2012;Walker et al., 2014;Wang et al., 2010). Our treatment of photosynthetic capacity could also be improved by incorporating species-specific mesophyll and stomatal conductance (Medlyn et al., 2011), by analyzing leaf properties such as leaf life span (Wright et al., 2004), or by considering soil nutrient, soil water availability and soil pH (Maire et al., 2015).
Measurement errors of V c,max25 and J max25 also affect negatively the ability of LUNA to describe perfectly the observational data. These errors can originate from many different sources but are rarely quantified in the literature. They can play a significant role in parameter fitting. Indeed, previous research has shown that the value of C i in the Farquhar et al. (1980) model used to differentiate between Rubisco and Ribulose-1,5-bisphosphate (RuBP) limitations, could be estimated from different methods in the literature (Miao et al., 2009). Furthermore, it is particularly difficult to obtain accurate and biologically realistic estimates of dark respiration (but see Dubois et al., 2007), and consequently, dark respiration is sometimes not reported (Medlyn et al., 2002b).

Importance of environmental control on V c,max25
and J max25 Our model predicts that higher temperatures generally lead to lower values of V c,max25 and J max25 (Fig. 3a, c). As temperature increases, the nitrogen use efficiencies of V c,max and J max also increase and thus plants need a lower amount of nitrogen allocated for carboxylation and electron transport. This is true for all the sites except for V c,max25 in the warmest regions of our planet when TRF2 is used (Fig. S9a). This is explained by a large increase in the nighttime temperature of LUNA (e.g., 22 to 30 • C) as the daytime temperature (e.g., from 31 to 33 • C) is constrained by the maximum temperature for optimization in TRF2 (i.e., 33 • C). To maximize the net photosynthetic carbon gain, the model predicts a higher proportion of nitrogen allocated to carboxylation to compensate for a higher nighttime respiration rate. Therefore, the LUNA model predicts higher values of V c,max25 . Yet, this may result from a deficiency of TRF2 in that this response function does not allow for thermal acclimation under global warming (Lombardozzi et al., 2015). Our model predicts that the future changes in atmospheric CO 2 concentration has a negligible effect on J max25 , a finding that is in agreement with results from other studies (e.g., Maroco et al., 2002). A meta-analysis of 12 Free-Air Carbon dioxide Enrichment (FACE) experiments demonstrated reductions in J max on the order of 5 % but with a 10 % reduction in V c,max25 under elevated CO 2 concentrations (Long et al., 2004). Our model also predicts that the relative humidity has a relatively minor effect on V c,max25 . This may be due to the fact that most of the values of V c,max25 and J max25 in our data set coincide with relatively high values of the humidity. As LUNA does not consider the effects of drought on photosynthesis, it may have underestimated the effects of water scarcity on V c,max25 under low humidity conditions (Xu and Baldocchi, 2003). Under prolonged drought, plants close their stomata and photosynthesis is greatly reduced (Breshears et al., 2008;McDowell, 2011). Without sufficient carbon input, photosynthetic enzymes may degenerate during the high temperatures of a drought, which could decrease V c,max25 substantially (Limousin et al., 2010;Xu and Baldocchi, 2003).
There are many different ways to incorporate environmental controls on V c,max25 and J max25 . One such approach is to use relatively simple empirical relationships between environmental variables and V c,max25 and J max25 (e.g., Ali et al., 2015;Verheijen et al., 2013). Such functions can improve the model performance (Verheijen et al., 2015), yet may exhibit rather poor extrapolative capabilities under future climate conditions. The optimality hypothesis used by LUNA is arguably better rooted in ecologic theory, and is therefore expected to exhibit a better predictive quality when confronted with novel future climate conditions. Indeed, the concept of optimality has been applied to the prediction of many different plant functions and structures under a wide array of  (Cowan and Farquhar, 1977). For photosynthetic capacity optimization, Haxeltine and Prentice (1996) have predicted V c,max25 based on a trade-off analysis of photosynthesis and respiration. This concept has been incorporated in different land surface models such as LPJ-GUESS (Smith et al., 2001) and LPJmL (Sitch et al., 2003). Both LUNA and the model of Haxeltine and Prentice (1996), hereafter conveniently referred to as HP, consider V c,max25 and respiration. The LUNA model is currently limited to prediction at the leaf level only while the HP model is applicable for both the leaf and canopy level. Nevertheless, key improvements of LUNA over HP include an explicit consideration of light capture, electron transport and storage. Furthermore, the parameters of LUNA have been derived from a much larger global data set, with many different environmental conditions.

Importance of changes in V c,max25 and J max25 to future photosynthesis estimation
Our model suggests that most regions of the world will likely experience reductions in V c,max25 (Figs. 4b and S8b) due to global warming. An increase of the temperature (Fig. S13) and atmospheric CO 2 concentration is expected to increase the nitrogen use efficiency of Rubisco and thus plants are able to reduce the amount of nitrogen allocated for Rubisco to reduce the carbon cost required for enzyme maintenance. Similarly, J max25 will also decrease globally, except in regions where the present temperatures of the growing season are relatively high (Fig. S12b). The increase of J max25 can be attributed to leaf temperature limitations and increased shortwave radiation (Figs. S14 and S15). Temperature will have a relatively small impact on nitrogen allocation in regions with historically high temperatures during the growing season because leaf temperature is already close to or higher than the upper limit of optimal nitrogen allocation (42 • C for TRF1 and 33 • C for TRF2). Based on Eq. (A11), higher levels of shortwave solar radiation will increase nitrogen allocation to electron transport (Evans and Poorter, 2001).
If we do not account for the potential acclimation of V c,max25 and J max25 under future climate conditions, our analysis based on the LUNA model indicates that ESM predictions of future global photosynthesis at the uppermost leaf layer will likely be overestimated by as much as 10-16 % if V c,max25 and J max25 are held fixed (Fig. 7). This overestimation is larger for TRF2 (16.3 %) than TRF1 (10.1 %) and can result from the fact that TFR2 does not account for thermal acclimation under future climate conditions. Consequently, LUNA predicts a large nitrogen allocation acclimation under climate change. In both cases, our results suggest that, to reliably predict global plant responses to future climate change, ESMs should take into explicit consideration environmental controls on V c,max25 and J max25 . It has been suggested recently that nitrogen-related factors are not well represented in ESMs (Houlton et al., 2015;Wieder et al., 2015). Our nitrogen partitioning scheme would help remove the prediction bias of future photosynthetic rates, which will also improve considerably related climate processes that are dependent on these predictions (Bonan et al., 2011;Knorr and Kattge, 2005;Rogers, 2014).

Code availability
The LUNA model has been implemented in CLM5.0 and will be made publicly available with its release in 2016. Standalone codes of LUNA are available in MATLAB, FORTRAN and C. These source codes can be obtained from the corresponding author upon request.

A. A. Ali et al.: A global scale mechanistic model of photosynthetic capacity 597 Appendix A: Leaf utilization of nitrogen for assimilation (LUNA) model
The LUNA model considers nitrogen allocation within a given leaf layer in the canopy that has a predefined leafarea-based plant leaf nitrogen content availability (LNC a ; g N m −2 leaf) to support its growth and maintenance. The structure of the LUNA model is adapted from Xu et al. (2012), where the plant nitrogen at the leaf level is divided into four pools: structural nitrogen (N str ; g N m −2 leaf), photosynthetic nitrogen (N psn ; g N m −2 leaf), storage nitrogen (N store ; g N m −2 leaf) and respiratory nitrogen (N resp ; g N m −2 leaf). Namely, The photosynthetic nitrogen, N psn , is further divided into nitrogen for light capture (N lc ; g N m −2 leaf), nitrogen for electron transport (N et ; g N m −2 leaf) and nitrogen for carboxylation (N cb ; g N m −2 leaf). Namely, The structural nitrogen, N str , is calculated as the multiplication of leaf mass per unit area (LMA; g biomass m −2 leaf), and the structural nitrogen content (SNC; g N g −1 biomass). Namely, where SNC is set to be fixed at 0.002 (g N g −1 biomass), based on data on C : N ratio from dead wood (White et al., 2000). The functional leaf nitrogen content (FNC a ; g N m −2 leaf) is defined by subtracting structural nitrogen content, N str , from the total leaf nitrogen content (LNC a ; g N m −2 leaf), We assume that plants optimize their nitrogen allocations (i.e., N store , N resp , N lc , N et , N cb ) to maximize the net photosynthetic carbon gain, defined as the gross photosynthesis (A) minus the maintenance respiration for photosynthetic enzymes (R psn ), under specific environmental conditions and given plant's strategy of leaf nitrogen use. Namely, the solutions of nitrogen allocations {N store , N resp , N lc , N et , N cb } can be estimated as follows, The gross photosynthesis, A, is calculated with a coupled leaf gas exchange model based on the Farquhar et al. (1980) model of photosynthesis and Ball-Berry-type stomatal conductance model (Ball et al., 1987) (See Appendix C for details). The maintenance respiration for photosynthetic enzymes, R psn , is calculated by the multiplication of total photosynthetic nitrogen (N psn ) and the maintenance respiration cost for photosynthetic enzyme (NUE rp , see Appendix D), namely, In the LUNA model, the maximum electron transport rate (J max ; µmol electron m −2 s −1 ) is simulated to have a baseline allocation of nitrogen and additional nitrogen allocation to change depending on the average daytime photosynthetic active radiation (PAR; µmol electron m −2 s −1 ), day length (hours) and air humidity. Specifically, we have The baseline electron transport rate, J max0 (µmol electron m −2 s −1 ), is calculated as follows, where J maxb0 (unitless) is the baseline proportion of nitrogen allocated for electron transport rate. NUE J max (µmol electron s −1 g −1 N) is the nitrogen use efficiency of J max (see Eq. D2 for details). J maxb1 (unitless) is a coefficient determining the response of the electron transport rate to amount of absorbed light (i.e. αPAR). f (day length) is a function specifies the impact of day length (hours) on J max in view that longer day length has been demonstrated by previous studies to alter V c,max25 and J max25 (Bauerle et al., 2012;Comstock and Ehleringer, 1986) through photoperiod sensing and regulation (e.g., Song et al., 2013). Following Bauerle et al. (2012), f (day length) is simulated as follows, where f (humidity) represents the impact of air humidity on J max . We assume that higher humidity leads to higher J max with less water limitation on stomata opening and that low relative humidity has a stronger impact on nitrogen allocation due to greater water limitation. When relative humidity (RH; unitless) is too low, we assume that plants are physiologically unable to reallocate nitrogen. We therefore assume that there exists a critical value of relative humidity (RH 0 = 0.25; unitless), below which there is no optimal nitrogen allocation. Based on the above assumptions, we have where H (unitless) specifies the impact of relative humidity on electron transport rate. Replacing Eq. (A7) with Eqs. (A8), (A9) and (A10), we have The efficiency of light energy absorption (unitless), α, is calculated depending on the amount of nitrogen allocated for light capture, N lc . Following Niinemets and Tenhunen (1997) we have, where 0.292 is the conversion factor from photon to electron. C b is the conversion factor (1.78) from nitrogen to chlorophyll. After we estimate J max , the actual electron transport rate with the daily maximum radiation (J x ) can be calculated using the empirical expression of Smith (1937), where PAR max (µmol m −2 s −1 ) is the maximum photosynthetically active radiation during the day. Based on Farquhar et al. (1980) and Wullschleger (1993), we calculate the electron-limited photosynthetic rate under daily maximum radiation (W J x ) and the Rubisco-limited photosynthetic rate (W c ) as follows, where K j and K c are the conversion factors from V c,max to W c and from J x to W J x , respectively (see Eqs. C4 and C6 in Appendix C for details of calculation). Based on Xu et al. (2012), Maire et al. (2012) and Walker et al. (2014), we assume that W c is proportional to W J x . Specifically, we have where t c,j is the ratio of W c to W J x . We recognize that this ratio may change depending on the nitrogen use efficiency of carboxylation and electron transport (Ainsworth and Rogers, 2007) and therefore introduce the modification as follows, where t c,j 0 (unitless) is the ratio of Rubisco limited rate to light limited rate, NUE c0 (µmol CO 2 s −1 g −1 N) and NUE j 0 (µmol CO 2 s −1 g −1 N) are the daily nitrogen use efficiency of W c and W j under reference climate conditions defined as the 25 • C leaf temperature and atmospheric CO 2 concentration of 380 ppm, with leaf internal CO 2 concentration set as 70 % of the atmospheric CO 2 concentration. NUE c (µmol CO 2 s −1 g −1 N) and NUE j (µmol CO 2 s −1 g −1 N) are the nitrogen use efficiency of W c and W j at the current climate conditions. See Eqs. (D6) and (D7) for details of calculation. The term NUE c NUE j NUE c0 NUE j 0 assumes that the higher nitrogen use efficiency of W c compared to that of W j will lead to a higher value of t c,j given the same value of W j . The exponent 0.5 is used to ensure that the response of V c,max to elevated CO 2 is down-regulated by approximately 10 % when CO 2 increased from 365 to 567 ppm as reported by Ainsworth and Rogers (2007). Replacing Eq. (A16) with Eqs. (A14), (A15) and (A17), we are able to estimate the maximum carboxylation rate (V c,max ; µmol CO 2 m −2 s −1 ) as follows, Following Collatz et al. (1991), the total respiration (R t ) is calculated in proportion to V c,max , Accounting for the daytime and nighttime temperature, we are able to estimate the daily respiration as follows, where D day and D night are daytime and nighttime durations in seconds. f r (T night ) and f r (T day ) are the temperature response functions for respiration (see Eq. B1 for details).
In summary, given an initial estimation of N lc , we are able to first estimate the efficiency of light energy absorption α using Eq. (A12). With that, we are able to estimate the maximum electron transport rate , J max , using Eq. (A11). The nitrogen allocated for electron transport can thus be calculated as follows, Then, based on Eq. (A18), we are able to estimate the corresponding the maximum carboxylation rate V c,max and the nitrogen allocated for carboxylation as follows, where NUE V c,max is the nitrogen use efficiency for V c,max . See Eq. (D1) for details of calculation. Using Eq. (A20), we are able to estimate R td and thus the nitrogen allocated for respiration as follows, where NUE r is nitrogen use efficiency of enzymes for respiration. See Eq. (D3) for details of calculation. Finally, the storage nitrogen is calculated as follows, Note that this storage nitrogen is mainly a remaining component of FNC a . Its formulation is different from the formulation of Xu et al. (2012) where N store is set as a linear function of net photosynthetic rate. This modification is based on the observation that the preliminary fitting to data using the linear function shows no dependence of N store on net photosynthetic rate. To make the solutions realistic, we set the minimum of N store as 5 % of FNC a in view of potential nitrogen for plant functionality that is not accounted for by photosynthesis and respiration. By exploring different values of nitrogen allocated for light capture N lc and using the Eqs. (A21)-(A23), we will find the optimal nitrogen allocations (N store,Nresp,Nlc ,N et,Ncb ) until the net photosynthetic rate is maximized (see Eq. A5) given a specific set of nitrogen allocation coefficients (i.e., J maxb0 , J maxb1 , H, and t c,j 0 ). The detailed optimization algorithms are implemented as follows: 1. increase the nitrogen allocated (N lc ) for light capture (from a small initial value of 0.05) and calculate the corresponding light absorption rate α with Eq. (A12); 2. calculate J max from Eq. (A11) and derive the nitrogen allocated to electron transport, N et , using Eq. (A21); 3. calculate V c,max from Eq. (A18) and derive the nitrogen allocated to Rubisco, N cb , using Eq. (A22); 4. calculate the total respiration R td from Eq. (A20) and derive the nitrogen allocated to respiration, N resp , using Eq. (A23); 5. calculate the total nitrogen invest in photosynthetic enzymes including nitrogen for electron transport, carboxylation and light capture using Eq. (A2); 6. calculate the gross photosynthetic rate, A, and the maintenance respiration for photosynthetic enzymes, R psn , by Eq. (A6); 7. repeat steps (1) to (6) until the increase from previous time step in A is smaller than or equal to the increase in R psn .
Since the response of V c,max and J max to increasing temperature shows a steady rise to an optimum followed by a relatively rapid decline (Bernacchi et al., 2003;Kattge and Knorr, 2007;Leuning, 2002;Medlyn et al., 2002a), we postulate that the detrimental heat stress on leaf enzymatic activity beyond this optimum (Crafts-Brandner and Law, 2000;Crafts-Brandner and Salvucci, 2000;Law and Crafts-Brandner, 1999;Spreitzer and Salvucci, 2002) will cause the leaf to fail to optimize its nitrogen allocation. Consequently, we hypothesized that plants only optimize nitrogen allocation up to their optimum enzymatic activity, which is 42 • C for TRF1 and 33 • C for TRF2. Regardless of whether plants acclimate to temperature or not, we assume that they do not optimally allocate nitrogen when leaf temperature is below 5 • C because low temperatures could substantially limit plant enzymes (Martin et al., 1978;Öquist et al., 1980;Strand and Öquist, 1988).
After we get the optimal nitrogen allocations (N store,Nresp,Nlc ,N et,Ncb ), we are able to estimate the V c,max25 and J max25 by rearranging Eqs. (A21) and (A22) as follows: where NUE V c,max25 and NUE J max25 are the nitrogen use efficiency for V c,max25 and J max25 . See Eqs. (D1) and (D2) in Appendix D for details of calculations.

B1 Temperature dependence of Rubisco properties and respiration
The temperature dependence of Rubisco kinetic parameters (K c , K o , τ ) and mitochondrial respiration in light (R d ) (Farquhar et al., 1980) is an Arrhenius function taken from Bernacchi et al. (2001). The temperature response functions of Rubisco kinetic parameters used are outlined below, which are the same irrespective of whether plants are assumed to acclimate to growth temperatures (Temperature response function one; TRF1) or not (Temperature response function two; TRF2). Community land model version 4.5 (CLM4.5) (Oleson et al., 2013) uses the partial pressures of oxygen, O as 20 900 Pa. The kinetic properties of Rubisco, which depend on temperature are Rubisco-specific factor, τ (Jordan and Ogren, 1984), K cc and K o , which are the Michaelis-Menten constants for CO 2 and O 2 , respectively. The temperature response function of R d and kinetic properties of Rubisco (K cc , K o , τ ) are described below, where the fixed coefficients of the equations are values at 25 • C.
In the above equations, R is the universal gas constant (8.314 J mol −1 K −1 ), T 1 is the leaf temperature (K) and the reference temperature, T 0 = 298.15 K.

B2 Temperature dependence of V c,max and J max
Temperature sensitivities of V c,max and J max are simulated using a modified Arrhenius function (e.g., Kattge and Knorr, 2007;Medlyn et al., 2002a;Walker et al., 2014). Because the temperature relationship could acclimate, we examined Kattge and Knorr (2007)'s formulation of with and without temperature acclimation to plant growth temperature. We use two temperature dependence functions of V c,max and J max , which are described below. Fundamentally, TRF1 is a temperature dependence function for V c,max and J max , which is based on the formulation and parameterization as in Medlyn et al. (2002a) but further modified by Kattge and Knorr (2007) to make the temperature optima a function of growth temperature (T g ; • C). It is specified for V c,max as follows: with where V c,max25 is the value of V c,max at the reference temperature (T 0 = 298.15 K). H a (J mol −1 ) is energy of activation and H d (J mol −1 ) is the energy of deactivation. R is the universal gas constant (8.314 J mol −1 K −1 ) and T 1 (K) is the leaf temperature. The entropy term, S v (J mol −1 K −1 ), is now a function of temperature , where a and b are acclimation parameters. TRF1 is implemented in CLM4.5 by Oleson et al. (2013), who use the form of temperature dependence function for V c,max and J max as shown in Eq. (B5), but with limited temperature acclimation, where S v = 668.39 − 1.07min max T g , 11 , 35 with T g representing the monthly mean air temperature ( • C). Other parameters present in CLM4.5 model include, H a = 72 000 J mol −1 and H d = 200 000 J mol −1 . The values of the acclimation parameters (a = 668.39 and b = −1.07) are taken from Table 3 of Kattge and Knorr (2007).
An equation similar to Eq. (B6), f J max T 1 , T g , is used to describe the temperature dependence of J max that considers temperature acclimation based on the S v term. The values of the acclimation parameters (a and b) for S v are taken from Table 3 of Kattge and Knorr (2007). Following Kattge and Knorr (2007) and CLM4.5, we set H a and H d for J max as 50 000 and 200 000 J mol −1 , respectively.

B4 Temperature response function two (TRF2)
TRF2 does not consider the thermal acclimation. The formulation of TRF2 is same as TRF1 except that in TRF2, the entropy term; S v (J mol −1 K −1 ) is fixed across our data set. The values of S v are taken from Table 3 of Kattge and Knorr (2007). S v is set as 649.12 J mol −1 K −1 and 646.22 J mol −1 K −1 for V c,max25 and J max25 , respectively.

C1 Overview
Photosynthesis is described using a system of three equations and three unknown variables. The three unknown variables include (1) the net rate of leaf photosynthesis (A), (2) the stomatal conductance (g s ) and (3) the intercellular partial pressure of CO 2 (C i ). All of the unknown variables influence one another. The three equations include (1) the Farquhar's non-linear equation (A vs. C i ), (2) the Ball-Berry equation (g s vs. A) and (3) the diffusion equation (A = g s (C a − C i )). We solved all of these equations simultaneously by taking an iterative approach (Collatz et al., 1991;Harley et al., 1992;Leuning, 1990). The detailed algorithm for modeling photosynthesis is described below.

C2 Modeling photosynthesis
The photosynthetic rate (A) depends upon (i) the amount, activity and kinetic properties of Rubisco, and (ii) the rate of ribulose-l,5 bisphosphate (RuBP) regeneration via electron transport (Farquhar et al., 1980). The minimum of these two limiting conditions yields the following expression, where W c is the Rubisco limited rate and W j is the electron transport limited rate. The Rubisco-limited carboxylation can be described by with where V c,max is the maximum rate of carboxylation, competitive with respect to both CO 2 and oxygen, and K cc and K o are Michaelis constants for carboxylation and oxygenation, respectively. τ is the specificity factor for Rubisco (Jordan and Ogren, 1984), while C i and O are the partial pressures of CO 2 and O 2 in the intercellular air space, respectively. Likewise, the electron-limited rate of carboxylation can be expressed by with where J is the potential rate of electron transport, and the factor 4 indicates that the transport of four electrons will generate sufficient ATP and NADPH for the regeneration of RuBP in the Calvin cycle (Farquhar and von Caemmerer, 1982). The potential rate of electron transport is dependent upon irradiance, I , according to the empirical expression of Smith (1937): where α, the efficiency of light energy conversion is considered as 0.292 (unitless) (Niinemets and Tenhunen, 1997) and J max is the maximum rate of electron transport.

C3 Ball-Berry model
The stomatal conductance (g, m s −1 ) is evaluated by the Ball-Berry empirical stomatal conductance model (Ball et al., 1987): where RH is the relative humidity (unitless) at the leaf surface, C a is the CO 2 concentration at the leaf surface, and g 0 (0.0005 m s −1 ) and m are the minimum stomatal conductance and slope (9; constant across all C 3 species), respectively.
The estimation of A could be sensitive to the choice of maximum stomatal conductance slope, as this parameter varies both within and across species (Harley and Baldocchi, 1995;Wilson et al., 2001). A recent synthesis provides the first analysis of the global variation in stomatal slope based on an alternative algorithm that considers representation of optimal stomatal behavior (Lin et al., 2015). However, following CLM4.5, which uses the Ball-Berry empirical stomatal conductance model (Ball et al., 1987), we fixed the value of stomatal slope (m) as 9 for all PFTs in our study.

C4 Calculation of photosynthesis and stomatal conductance
We solved Farquhar's non-linear equation (A vs. C i ), the Ball-Berry equation (g s vs. A) and the diffusion equation (A = g s (C a − C i ) simultaneously by taking an iterative approach (Collatz et al., 1991;Harley et al., 1992;Leuning, 1990) until values of A, g s and C i are obtained. The three equations are solved in two phases; the first phase included solving the equations for which Rubisco is limiting while the second phase considered light limitation. The following steps are followed: 1. Given the initial values of C i (where initial value of C i is assumed 0.7 × ambient CO 2 concentration), the temperature dependence functions of V c,max and J max (see Appendix B), and the temperature dependence of Rubisco kinetics (O, τ , K c and K o , Appendix B), A is calculated from Eq. (C2).
2. CO 2 concentration at the leaf surface (C a ) is determined by calculating the difference between C i and the partial pressure due to A, wind speed and the dimension of the leaf.
3. Given A and C a , the stomatal conductance (g s ) is determined using Eq. (C8).
4. C i is determined by calculating the difference between C a and partial pressure due to A and boundary conditions of the stomata.
5. Using the leaf energy balance based on absorbed shortwave radiation, molar latent heat content of water vapor, air temperature, and a parameter that governs the rate of convective cooling (38.4 J m −2 s −1 K −1 ) (Jarvis, 1986;Moorcroft et al., 2001), leaf temperature is calculated.
The above five steps are repeated in a systematic way until C i is equilibrated and the final value of A is then recorded. The net photosynthetic rate, A net , is then calculated by subtracting the respiration from A as follows: where R t is respiration calculated by Eq. (A19).

Appendix D: Nitrogen use efficiencies
The nitrogen use efficiency for V c,max (NUE V c,max , µmol CO 2 g −1 N s −1 ) is estimated from a baseline nitrogen use efficiency 25 • C (NUE V c,max25 ) and a corresponding temperature response function at as follows: with NUE V c,max25 = 47.3 × 6.25, where the constant 47.3 is the specific Rubisco activity (µmol CO 2 g −1 Rubisco s −1 ) measured at 25 • C and the constant 6.25 is the nitrogen binding factor for Rubisco (g Rubisco g −1 N) (Rogers, 2014). f V c,max (T , T g ) is the function specifying the temperature dependence of V c,max with T as the leaf temperature (K) and T g as the growth air temperature (see Appendix B for details of the temperature dependence of V c, max ). The nitrogen use efficiency for J max (NUE J max , µmol electron g −1 N s −1 ) is estimated based on a characteristic protein cytochrome f (Evans and Poorter, 2001), with NUE J max25 = 8.06 × 156, where the coefficient 156 is the maximum electron transport rate for cytochrome f at 25 • C (µmol electron/µmol cytochrome f ); 8.06 is the nitrogen binding coefficient for cytochrome f (µmol cytochrome f g −1 N in bioenergetics). f J max (T , T g ) is a function specifies the dependence of J max on temperature (see Appendix B for details of the temperature dependence of J max ).
The nitrogen use efficiency of enzymes for respiration (µmol CO 2 g −1 N day −1 ), NUE r , is assumed to be temperature-dependent. Specifically, it is calculated as follows, where 33.69 is the specific nitrogen use efficiency for respiration at 25 • C (µmol CO 2 g −1 N s −1 ) (Makino and Osmond, 1991) and f r (T ) specifies the dependence of respiration on temperature. D day and D night is the daytime and nighttime length in seconds. The maintenance respiration cost for all photosynthetic enzymes (NUE rp , µmol CO 2 g −1 N s −1 ) is calculated as follows: where NUE rp25 is the nitrogen use efficiency at 25 • C. NUE rp25 is estimated from the observation of J max25 and V c,max25 as follows: where the total respiration is set as 1.5 % of V c,max (Collatz et al., 1991). We assume that 50 % of the total respiration is used for maintenance respiration (Van Oijen et al., 2010) and 80 % of the maintenance respiration is used for photosynthetic enzyme. In view that the light absorption rate is generally around 80 % (Evans and Poorter, 2001), we set the nitrogen for light capture as 0.2 based on Eq. (A12) in Appendix A. NUE J max25 and NUE V c,max25 are the nitrogen use efficiency for J max25 and V c,max25 estimated from Eqs. (D1) and (D2). In this study, we use the estimated mean value of 0.715 for NUE rp25 based on the data of Ali et al. (2015).
The nitrogen use efficiency for carboxylation (NUE c ) is calculated as the multiplication of conversion factor K c and the nitrogen use efficiency for V c,max follows: where K c is calculated based on the actual internal CO 2 concentrations and leaf temperature (see Eq. C4 for details). Correspondingly, the reference nitrogen use efficiency for carboxylation (NUE c0 ) is calculated using the Eq. (D5) except that K c is calculated based on the reference internal CO 2 concentration of 26.95 Pa and the reference leaf temperature of 25 • C. The reference internal CO 2 concentration is estimated by assuming 70 % of the atmospheric CO 2 concentration of 380 ppm and an air pressure of 101, 325 Pa. The nitrogen use efficiency for electron transport (NUE j ) is calculated as the multiplication of conversion factor K j and the nitrogen use efficiency for J max follows: where K j is calculated based on the actual internal CO 2 concentrations and leaf temperature (see Eq. C6 in Appendix C for details). Correspondingly, the reference nitrogen use efficiency for electron transport (NUE j 0 ) is calculated using the Eq. (D6) except that K j is calculated based on the reference internal CO 2 concentration of 26.95 Pa and the reference leaf temperature of 25 • C. The reference internal CO 2 concentration is estimated by assuming 70 % of the atmospheric CO 2 concentration of 380 ppm and an air pressure of 101, 325 Pa.