Effects of mesophyll conductance on vegetation responses to elevated CO2 concentrations in a land surface model

Abstract Mesophyll conductance (g m) is known to affect plant photosynthesis. However, g m is rarely explicitly considered in land surface models (LSMs), with the consequence that its role in ecosystem and large‐scale carbon and water fluxes is poorly understood. In particular, the different magnitudes of g m across plant functional types (PFTs) are expected to cause spatially divergent vegetation responses to elevated CO2 concentrations. Here, an extensive literature compilation of g m across major vegetation types is used to parameterize an empirical model of g m in the LSM JSBACH and to adjust photosynthetic parameters based on simulated A n − C i curves. We demonstrate that an explicit representation of g m changes the response of photosynthesis to environmental factors, which cannot be entirely compensated by adjusting photosynthetic parameters. These altered responses lead to changes in the photosynthetic sensitivity to atmospheric CO2 concentrations which depend both on the magnitude of g m and the climatic conditions, particularly temperature. We then conducted simulations under ambient and elevated (ambient + 200 μmol/mol) CO2 concentrations for contrasting ecosystems and for historical and anticipated future climate conditions (representative concentration pathways; RCPs) globally. The g m‐explicit simulations using the RCP8.5 scenario resulted in significantly higher increases in gross primary productivity (GPP) in high latitudes (+10% to + 25%), intermediate increases in temperate regions (+5% to + 15%), and slightly lower to moderately higher responses in tropical regions (−2% to +5%), which summed up to moderate GPP increases globally. Similar patterns were found for transpiration, but with a lower magnitude. Our results suggest that the effect of an explicit representation of g m is most important for simulated carbon and water fluxes in the boreal zone, where a cold climate coincides with evergreen vegetation.


| INTRODUC TI ON
The representation of photosynthesis in land surface models (LSMs) is critical for simulating the response of the terrestrial biosphere to global environmental change (Booth et al., 2012;Rogers et al., 2017), the land uptake of CO 2 , as well as the coupling of the water and carbon cycles. The photosynthesis schemes embedded within state-of-the-art LSMs commonly assume that the CO 2 concentration available for carboxylation equals the CO 2 concentration in the sub-stomatal cavity, that is the intercellular CO 2 concentration (C i ). This corresponds to the assumption that the conductance to CO 2 transfer within the leaf (mesophyll conductance, g m ) is infinite, and that the CO 2 concentration at the actual place of carboxylation in the chloroplast stroma (chloroplastic CO 2 concentration, C c ) equals C i . However, evidence has clearly shown that g m is finite (Flexas, Ribas-Carbó, Diaz-Espejo, Galmés, & Medrano, 2008;Warren, 2008) and that it causes a clear drawdown of the CO 2 concentration between the sub-stomatal cavity and the chloroplast stroma. The magnitude of this drawdown depends both on g m and the photosynthetic capacity of the leaf, which is reflected in the definition of g m : g m = A n /(C i − C c ), where A n is net assimilation. Replacing C i with C c as the available CO 2 concentration for photosynthesis has been shown to change the response of simulated photosynthesis to environmental drivers (Niinemets, Díaz-Espejo, Flexas, Galmés, & Warren, 2009), which has important implications for large-scale simulations of land carbon uptake (Sun, Gu, Dickinson, Norby et al., 2014). g m is a complex physiological property which integrates several leaf-internal sub-conductances in both the gaseous and liquid phase, including the intercellular airspace, cell walls, plasma membranes, cytoplasm, and the chloroplast envelopes and stroma (Evans, Kaldenhoff, Genty, & Terashima, 2009). g m is known to change dynamically in response to several environmental stimuli at the time scale of minutes (Warren, 2008). At the same time, its absolute magnitude is constrained by leaf anatomical and structural traits (e.g. cell wall thickness, chloroplast surface area attached to the intercellular airspaces (Tomás et al., 2013)), with the consequence that the values of g m differ considerably among vegetation types (Flexas et al., 2008).
Despite its important role in photosynthesis and the distinct differences across plant functional types (PFTs), g m is at present not explicitly considered in the vast majority of LSMs for two main reasons: (1) the current process understanding of g m is severely limited (Rogers et al., 2017) as its response to environmental drivers, foremost light and CO 2 concentration but also temperature, is largely unknown and currently an area of intensive research (von Caemmerer & Evans, 2015;Tazoe, Caemmerer, Badger, & Evans, 2009;Théroux-Rancourt & Gilbert, 2017;Xiong et al., 2015), and (2) the effects of g m are implicitly included in current models since the overestimation of CO 2 available for photosynthesis is compensated for by an underestimated (apparent) photosynthetic capacity. This means that parameters representing photosynthetic capacity, which are currently estimated on a C i -basis, would need to be re-estimated on a C c -basis if g m were to be explicitly considered in models .
It is likely for these two complications that so far only one study (Sun, Gu, Dickinson, Norby et al., 2014) focused on the effects of an explicit representation of g m in a LSM (the Community Land Model 4.5). Sun, Gu, Dickinson, Norby et al. (2014) showed that the overestimation of the available CO 2 concentration for photosynthesis due to the assumption of an infinite g m leads to an underestimation of the photosynthetic sensitivity to rising atmospheric CO 2 concentrations (C a ). As a consequence, replacing the implicit simulation of g m with an explicit one significantly increased the responsiveness of GPP (+16% from 1901 to 2010) to rising C a as long as C a was not saturating.
The stronger response of photosynthesis to rising atmospheric CO 2 concentrations with an explicitly modeled g m as shown in the study by Sun, Gu, Dickinson, Norby et al. (2014) implies that the physiological responses to rising atmospheric CO 2 concentrations will vary among plant groups that have intrinsically different values of g m . Consequently, it might be hypothesized that photosynthesis of plants with a low g m (e.g. evergreen species) are more responsive to rising atmospheric CO 2 concentrations than plants with a higher g m (e.g. herbaceous plants), which potentially gives the former plant group a relative advantage over the latter in a high CO 2 world (Niinemets, Flexas, & Peñuelas, 2011). A stronger response of photosynthesis to C a is likely to also affect stomatal conductance (g s ) given that g s and A n are tightly coupled (Wong, Cowan, & Farquhar, 1979). As a consequence, the consideration of g m is expected to have important implications for both terrestrial carbon and water fluxes, as well as their coupling (e.g. water-use efficiency, . Such plant type-specific physiological responses would thus not only have important implications for the future global distribution of vegetation types, but also for large-scale patterns of biogeochemical cycles and associated physical climate feedbacks (e.g. evaporative cooling).
In this paper, we explore whether g m has implications for simulations of future global carbon and water fluxes, and to what extent the effects are expected to differ among vegetation types and climatic conditions. In the following, we (1) compile a global database of g m measurements, (2) describe the g m model and its incorporation into the LSM JSBACH (Knauer, Werner, & Zaehle, 2015;Reick, Raddatz, Brovkin, & Gayler, 2013), (3) outline the model parameterization and the necessary adjustment of photosynthetic parameters, (4) analyze the effects of an explicit g m on the photosynthetic sensitivity to C a at the leaf-and ecosystem level, and (5) investigate its relevance for future carbon and water fluxes globally.

| ME THODS
To investigate the effects of g m on simulations of water and carbon fluxes, we tested two different approaches in the LSM JSBACH: Implicit g m : Effects of g m are considered implicitly by employing apparent (C i -based) photosynthetic parameters. This represents the current scenario in most LSMs.
Rubisco kinetic parameters were taken from Bernacchi et al. (2001). This model version is denoted as Imp.
Explicit g m : g m is modeled explicitly as described in Section 2.1. Rubisco kinetic parameters were taken from Bernacchi et al. (2002), and were determined on a C c -basis. Four sub-versions (denoted as Exp, ExpC, ExpL, ExpCL) were implemented, which differ with respect to whether g m is affected by C i and/or light ( Table 1). The effect of these two factors is contentious in the literature Théroux-Rancourt & Gilbert, 2017), hence it is relevant to investigate their potential sensitivities to simulations of photosynthesis at the leaf to the global scale.
Note that the two approaches differ only in the consideration of g m (included implicitly or explicitly) and the Rubisco kinetic parameters (Michaelis-Menten constants for CO 2 (K c ) and O 2 (K o ), photorespiratory CO 2 compensation point (Γ*)) as well as their temperature responses (see Appendix S1 for model formulations and Table S1 for parameter values).

| Mesophyll conductance model
The g m model implemented here is a multiplicative formulation, in which a PFT-specific maximum (i.e. unstressed) value of g m at the reference temperature of 25°C (g m,max25 ) is modified by environmental factors: where N is leaf nitrogen content, T l is leaf temperature, is soil moisture content, C i is intercellular CO 2 concentration, Q a is absorbed photosynthetic photon flux density, and f denotes "function of". g m,min is defined as g m,min = f min × g m,max25 , and accounts for the fact that g m does not decrease to zero even under unfavorable conditions such as severe water stress (e.g. . f min was parameterized from data presented in Delfine, Loreto, Pinelli, Tognetti, and Alvino (2005), Galmés, Abadía, Medrano, and Flexas (2007), and  as f min = 0.15. In Equation (1), g m,max25 and f 1 represent leaf structural determinants of g m , whereas f 2 − f 5 describe instantaneous physiological responses.
Note that the last two terms in Equation (1) (f 4 and f 5 ) are only considered in some model versions (Table 1). Acclimation of g m to elevated CO 2 was not considered in the model as measured g m of plants grown under ambient and elevated CO 2 concentrations did not show consistent differences (Kitao et al., 2015;Mizokami, Noguchi, Kojima, Sakakibara, & Terashima, 2018;Singsaas, Ort, & Delucia, 2003).

| Canopy profile
g m generally declines with depth through the canopy, and is usually higher in sun than in shade leaves (Hanba, Kogami, & Terashima, 2002;Piel, Frak, Roux, & Genty, 2002). It has been found that g m varies in a similar manner to photosynthetic capacity (or N) across the canopy profile (Montpied, Granier, & Dreyer, 2009). This decline with canopy depth might be related to the relatively low mesophyll thickness and the lower chloroplast surface area exposed to the intercellular airspaces in shade-adapted leaves (Evans, Caemmerer, Setchell, & Hudson, 1994;Hanba et al., 2002). Here, we implemented the following canopy profile of g m : where k n is the canopy nitrogen extinction coefficient and L is the leaf area index (LAI). k n was assumed to be 0.11 following Zaehle and Friend (2010). Thus, the canopy gradient of g m equals the one of V cmax and J max in the model. Such a behavior was confirmed by several studies (Han, Iio, Naramoto, & Kakubari, 2010;Montpied et al., 2009;Warren, Löw, Matyssek, & Tausz, 2007), but also higher (Zhang & Yin, 2012) and lower (Cano et al., 2011;Niinemets, Cescatti, Rodeghiero, & Tosens, 2006) gradients for g m compared to V cmax have been found, suggesting that k n is site-and probably PFT-specific (Warren et al., 2007).

| Temperature response
The temperature response of g m is the result of different physical and physiological processes in mesophyll cells (e.g. solubility and diffusivity of CO 2 in water), and the response is likely to differ across cell compartments, for example, membranes and cell walls (Evans & von Caemmerer, 2013). The overall response of g m to leaf temperature can be described by a modified Arrhenius function (Johnson, Eyring, & Williams, 1942): where H a is the activation energy (J/mol), H d is the deactivation energy (J/mol), ΔS is the entropy term (J mol −1 K −1 ) (see Table S1 for parameter values), T l is the leaf temperature (K), T ref is the reference temperature (298.15 K), and R is the universal gas constant (8.314 J mol −1 3) was parameterized according to Bernacchi, Portis, Nakano, Caemmerer, and Long (2002) and shows a temperature optimum close to 35.5°C. The use of the parameter values reported in Bernacchi et al. (2002) is consistent with the C c -based Rubisco kinetic parameters used in this study (K o,Cc , K c,Cc , Γ * Cc , Appendix S1), which were derived assuming the same temperature response of g m (Equation (3)). Published temperature responses of g m differ with respect to the behavior at high temperatures, and both hump-shaped (Egea, González-Real, Baille, Nortes, & Diaz-Espejo, 2011), as well as monotonously increasing responses (Scafaro, Caemmerer, Evans, & Atwell, 2011)

| Soil moisture response
The decline of g m with increasing soil water stress has been widely reported (e.g. Misson, Limousin, Rodriguez, & Letts, 2010;Varone et al., 2012), and has been attributed to the role of aquaporins in leaf-internal CO 2 transport (Miyazawa, Yoshimura, Shinzaki, Maeshima, & Miyake, 2008;Perez-Martin et al., 2014). Here, we implemented the following soil moisture dependence of g m : where is soil moisture (m), wilt is the permanent wilting point (m), below which water stress is at its maximum, and crit is the critical soil moisture content (m), which marks the onset of soil water stress. wilt and crit are constant fractions (0.32 and 0.7 for wilt and crit , respectively) of the field capacity, which is calculated depending on the grain size distribution of the soil. Equation (4) is applied to g m , g s , and leaf biochemistry (V cmax and J max ) but with different sensitivities (i.e. different values of the exponent q, i.e. q m , q s , and q b for g m , g s , and leaf biochemistry, respectively). Using the same formulation, Egea, Verhoef, and Vidale (2011) found that imposing the highest sensitivity to g m , then to g s , and finally to V cmax and J max best captured the behavior of photosynthesis under water stressed conditions for a variety of species from different PFTs.
The q parameters were defined accordingly as q m = 0.75, q s = 0.50, and q b = 0.25 for all PFTs.

| Response to intercellular CO 2 concentration
Most studies investigating the response of g m to C i have found a continuous decrease of g m with increasing C i under field conditions (i.e. C i above c. 200 μmol/mol) Hassiotou, Ludwig, Renton, Veneklaas, & Evans, 2009;Xiong et al., 2015), but see Tazoe et al., 2009). However, there is currently no physiological explanation to link the response of g m to C i , and some concerns on the reliability of these measurements have been raised . We have implemented a C i response function which was derived empirically based on leaf-level measurements as shown in Figure S1.
Equation (5) describes an abrupt increase of g m at low C i (until approx. 100 μmol/mol), and an exponential decline thereafter ( Figure   S1).

| Light response
The effect of absorbed radiation on g m and the underlying mechanisms driving this potential response are currently unresolved.
Studies in which g m was measured at different light levels have reported either no clear responses of g m to variations in light (Loucos, Simonin, & Barbour, 2017;Tazoe et al., 2009;Yamori, Evans, & Caemmerer, 2010), or clear increases with light (Cai, Yang, Li, Wang, & Huang, 2017;Douthe, Dreyer, Brendel, & Warren, 2012;Yin et al., 2009) (see Figure S2). The following function was used to simulate the potential effects of light on g m : where Q a is absorbed photosynthetic photon flux density (μmol m −2 s −1 ). Equation (6) corresponds to a light response curve that takes the values of approximately f min and approximately 1 at Q a values of 0 and 1,500 μmol m −2 s −1 , respectively (see Figure S2).
The function corresponds to a steep increase of g m at low Q a , a value of 0.8 (relative to g m,max25 ) at approximately 500 μmol m −2 s −1 and a shallow response at high Q a .

| C4 plants
In C4 plants, g m describes the conductance to CO 2 transfer from the intercellular airspace to the cytosol of the mesophyll cells, where the first binding of CO 2 occurs. Thus, the main difference to C3 plants is that the chloroplast components are not part of the diffusion pathway (von Caemmerer & Furbank, 1999). This means that g m in C4 plants causes a CO 2 concentration drawdown from C i to C m , the CO 2 concentration in the mesophyll cytosol (i.e. g m = A n /(C i − C m )). Recent methodological advances have enabled measurements of g m in C4 plants (Barbour, Evans, Simonin, & Caemmerer, 2016;Ubierna, Gandin, Boyd, & Cousins, 2017).
These measurements indicate that g m is higher in C4 plants than in C3 plants (see Figure 1). The response of g m to environmental factors was assumed to be identical to that in C3 plants (Equations (2)- (6)). This assumption could not be confirmed due to the scarcity of g m measurements in C4 plants, but recent studies indicate that the temperature as well as the C i response are qualitatively similar to that in C3 plants (Kolbe & Cousins, 2018;Ubierna et al., 2017). The relationship among C4 photosynthetic parameters was kept as in von Caemmerer and Furbank (1999) (see Table S1). (4)

| Implementation into the LSM JSBACH
The developed g m model was incorporated into the LSM JSBACH (Knauer et al., 2015;Reick et al., 2013), which is the land component of the MPI Earth system model (Giorgetta et al., 2013).
Vegetation in JSBACH is classified into PFTs, which may co-occur in model grid cells as tiles, and which differ with respect to key physiological and biophysical properties. Fluxes and conductances are scaled to canopy-level with the LAI for each PFT, and the cover fraction-weighted mean of all tiles gives the respective grid cell value. Land-atmosphere water fluxes are calculated with a bulk transfer approach (Schulz, Dümenil, & Polcher, 2001). Canopy radiative transfer is modeled as described in Wang (2003) based on the model of Goudriaan (1977) and considers sun-lit and shaded canopy fractions in nine vertical layers. g s is modeled according to Medlyn et al. (2011) with PFT-specific stomatal slope parameters (g 1 ) taken from Lin et al. (2015) and a constant residual stomatal conductance (g 0 ) of 0.005 mol m −2 s −1 . A n is simulated according to Farquhar, Caemmerer, and Berry (1980) and von Caemmerer and Furbank (1999) for C3 and C4 vegetation, respectively, and photosynthetic capacity is taken from Kattge, Knorr, Raddatz, and Wirth (2009), and if applicable re-calculated based on N a (leaf nitrogen per area) data in Kattge et al. (2011). In the photosynthesis routine of the model, g m is calculated first according to Equation (1), and A n , g s , C i , and C c are subsequently solved iteratively. In the model versions where g m depends on C i (ExpC and ExpCL), g m is first calculated according to Equation (1) with f 4 set to 1 and then iteratively adjusted for C i (Equation (5)) in the same loop where A n , g s , C i , and C c are solved.

| Maximum mesophyll conductance values (g m,max25 )
To parameterize the key parameter in the model, g m,max25 (Equation 1), we compiled an extensive literature review of leaf-level g m -measurements as described in Appendix S2. This dataset (Appendix S3) adds substantial new data to previous databases (e.g. Flexas et al., 2008) and comprises 609 individual g m measurements of 319 species from 295 studies. Measurements were performed using all common methods used to estimate g m (see e.g. Pons et al., 2009) and represent unstressed, fully expanded, and sun-exposed leaves.
If necessary, measurements were converted to units of mol m −2 s −1 and standardized to 25°C using Equation (3). If g m was assumed to be light-dependent (model versions ExpL and ExpCL), measurements were standardized to high light (1,500 μmol m −2 s −1 ) according to Equation (6). If g m was assumed to be C i dependent (model versions ExpC and ExpCL), g m,max25 was adjusted according to the C i measured along with g m (Equation 5). This adjustment accounts for the fact that different vegetation types operate at different C i /C a (depending on the stomatal behavior in the model (g 0 and g 1 parameters) and the vapor pressure deficit (VPD)). The g m values were assigned to PFTs and the mean, median and standard error of the median were calculated (see Figure 1, Table 3).

| Adjustment of C i -based to C c -based photosynthetic parameters
The explicit representation of g m in photosynthesis models requires that photosynthetic parameters represent C c -based F I G U R E 1 Maximum (unstressed) mesophyll conductance values for different plant functional types (PFTs), standardized to 25°C. Horizontal lines within boxes represent medians, red dots represent means, the lower and upper boundaries of the boxes represent the first and third quantile, respectively, and whiskers represent 1.5 times the interquartile range. PFT abbreviations are: DNF = deciduous needleleaf trees, TDF = tropical deciduous trees, ENF = evergreen needle-leaf trees, DSH = deciduous shrubs, EBF = evergreen broadleaf trees/ shrubs, TRF = tropical evergreen trees, DBF = deciduous broadleaf trees, C3G = C3 herbs and grasses, RSH = raingreen shrubs, C3C = C3 crops, C4G = C4 grasses and herbs, C4C = C4 crops. Data presented here were not standardized to a given C i or to high light, and can be found in Appendix S3 Woody Herbaceous (n = 2) (n = 3) (n = 7) (n = 163)(n = 16) (n = 93) (n = 99) (n = 5) (n = 174) (n = 5) (n = 2) (n = 40) rather than C i -based values, as the latter implicitly include the effects of g m (Ethier & Livingston, 2004). This typically requires that existing (i.e. C i -based) parameters are adjusted to C c -based parameters. Previous approaches for this parameter adjustment focused on the simultaneous derivation of g m , V cmax and J max from A n − C i curves using curve fitting techniques (Gu, Pallardy, Tu, Law, & Wullschleger, 2010;). An alternative approach as applied in this study makes use of independent g m estimates which allow the conversion of A n − C i curves to A n − C c curves and the subsequent re-estimation of photosynthetic parameters on a C c basis. This alternative approach consists of three main steps (illustrated in Figure   S3; R code available at https://bitbucket.org/juergenknauer/ mesophyll_conductance): 1. Simulation of a PFT-specific A n − C i curve under unstressed conditions, saturating light, and 25°C using the current (implicit g m ) photosynthesis routine of the model with C i -based Rubisco parameters from Bernacchi et al. (2001). Under these conditions, g m is assumed to equal g m,max25 .

2.
Calculation of C c from Fick's first law: C c = C i -A n /g m and construction of the corresponding A n − C c curve. Depending on whether g m is assumed to be independent of C i or not, g m is either assumed to be constant or a function of C i (Equation (5)).

3.
Simultaneous fitting of V cmax25 and J max25 to the A n − C c curve calculated in Step 2 using the same model as in step 1, but with C cbased Rubisco parameters taken from Bernacchi et al. (2002). The fitting is done with a non-linear regression routine.

Compared to parameter adjustments based on measured
A n − C i curves, this approach has the advantage of being universally applicable across model types and model structures, and to both C3 and C4 photosynthesis models. This flexibility is achieved by an internally consistent parameter adjustment which is ensured by the employment of the exact same photosynthesis model and parameter values (e.g. leaf day respiration, Rubisco kinetic parameters) for both the parameter adjustment and the actual model simulations. In addition, this approach circumvents uncertainties associated with the determination of g m from A n − C i curves (e.g. assignment of limitation states) by taking independent g m measurements. It follows that no raw data (i.e. A n − C i curves) are required, but instead a sufficient number of g m measurements, from which representative estimates of g m can be inferred.

| Site-level simulations
The JSBACH model was run for eight eddy covariance sites within the FLUXNET network. The sites were selected to cover different PFTs and contrasting hydro-climates (Table 2). Meteorological data for all sites was downloaded from the FLUXNET2015 webpage i -to C c -based parameters were performed as described in Section 2.5. V cmax25,Ci values were taken from Kattge et al. (2009), and if applicable re-calculated based on N a (leaf nitrogen per area) data in Kattge et al. (2011). PFT abbreviations are as in Figure 1.

| Global simulations
To investigate the large-scale implications of an explicit representation of g m in JSBACH, we conducted global runs for the Imp, Exp, and  Table 3.  Figure 1 show g m values that

| Unstressed g m values across PFTs
were not standardized to a given C i or to high light. Accounting for a potential response of g m to light or C i led to only minor changes in the magnitude of g m and its pattern across PFTs (Table 3 and Table S2).

| Parameter adjustment
The required parameter adjustment procedure as described in Section 2.5 led to significant changes to the two key photosynthetic parameters in the model, V cmax25 and J max25 ( Table 3). The C c -based parameters (V cmax25,Cc and J max25,Cc ) account for the lower available CO 2 concentration due to the effects of g m , and are thus usually higher than their C i -based counterparts. For all PFTs, V cmax25 was more strongly affected than J max25 , which resulted in a decrease of the J max25 /V cmax25 ratio. The difference between the C i -based and C c -based parameters depends both on the magnitude of g m and the magnitude of V cmax25,Ci and J max25,Ci , and is highest when g m is low and photosynthetic capacity is high (as e.g. in ENF). Thus, effects are strongest when the CO 2 drawdown from the intercellular airspaces to the chloroplasts (C i − C c ) is high. The decrease of the J max25,Cc /V cmax25,Cc ratio led to a shift of the inflection point (the C i where photosynthetic limitation changes from Rubisco-limited to RuBP (ribulose-1,5-bisphosphate)-limited) to lower C i , which is associated with a higher fraction of photosynthesis occurring in the electron transport-limited domain (see e.g. Figure S3).
In the model versions where g m is affected by C i (ExpC and ExpCL), g m was assumed to change throughout the A n − C i curve according to Equation (5), that means it increases sharply at low C i and decreases continuously thereafter. When accounting for this potential response, the re-adjusted photosynthetic parameters, in particular J max25,Cc , were considerably higher compared to the default version (Exp). The higher J max25,Cc compensates for the low g m simulated at higher C i where RuBP-regeneration is limiting ( Figure S1), and thus maintains the same A n at high C i as in the implicit case. ExpC and ExpCL were thus characterized by significantly higher J max25,Cc /V cmax25,Cc ratios compared to the Exp model version (  (Table S2).
In the C4 photosynthesis model described by von Caemmerer and Furbank (1999), the only parameter affected by the parameter adjustment is the maximum PEP-carboxylation rate (V pmax25 ) ( Figure   S4). In case of a g m,max25 of 0.739 mol m −2 s −1 , the median value observed in C4 crops, V pmax25 increased strongly from 60 (C i -based) to approximately 145 μmol m −2 s −1 (C m -based; Table 3).  (Table S3).

| Effects on simulated leaf-level photosynthesis
Importantly, when temperature and light deviate from the reference conditions, the agreement between the implicit and explicit model deteriorates. This is especially relevant when temperature changes, because g m exhibits a strong temperature response (Equation (3) The ExpC model led to similar curves as shown in Figure 2 ( Figure   S6). Assuming that g m responds to light (ExpL) led to much lower simulations of A n under low light, as well as to higher sensitivities to rising CO 2 throughout the whole C i range ( Figure S7). The deviations between the implicit and explicit model versions caused changes in the relative sensitivity of A n to changes in C i compared to the reference conditions (Figure 2c). In general, A n showed a stronger relative response to C i in the explicit compared to the implicit model at lower temperatures, but the opposite behavior at high temperatures. Note that Figure 2c depicts changes in the sensitivity of A n to C i in relative terms (see Figure caption), which was decreased at high temperatures despite similar slopes of A n − C i in Figure 2a. These contrasts were more pronounced at lower light conditions. It has to be noted that the sensitivities and their relation between the implicit and explicit model version depend on the C i range of interest (shaded areas in Figure 2a,b), which explains the fact that the CO 2 effect of g m as shown in Figure 2c becomes negative already at temperatures lower than 25°C at high Q a .

| Site-level simulations
The integrated response of ecosystem-level A n (A n,canopy ) and g s (canopy conductance, G c ) to changes in atmospheric CO 2 concentrations are analyzed in Figure 3 for an exemplary set of ecosystems from the FLUXNET2015 database. In the implicit model version (Imp), A n,canopy increased under eCO 2 at all sites with C3 vegetation. The relative increases depend on temperature as described previously (Kirschbaum, 1994), and were more pronounced in warm (e.g. GF-Guy) than in cold climates (e.g. FI-Hyy). The explicit model version that does not consider a light and C i response (Exp) showed higher sensitivities of A n,canopy to eCO 2 for most sites, but a significantly lower sensitivity for GF-Guy, which can be explained by the lower photosynthetic sensitivity to CO 2 at higher temperatures ( Figure 2).

F I G U R E 2 (a)
A n − C i curves for the implicit (Imp, blue) and explicit (Exp, orange) model versions for three different temperatures (b) and light conditions and (c) the resulting differences in photosynthetic sensitivity to CO 2 between the implicit and the explicit model version for the grey shaded C i regions in (a) and (b). ΔA n in (c) is defined as ΔA n = (A n,eCO2 − A n,aCO2 )/A n,aCO2 × 100, where aCO 2 denotes the intercellular CO 2 concentration (C i ) range between 240 and 320 μmol/mol (=0.6 × 400 -0.8 × 400 μmol/mol) and eCO 2 denotes the C i range between 360 and 480 μmol/mol (=0.6 × 600 -0.8 × 600 μmol/mol). Positive values in (c) indicate that A n in the explicit model is more sensitive to CO 2 than A n in the implicit model, negative values indicate the opposite. Shown are leaf-level simulations using the C3 photosynthesis model described in Farquhar et al. (1980) with the following parameters: V cmax25,Ci = 40 μmol m −2 s −1 ; J max25,Ci = 76 μmol m −2 s −1 ; g m,max25 = 0.1 mol m −2 s −1 ; V cmax25,Cc = 50.8 μmol m −2 s −1 ; J max25,Cc = 76.86 μmol m −2 s −1 ; R l (respiration rate in light) = 0.44 μmol m −2 s −1 ; and Rubisco kinetic parameters as listed in Table S1 1,500 The model configuration that responds to C i , but not to Q a (ExpC) showed similar or slightly lower responses compared to the Exp model version, which is likely due to the fact that the lower simulated g m was largely compensated by the higher C c -based photosynthetic capacity (Table 3). In contrast, model configurations that simulate a response of g m to light (ExpL and ExpCL) showed the highest responsiveness of A n,canopy to CO 2 , which is a consequence of the continuously higher sensitivity under low light in the ExpL and ExpCL versions due to the marked decrease of g m at low light (Equation (6), Figure S7). This effect is amplified at the canopy level, as a considerable fraction of a closed canopy continuously operates at low light conditions.
The positive responses in A n,canopy were accompanied by negative responses in G c , that is stomatal closure. A consistent pattern in Figure 3 is the opposite response of G c compared to A n,canopy in the sense that a stronger response of A n,canopy was associated with a weaker response of G c , with the result that the response of ecosystem-level intrinsic water-use efficiency (iWUE ,canopy = A n,canopy /G c ) to eCO 2 did not vary among the model runs (i.e. it always increased by the same amount). This can be explained as an intrinsic property of the stomatal model employed here (Medlyn et al., 2011), in which C i /C a is assumed to stay constant with rising As shown in Figure 2, the fact that g m responds to temperature leads to a significantly different temperature response of A n,canopy .
It follows that the photosynthetic sensitivity to CO 2 shows a different response to temperature in the explicit compared to the implicit model versions (Figure 4) As demonstrated in Figures 2-4, the effects of g m on the photosynthetic responses to eCO 2 depend not only on the magnitude of g m (and thus PFT) but also on the environmental conditions, foremost temperature and radiation. To investigate the isolated effects of g m without any confounding meteorological factors, we conducted additional ecosystem-level simulations for the sites US-Ha1 and GF-Guy, in which g m,max25 was varied while keeping the climate forcing unchanged. For these simulations, g m,max25 was reduced F I G U R E 3 Relative responses of ecosystem-level net assimilation and canopy conductance to elevated atmospheric CO 2 concentrations for the five main model versions tested in this study (Imp = implicit g m ; Exp = explicit g m ; C = C i response of g m ; L = light response of g m ). CO 2 response is calculated as (X eCO2 -X aCO2 )/X aCO2 × 100, where X denotes either canopy net assimilation or canopy conductance, and aCO 2 and eCO 2 denote ambient and elevated (ambient + 200 μmol/mol) atmospheric CO 2 concentrations, respectively. PFT abbreviations are as in Figure 1. T air represents the mean annual temperature (see also stepwise from 10,000 (i.e. non-limiting) to 0.075 mol m −2 s −1 , and V cmax25 and J max25 were re-adjusted for each change in g m,max25 as described in Section 2.5. The results demonstrate that the effects of g m on simulations of photosynthesis strengthen when its magnitude decreases ( Figure 5). This is a consequence of the increasing mismatch between the implicit and explicit model versions when g m decreases (Table S3), an effect that amplifies when conditions deviate from those that were used for the parameter adjustment At both sites, the proportion of Rubisco-limited A n,canopy decreased when g m,max25 decreased. Again, this is a consequence of the parameter adjustment (see Section 2.5), in which the stronger changes in V cmax25 compared to J max25 lead to a shift of the inflection point to a lower C i , which is associated with a higher fraction of photosynthesis occurring in the RuBP regeneration-limited domain. This effect (i.e. the change in J max25 /V cmax25 ) increases with a decrease in g m under otherwise equal conditions. In general, this shift towards lower proportions of Rubisco-limited photosynthesis on total canopy level photosynthesis counteracts the higher photosynthetic sensitivity to CO 2 caused by an explicit g m , as photosynthesis in the RuBP regeneration limited region shows a lower sensitivity to rising CO 2 concentrations.

| Global simulations
At the global scale, the widespread substantial increases in mean annual GPP from the historical  to the future RCP8.5 (2070-2099) simulation illustrate the commonly observed CO 2fertilization effect (Figure 6a). Exceptions from this upward trend were found in some semi-arid regions, as well as in parts of the Amazon basin, which experience a drying trend in the climate projections by HadGEM2-ES. Transpiration (Figure 6b), showed weaker absolute responses and a more diverse pattern throughout the globe. In contrast to GPP, transpiration tended to be reduced due to stomatal closure, but this reduction may be offset by increasing VPD in some regions of the earth (Kala et al., 2016). The more moderate RCP4.5 future scenario showed similar patterns, but smaller absolute differences ( Figure S8). Figure 6 further reveals that g m had spatially contrasting effects on the photosynthetic sensitivity to CO 2 . The differences in the ∆ values between the Exp and the Imp version largely reflect both the magnitude of g m (and thus vegetation type), and the environmental conditions as described earlier. It follows that the largest changes could be found in high latitudes, in particular in boreal forests, which show a combination of vegetation with a low g m (ENF and DNF) and a cool climate, both of which increased the photosynthetic sensitivity to CO 2 when g m is modeled explicitly. Changes were moderately positive throughout large parts of the temperate (+5% to + 15%) and semi-arid regions of the earth (+0% to + 5%) and slightly F I G U R E 4 Sensitivity of canopy-level net assimilation (A n,canopy ) to elevated CO 2 concentrations in the implicit (Imp) and explicit model versions (Exp, ExpC, ExpL) for the Mediterranean pine forest site FR-LBr. CO 2 response of A n,canopy is defined as (A n,canopy,eCO2 -A n,canopy,aCO2 )/A n,canopy,aCO2 × 100, where aCO 2 and eCO 2 denote ambient and elevated (ambient + 200 μmol/mol) atmospheric CO 2 concentrations, respectively. Data were filtered to represent periods in the growing season (four most productive months), at daytime (Q a > 200 μmol m −2 s −1 ), and in the absence of soil water stress (β > 0.95; Equation (4)). Points are half-hourly simulation results, and lines indicate local polynomial regression fits (loess) to the points The differences among plant types are more clearly demonstrated in Figure 7 (for the RCP8.5 scenario; see Figure S9 for the RCP4.5 scenario). As stated earlier, the differences among the PFTs are not only caused by differences in g m , but also by differences in the prevailing climatic conditions. For example, the lower response of TDF compared to DNF, despite similar g m,max25 values, can be attributed to the higher temperatures the TDF are exposed to ( Figure 4). Nonetheless, the comparison of co-occurring PFTs in the same model grid cells (i.e. PFTs experiencing identical climate forcing), showed significant differences in the simulated photosynthetic sensitivity, indicating that changes therein can primarily be attributed to differences in g m , and not to climate ( Figure S10). and J max25,Cc were taken directly from leaf-level measurements, as these values are often derived assuming different photosynthetic parameters than the model (see Table S4 for a sensitivity analysis).

The widespread increases in
It should be noted that the original C i -based estimates of V cmax and J max might not represent A n − C i curves well due to the assumption of an infinite g m (Ethier & Livingston, 2004). Any potential bias inherent in the C i -based parameters will be propagated to their C c -based values (Table 3). However, the degree of bias in the C ibased parameters as well as the actual implications for the derived C c -based parameters still needs to be investigated.
The adjustment from apparent to true values resulted in changes to the key parameters V cmax25 and J max25 that are qualitatively comparable to the results of previous adjustments , and that compare well with independently adjusted parameters by Bahar, Hayes, Scafaro, Atkin, and Evans (2018) F I G U R E 6 (a-b) Simulated differences between the RCP8.5 future scenario (2070-2099) and the historical (hist) runs  in mean annual gross primary productivity (GPP) and Transpiration (Transp) in the implicit (Imp) model version.
(c-f) Relative differences between the Imp and the Exp (c-d) and ExpCL (e-f) model versions. Δ is defined as Δ = (X RCP8.5 − X hist )/X hist × 100, where X is either GPP or Transpiration. Regions with an average annual GPP < 200 g C m −2 yr −1 were masked out. Red colors in panels c-f denote stronger increases in GPP or Transpiration in the g m -explicit simulations compared to the g m -implicit simulations an effect that was previously asserted by Sun, Gu, Dickinson, Norby et al. (2014). Most relevant in this context is the strong temperature response of g m (Equation (3)), which leads to a significant deviation of simulated photosynthesis under higher and lower temperatures in the explicit model version. It may be noted that these introduced deviations could be avoided by additionally re-adjusting the activation energy of V cmax,Cc . This would, however, not be in accordance with the adjustment of the Rubisco kinetic parameters as performed in Bernacchi et al. (2002), where changes in the temperature response of A n were entirely attributed to K c and K o , thereby assuming an unchanged activation energy of V cmax . This approach is also justified theoretically since V cmax,Cc , the substrate-saturated photosynthesis rate, is by definition unaffected by g m . We thus argue that the observed changes in the photosynthesis response to temperature are not an artifact.
In this study (as in many others), the assumption was made that Rubisco kinetic parameters as determined in tobacco (i.e. following Bernacchi et al., 2002) adequately represent all PFTs. Recent studies have found notable differences in Rubisco kinetic parameters across plant types and species (Galmés, Hermida-Carrera, Laanisto, & Niinemets, 2016;Walker et al., 2013), and differences in Rubisco properties across plant types and climate conditions, as outlined in Galmés et al. (2016), should be included in future LSMs to better represent the temperature response of photosynthesis across the globe. However, so far no studies have determined Rubisco kinetic parameters on a C c -basis across PFTs, which could be used in LSMs where g m is included explicitly. For use in models, it is essential that C c -based Rubisco kinetic parameters, g m , as well as their temperature responses, are measured on the same set of leaves (as in Bernacchi et al., 2002), in order to ensure consistency across photosynthetic parameters.

| Implications for water and carbon fluxes at ecosystem level
The adjustments to the photosynthesis model cause modest changes to the CO 2 sensitivity of A n,canopy and G c . However, the responses depend on the type of g m model that is implemented.
In the Exp version (no light and C i response), the sensitivity of A n,canopy and G c to CO 2 depends both on the magnitude of g m and the climatic conditions, foremost temperature. Strongest effects were found in cold ecosystems with a low g m (FI-Hyy), but this response does not hold across all climate types, and the tropical site GF-Guy showed the reverse response and a decreasing responsiveness to eCO 2 .
The ExpC version (C i response) does not differ markedly from the Exp version described above for any of the ecosystems investigated here. This indicates that the parameter adjustment is capable of completely offsetting the g m response to C i by a concomitant increase in J max25,Cc . This is an important implication for models as our results indicate that a potential response of g m to C i is not expected F I G U R E 7 Photosynthetic sensitivity to future climate conditions for different plant functional types (PFTs), expressed as the differences between the ∆GPP of the explicit (Exp or ExpCL) and implicit (Imp) g m model versions. ∆GPP was calculated as (∆GPP = GPP RCP8.5 -GPP hist )/ GPP hist × 100, where GPP RCP8.5 and GPP hist denote GPP simulated in the RCP8.5 future scenario (2070-2099) and the historical runs , respectively. Shown are results at tile-level (i.e. GPP is simulated only for the respective PFT and not for the entire grid cell) and only for grid cells where the cover fraction of the respective PFT was at least 30%. PFT abbreviations are as in Figure 1 −10  Figure S7). This effect is amplified at the ecosystem level, where a certain fraction of the canopy operates in sub-saturating light conditions regardless of the amount of incident radiation. This potential light response of g m thus significantly increases the CO 2 sensitivity of all C3 ecosystems investigated here and amplifies the strong positive changes in photosynthetic responsiveness to CO 2 in cold climates, or compensates for the negative response in warmer climates.
The explicit representation of g m did not change simulations of iWUE ,canopy (A n,canopy /G c ), which increased at the same rate as in the implicit model version regardless of the g m model employed. This behavior is a consequence of the implemented stomatal conductance model (Medlyn et al., 2011), which is based on the strong coupling between A n and g s , that causes the C i /C a ratio to stay constant regardless of the atmospheric CO 2 concentration. Since most LSMs employ similar g s models as the one used here (see e.g. Sato, To,  Figure S5).

| Global implications
Global simulations under anticipated future climate suggest clear and regionally contrasting effects of g m on GPP and transpiration.
The differences between the g m -implicit and g m -explicit simulations depend on the projected climate, and on the PFT distribution through vegetation-type differences in the magnitude of g m . The fact that plant groups with a low g m showed stronger responses to eCO 2 than those with a high g m under the same climate generally supports an earlier hypothesis that evergreen species are more likely to have a competitive advantage over other plant types in a high CO 2 world (Niinemets et al., 2011). However, our analysis suggests that this hypothesis does not hold in the tropics, where a low g m led to a decrease in the photosynthetic CO 2 responsiveness (Figures 5b and   6). However, the actual relevance of g m in present and future vegetation dynamics must still be investigated using experimental and modeling approaches.
The replacement of the g m -implicit with a g m -explicit model caused significant changes to simulations of GPP, ranging from 2.3 Pg C yr −1 in the Exp model and the RCP4.5 scenario to 6.6 Pg C yr −1 in the ExpCL model and the RCP8.5 scenario ( Figure S11). About two thirds of this increase were caused by regions north of >30°N, where it mostly occurred in regions covered by boreal forests.
Changes of this magnitude are likely large enough to significantly affect the amplitude of atmospheric CO 2 in the high latitudes, hence g m , which is so far neglected in this context (Forkel et al., 2016;Zeng et al., 2014), should be considered as an additional explanatory factor.
Although our results are broadly consistent with those of Sun, Gu, Dickinson, Norby et al. (2014), our estimates of the GPP response to CO 2 are more moderate. Compared to the 16% increase in the cumulative GPP found by Sun, Gu, Dickinson, Norby et al. (2014), our results (calculated from Figure S13 using Equation (2)

| Future model developments and research needs
Our results emphasize that the absolute value of g m is important for the adjustment of photosynthetic parameters and the associated effects on simulations of photosynthesis. The magnitude of g m is relatively robust for well-sampled PFTs, and similar to the results of earlier data compilations (Flexas et al., 2008), but more measurements are needed for tropical species, deciduous needle-leaf species, and C4 plants. The parameterization of these plant types is important for large-scale simulations, but their maximum g m values can at the moment not be confidently parameterized due to a lack of data. In addition, other plant groups such as ferns should be investigated in future LSMs. These plant groups are characterized by a low g m (Carriquí et al., 2015;Tosens et al., 2016), thus their consideration may have important effects on simulated water and carbon fluxes in models that explicitly simulate g m .
It is clearly desirable to bring empirical formulations of g m as used here and in previous studies (Suits et al., 2005;Sun, Gu, Dickinson, Norby et al., 2014) to a more process-based representation. While existing leaf-level models of g m (Tholen & Zhu, 2011;Tomás et al., 2013) are likely too complex to be parameterized at large scales, global models of g m could be readily improved by relating key model parameters (e.g. g m,max25 ) to both anatomical (e.g. cell wall thickness, mesophyll porosity) (Peguero-Pina et al., 2017;Tomás et al., 2013), and biochemical plant traits (e.g. leaf nitrogen content) (von Caemmerer & Evans, 1991;Xue, Ko, Werner, & Tenhunen, 2017;Yamori, Nagai, & Makino, 2011) within a parsimonious model framework that can be applied across plant types.
Currently, one factor hampering future model development is the poor process understanding of g m , which is associated with the fact that measurements of g m are challenging (Pons et al., 2009). It is particularly critical that the role of environmental factors such as C i and light is unresolved. Here, we tested the potential effects of these two drivers on large-scale simulations of carbon and water fluxes.
We found that a potential C i response does not change model predictions, as its effects would be offset by the adjustment of V cmax25 and J max25 from their apparent to true values. A potential light response of g m , however, would be amplified at canopy level and lead to a significantly higher responsiveness of A n to rising atmospheric CO 2 concentrations. Extrapolated to the global scale, such a leaflevel response would significantly increase global carbon uptake and water loss. It is thus highly relevant that potential measurement artifacts are ruled out , and that the recently put forward hypothesis of an apparent light response (Théroux-Rancourt & Gilbert, 2017) is investigated further, as its existence would mean that the light response of g m as observed at the leaf level should not be implemented in models.