Implementation and evaluation of the unified stomatal optimization approach in the Functionally Assembled Terrestrial Ecosystem Simulator (FATES)

Stomata play a central role in regulating the exchange of carbon dioxide and water vapor between ecosystems and the atmosphere. Their function is represented in land surface models (LSMs) by conductance models. The Functionally Assembled Terrestrial Ecosystem Simulator (FATES) is a dynamic vegetation demography model that can simulate both detailed plant demographic and physiological dynamics. To evaluate the effect of stomatal conductance model formulation on forest water and carbon fluxes in FATES, we implemented an optimality-based stomatal conductance model – the Medlyn (MED) model – that simulates the relationship between photosynthesis (A) and stomatal conductance to water vapor (gsw) as an alternative to the FATES default Ball–Woodrow–Berry (BWB) model. To evaluate how the behavior of FATES is affected by stomatal model choice, we conducted a model sensitivity analysis to explore the response of gsw to climate forcing, including atmospheric CO2 concentration, air temperature, radiation, and vapor pressure deficit in the air (VPDa). We found that modeled gsw values varied greatly between the BWB and MED formulations due to the different default stomatal slope parameters (g1). After harmonizing g1 and holding the stomatal intercept parameter (g0) constant for both model formulations, we found that the divergence in modeled gsw was limited to conditions when the VPDa exceeded 1.5 kPa. We then evaluated model simulation results against measurements from a wet evergreen forest in Panama. Results showed that both the MED and BWB model formulations were able to capture the magnitude and diurnal changes of measured gsw and A but underestimated both by about 30 % when the soil was predicted to be very dry. Comparison of modeled soil water content from FATES to a reanalysis product showed that FATES captured soil drying well, but translation of drying soil to modeled physiology reduced the models’ ability to match observations. Our study suggests that the parameterization of stomatal conductance models and current model response to drought are the critical areas for improving model simulation of CO2 and water fluxes in tropical forests.

Abstract. Stomata play a central role in regulating the exchange of carbon dioxide and water vapor between ecosystems and the atmosphere. Their function is represented in land surface models (LSMs) by conductance models. The Functionally Assembled Terrestrial Ecosystem Simulator (FATES) is a dynamic vegetation demography model that can simulate both detailed plant demographic and physiological dynamics. To evaluate the effect of stomatal conductance model formulation on forest water and carbon fluxes in FATES, we implemented an optimality-based stomatal conductance model -the Medlyn (MED) model -that simulates the relationship between photosynthesis (A) and stomatal conductance to water vapor (g sw ) as an alternative to the FATES default Ball-Woodrow-Berry (BWB) model. To evaluate how the behavior of FATES is affected by stomatal model choice, we conducted a model sensitivity analysis to explore the response of g sw to climate forcing, including atmospheric CO 2 concentration, air temperature, radiation, and vapor pressure deficit in the air (VPD a ). We found that modeled g sw values varied greatly between the BWB and MED formulations due to the different default stomatal slope parameters (g 1 ). After harmonizing g 1 and holding the stomatal intercept parameter (g 0 ) constant for both model formulations, we found that the divergence in modeled g sw was limited to conditions when the VPD a exceeded 1.5 kPa. We then evaluated model simulation results against measurements from a wet evergreen forest in Panama. Results showed that both the MED and BWB model formulations were able to capture the magnitude and diurnal changes of measured g sw and A but underestimated both by about 30 % when the soil was predicted to be very dry. Compari-son of modeled soil water content from FATES to a reanalysis product showed that FATES captured soil drying well, but translation of drying soil to modeled physiology reduced the models' ability to match observations. Our study suggests that the parameterization of stomatal conductance models and current model response to drought are the critical areas for improving model simulation of CO 2 and water fluxes in tropical forests.

Introduction
Global climate change has resulted in significant modifications to Earth's ecosystems through changing weather patterns, including an increased frequency and severity of extreme drought and heatwaves, which has resulted in increased risk for terrestrial vegetation (Pachauri et al., 2014;Reichstein et al., 2013;Gatti et al., 2021). The exchange of water vapor and carbon dioxide between plants and the atmosphere is dominated by transport through stomata (Hetherington and Woodward, 2003;Kala et al., 2016). The mechanisms regulating stomatal opening involve complex biochemical and biophysical processes that are currently not represented in land surface models (LSMs) (Lawson et al., 2014;Buckley and Mott, 2013;Blatt, 2000;Davies et al., 2002). However, a range of much simpler, largely empirical formulations that describe the responses of stomata to their environment have been successfully used by LSMs for many years. Most of them require only two parameters, i.e., the intercept parameter (g 0 ), which is the conductance when pho-tosynthesis (A) is zero, and the slope parameter (g 1 ) that describes the relationship between stomatal conductance to water vapor (g sw ) and a regressor term that includes A and environmental drivers (Damour et al., 2010;Berry et al., 2010). The most widely used representation of g sw , and the default formulation used in the Functionally Assembled Terrestrial Ecosystem Simulator (FATES), is the Ball-Woodrow-Berry model (BWB, Ball et al., 1987), wherein g sw is based on an empirical relationship with leaf net photosynthesis (A net ), carbon dioxide concentration at the leaf surface (C s ), and relative humidity at the leaf surface (H s ) as Eq. (1): The optimality-based unified stomatal conductance model (Medlyn model, MED, Medlyn et al., 2011) is based on the assumption that plants will attempt to maximize carbon gain while minimizing water loss (Cowan and Farquhar, 1977). The MED model has been proposed as an alternative representation of g sw in LSMs (De Kauwe et al., 2015;Lawrence et al., 2019). The basic functional form of the MED model is shown in Eq. (2). One important difference between the BWB and MED formulations is that g sw responds to vapor pressure deficit at the leaf surface (VPD s ) instead of H s : Although the functional form of the MED model is similar to the BWB model, g 1 is based on underlying optimization theory and has a strong theoretical link to plant water use efficiency. The MED model is also favored by many plant physiologists given that g sw responds to VPD s rather than H s (Rogers et al., 2017a). Importantly, the g 1 parameter has also been found to vary significantly across a wide range of different plant functional types (PFTs) and climate regions . Better representation of g sw in LSMs requires efforts to improve the fidelity of g 1 parametrization by PFT. The g 1 parameter can be estimated through field measurement campaigns Wu et al., 2020) or model inversion (Bonan et al., 2011;Fer et al., 2018). For example, De Kauwe et al. (2015) derived a PFT-specific g 1 parameterization from Lin et al. (2015) for the CABLE model and found a significant reduction in annual fluxes of transpiration using MED compared with the original model formulation of CABLE (Leuning, 1995). Despite considerable analysis supporting the adoption of the MED model in LSMs (De Kauwe et al., 2015;Lawrence et al., 2019), the formulation has not been widely adopted and is not common in dynamic vegetation models (Fisher et al., 2018). Exploring plant physiological responses to key environmental variables is emerging as a promising way to understand model representation and evaluate model behaviors (Rogers et al., 2017a;Bonan et al., 2011). Stomatal conductance in the MED and BWB models responds to direct environmental drivers including atmospheric CO 2 concentration and VPD s or H s , as well as indirect drivers like radiation and leaf surface temperature (via photosynthesis) (Franks et al., 2017). Evaluating the BWB and MED formulations in response to changing climate in a complete LSM, wherein atmospheric, ecological, and hydrologic processes are highly coupled, is urgently needed to understand model responses within a larger domain.
Improved projection of the response of ecosystems to global climate change requires improved understanding and model representations of plant responses to a hotter, drier, and CO 2 -enriched future (Sullivan et al., 2020). To mimic the drought effects on ecosystems, some models have included a soil water stress factor (often denoted as β), which is used to reduce the "base rate" of stomatal model parameters: either g 0 (e.g., CLM, Lawrence et al., 2019), g 1 (e.g., G'DAY, Comins and McMurtrie, 1993;O-CN, Zaehle and Friend, 2010;CABLE, De Kauwe et al., 2015), or both (e.g., ORCHIDEE, Guimberteau et al., 2018). In some cases, it is also used to lower the maximum carboxylation rate of Rubisco (V cmax ) (e.g., CLM; O-CN; SIBCASA, Schaefer et al., 2008), both V cmax and the maximum rate of electron transport (J max ) (e.g., G'DAY), or directly A (e.g., JULES, Best et al., 2011;Clark et al., 2011). Reduction in A will further reduce g sw . Some models also consider the soil water stress on mesophyll conductance (e.g., SIBCASA; ORCHIDEE). However, the application of a β factor to different physiological parameters has not been evaluated against measurements for g sw models. Therefore, evaluating different g sw schemes and parameterization with data collected under normal and stressed conditions may help reveal areas for model improvement.
In this study, we explored the impact of stomatal behavior under simulated and realistic environmental conditions in the FATES model (Koven et al., 2020), into which we implemented the MED formulation as an alternative approach to the default BWB formulation. The FATES model is a dynamic vegetation demography model that simulates leaf to ecosystem-scale carbon, water, and energy fluxes, as well as cohort-level plant growth, competition, and mortality processes, enabling FATES to predict the distribution, structure, and composition of vegetation (Fisher et al., 2015;Koven et al., 2020). FATES itself is not a stand-alone model, but instead is used in conjunction with a host land model, and is currently coupled with the Community Land Model (CLM, Lawrence et al., 2019) and the Energy Exascale Earth System Model (E3SM) Land Model (ELM, Holm et al., 2020). Using FATES and the MED and BWB representations we addressed the following questions: (1) how do projected leaflevel and canopy-level CO 2 and water vapor fluxes differ between the BWB and MED formulations in response to key meteorological forcing variables? (2) How do the two model outputs of stomatal conductance and photosynthesis compare to leaf-level gas exchange measurements collected through a dry season in a tropical forest? (3) How does the application of a soil water stress factor affect the simulation of water and carbon cycles during dry periods in a tropical forest?

Implementation of the Medlyn model into FATES
In FATES, leaf-level photosynthesis (A) in C 3 plants is based on the model of Farquhar et al. (1980) as modified by Collatz et al. (1991). A is calculated as the minimum of the RuBPcarboxylase-limited (Rubisco-limited) rate and RuBP regeneration rate (i.e., the light-limited rate). The net photosynthesis rate (A net ) is the difference between A and leaf respiration. Gross primary productivity (GPP) is calculated as the weighted average of the photosynthetic rate from sunlit and shaded leaves, which is integrated through the vertical profile, and finally across all the leaf layers by multiplying exposed leaf area for a given cohort. A cohort is a group of plants with similar disturbance history, height, and PFT type. The leaf area of each cohort is calculated from leaf biomass and specific leaf area (SLA). Leaf biomass is controlled by the processes of phenology, allocation, and turnover. SLA is a PFT-specific parameter.
We implemented the MED stomatal conductance model as an alternative to the BWB model for the calculation of g sw in FATES. Leaf-level g sw is central to the water, CO 2 , and energy cycles in forests. It not only controls the water and CO 2 exchange, but also modifies the energy balance and biochemical processes. Similarly, in FATES, the variable g sw is used to model several processes such as the heat and water transfer and photosynthesis. The calculation of this variable is therefore complex and uses both analytical and numerical solutions to couple the equations describing each process. A detailed description of the implementation can be found in online documentation (FATES Development Team, 2020b). It should be noted that parameters g 0 and V cmax used to calculate A net in Eqs. (1) and (2) are multiplied by an empirical soil moisture stress factor (β) by default in the FATES model. The β factor ranges from 1 when the soil is wet to zero when the soil is dry. The β factor depends on the soil water potential of each soil layer, the root distribution of the PFT, and a plant-dependent response to soil water stress as shown in Eq. (3): where w j is a plant wilting factor for layer j and r j is the fraction of roots in layer j . The soil wilting factor is a bounded linear function of soil matric potential defined by two parameters, the soil water potential at (and above) which stomata are fully open and the value at which stomata are fully closed. The soil matric potential is related to the soil water content, soil texture, and organic matter content. The root fraction is determined by PFT-specific root distribution parameters. For more details on the calculation of the plant wilting factor and the fraction of roots, see the CLM version 4.5 (CLM4.5) technical note (Oleson et al., 2013).

The San Lorenzo, Panama, model test-bed site
Our model simulations were made for a single tropical forest located in Bosque Protector San Lorenzo, Panama (9 • 16 51.71 N, 79 • 58 28.27 W, elevation 25 m), which is a part of the Center for Tropical Forest Science (CTFS) -Forest Global Earth Observatory (ForestGEO). The Smithsonian Tropical Research Institute canopy crane provides access to the top of the forest canopy and allows us to compare our simulations with previous measurements of stomatal conductance and net photosynthesis rate Rogers et al., 2017c). The site is characterized as a moist tropical forest, with mean annual temperature of 26 • C, with only small seasonal variation. The mean annual precipitation is 3300 mm, and 90 % of this precipitation falls during the wet season (May-December). More details about the site can be found in Wright et al. (2003). For our study, we conducted all model simulations using the FATES model coupled with the CLM version 5.0.34 (CLM5). For all simulations, we initialized the FATES model using real-world forest inventory data that provided information on tree size distribution for the whole forest area (Condit et al., 2009) and enabled us to better compare model outputs with field measurements by matching the internal cohort structure with that observed in inventory data. Inventory data from the most recent census (1999) were used as the initial state for the simulations. For simplicity, in our FATES simulations we assumed that the site is populated entirely by the broadleaf evergreen tropical (BET) tree plant functional type.

Sensitivity simulations of FATES with synthetic forcing
The MED and BWB stomatal conductance models differ in the representation of atmospheric dryness as well as the g 1 values. To isolate the influence of structural and parametric differences on FATES simulations using the MED and BWB stomatal models, we employed three model ensemble simulations associated with a BET tree PFT.
For the BWB configuration we used the BWB model with a default g 1 value of 8 (unitless) for the BET tree PFT in FATES. In our MED-default setup, the MED model was parameterized with g 1 set to 4.1 kPa 0.5 to match the best estimate from Lin et al. (2015). To constrain the model difference to structural difference we also ran FATES with the MED model with a g 1 value that was harmonized with the BWB model in FATES, which was abbreviated as MED-B. Here, we assumed g 1 for BWB (g 1b ) to be 8, air temperature to be 25 • , and relative humidity in the air (H a ) to be 0.8 following Franks et al. (2017) in Eq. (4) to obtain a BWBequivalent g 1 = 2.39 kPa 0.5 in the MED-B simulation (g 1m , Eq. 4), wherein VPD a is VPD in the air. For all simulations, 4316 Q. Li et al.: Implementation of the Medlyn model in FATES g 0 was fixed at 1000 µmol m −2 s −1 .
The FATES model is driven by half-hourly longwave radiation, shortwave radiation, air temperature, specific humidity, precipitation, surface pressure, wind speed, and atmospheric CO 2 concentration. These variables modify the leaf conductance by changing the environment at the leaf surface (H s , VPD s , and C s in Eqs. 1 and 2). Following model initialization, the model was run with our synthetic climate forcing data in order to reveal model responses to specific climate forcing. For our synthetic climate forcing, each represented the scenarios of a linear increase in VPD a , air temperature (T a ), photosynthetically active radiation (PAR), and atmospheric CO 2 concentration (CO 2 ), respectively, while other climate forcing data were kept as constant. The details for these scenarios are listed in Table 1. In addition, we set the precipitation to 1.3 mm d −1 , surface pressure to 99 626 Pa, wind speed to 4.8 m s −1 , and longwave radiation to 407.4 W m −2 for all these scenarios, which represent annual average conditions at our field site. Given the physical dependence of saturated water vapor on T a (Ficklin and Novick, 2017), it was necessary to adjust the specific humidity together with T a to keep the VPD a fixed at 1 kPa for the MED-default and MED-B simulations. For the T a scenario in the BWB model, H a was kept at 80 % as H a rather than VPD a was used in the BWB model (Franks et al., 2017). We then studied the responses of g sw , net photosynthesis (A net ), gross primary productivity (GPP), and evapotranspiration (ET) to these drivers for the top layer (averaged across sunlit and shaded leaves) of the canopy. We also checked the number of plants per hectare (nplant) to ensure that cohort density did not change during our simulations.

Evaluation of FATES against in situ measurements
We compared the modeled diurnal g sw and A net of upper canopy layers with measured values (Rogers et al., 2017c;Wu et al., 2020). The data were collected at the San Lorenzo field site at monthly intervals across the dry season and the beginning of the wet season during the strong El Niño-Southern Oscillation (ENSO) year of 2016. The g sw and A net of top-of-canopy leaves of eight species, which all belonged to the BET PFT, were measured across 4 months starting in February 2016 and ending in May 2016. Measurements of g sw and A net were made with an LI-6400 (LI-COR Biosciences, Lincoln, NE, USA) for which the conditions of radiation, humidity, CO 2 concentration, and temperature surrounding the leaves were closely matched to the ambient conditions. Using this dataset, the g 1 values of the BWB and MED models were estimated for each species (see Table 2 in Wu et al., 2020). We used those estimations to parametrize the g sw model in FATES, and we used their g sw and A net measurements to compare with FATES simulation results. It should be noted that the g 1 values we used were varied between 4.43 to 8.3 for the BWB model and between 1.14 and 2.85 kPa 0.5 for the MED model; most values were lower than the defaults for evergreen tropical trees in both of the models, as discussed in Wu et al. (2020). Because g 1 was estimated for BWB and MED models based on the same measurements, g 1 was equivalent for the two models and the simulations resemble MED-B and BWB in Sect. 2.3. An ensemble of simulations with varying measured species-specific g 1 values was carried out to evaluate the impact of stomatal slope parameterization on FATES-simulated g sw and A net . In addition, V cmax at 25 • was set to 63 µmol m −2 s −1 based on the A/C i curves measured at the same time during the 2016 campaign (Rogers et al., 2017b). Other parameters such as J max and leaf dark respiration rate (R dark ) at 25 • were directly calculated by FATES based on their relationship with V cmax or leaf nitrogen content. For our simulations, we used the observed half-hourly weather data including precipitation, air temperature, and humidity from the field site meteorological station as the atmospheric forcing data to drive the FATES simulations (Faybishenko et al., 2019). Atmospheric CO 2 concentration was set to a background level of 403.3 ppm based on data from the NOAA Mauna Loa observatory, which is also very close to the CO 2 concentration inside the leaf chamber of the gas exchange equipment.

Drought effects on physiological parameters
In FATES, a soil water stress factor (β) is used to adjust g 0 and V cmax in the original form of the BWB model (Bonan et al., 2011). For the MED approach we implemented, we also applied the β factor in the same manner as the default setting (Sect. 2.1). However, whether the calculation of the β factor can truly reflect soil water conditions is unclear. To the best of our knowledge, the relevance of the β factor has not been rigorously tested for tropical ecosystems or in comparison with measured g sw and A net . We therefore first compared the modeled soil water content and β factor against soil moisture products from the ECMWF Reanalysis data version 5 (ERA5) (Hersbach et al., 2018). Then we explored whether this formulation of the β factor accurately represents observed physiological responses to soil water stress and whether the stress factor should also be applied to g 1 for both the BWB and MED models. To test this, we designed model simulations (Table 2) to assess how the inclusion of the β factor modifies modeled g sw and A net and the comparison with the measurements. In these simulations, g 1 and V cmax were set as the averages across all the species measured for the BET PFT. The notation "on" indicates that the soil water stress effect is turned on, and "off" indicates that the soil water stress effect is turned off in the simulation.

Model responses to climatic drivers
The responses of g sw and A net to climatic forcing as modeled by the BWB (BWB model with default g 1 ) and MEDdefault (MED model with default g 1 ) simulations had similar shapes (Figs. 1 and 2, blue and black lines). The MEDdefault yielded markedly higher g sw than BWB for all climatic drivers considered with an average difference of 75 % (Fig. 1). The MED-default simulation also resulted in higher estimations of A net , but the increase over the BWB simulation was much smaller (around 15 % on average) than the increase in g sw (Fig. 2). When the VPD a increased above 1.5 kPa, the two models showed strong additional divergence. At a VPD a of 2.0 kPa projected g sw and A net from the BWB simulation were 316 % and 86 % lower than projections from the MED-default simulation (Figs. 1c and 2c). The particularly large divergence between the BWB and the MED-default simulations can be explained by a combination of parametric and structural differences. Comparison of the MED-B (for which we used a parameterization equivalent to that of BWB in the MED model) with the BWB limited potential model deviation to structural difference between the two approaches. Both simulations yielded similar responses of g sw and A net to radiation, temperature, and CO 2 (Figs. 1a, b, d and 2a, b, d; blue and red lines), demonstrating that the differences between the BWB and MED-default settings were attributable to the difference in parameterization associated with g 1 . With a harmonized parameterization of g 1 the divergence between the two models above a VPD a of 1.5 kPa was still readily apparent (Figs. 1c and 2c, blue and red lines). The MED-B simulation showed a slight decrease in g sw with high VPD a , while g sw modeled with the BWB simulation decreased more markedly when VPD a was beyond 1.5 kPa. At 2.0 kPa the BWB simulation projected g sw and A net that were 126 % and 53 % lower than the MED-B simulation. For the temperature response of g sw , BWB and MED-B were very similar although BWB had slightly higher g sw values than MED-B (Fig. 1d, blue and red lines).
In contrast to the g sw response, the differences between BWB and MED-default were generally smaller for A net , except when VPD a was above 1.5 kPa (Fig. 2, blue and black lines). The use of measured g 1 for the MED model (MEDdefault) did not markedly change the magnitude of A net compared with MED-B (Fig. 2, red and black lines). When we explored the ecosystem-scale responses (Figs. 3 and 4), we found that the patterns of ET and GPP mirrored the leaf-level responses described above when using our synthetic climatic drivers. The difference between BWB and MED-B was also apparent when VPD a was above 1.5 kPa (Figs. 3c and 4c, blue and red lines).
To rule out the possibility that these differences were related to changes in underlying plant community structure, we looked for any significant changes in cohort density (number of plants per hectare). Our results showed that there was no significant change (Fig. S1); thus, these ecosystem-scale responses were primarily related to changes in underlying leaflevel physiology.

Model evaluation against field measurements
Before comparing the results of the BWB and MED model representations within FATES against field measurements, we first evaluated the consistency of the site meteorological measurements used to drive FATES simulations with those measured with our gas exchange instruments during the campaign. We anticipated that the two conditions would be comparable as the environmental controls in the gas exchange instruments were set to mimic the ambient conditions just before the leaf measurements. We found that for PAR, H a , and CO 2 concentration, the atmospheric and leaf chamber conditions at the time of measurements were in reasonably close agreement, while the in situ measured T a and VPD a were higher than climate data in all months (Fig. 5a-d).  To account for measurement and natural variability of g 1 across different species, we ran a series of FATES simulations driven by meteorological forcing data with different g 1 values. These experiments showed that FATES A net and g sw were sensitive to different g 1 values for both model formulations ( Fig. 5e-l). The MED model ensemble results for A net and g sw with different g 1 values, represented as the envelopes in Fig. 5e-l, generally overlapped with those from the BWB model, with comparable averages. Compared with field measurements, both models captured the diurnal patterns well (Fig. S2) but tended to underestimate A net and g sw , notably in the month of April by about 30 %, at the peak of the dry season ( Fig. 5g and k).

Water stress factor on physiological parameters
Compared with the ERA5 soil moisture products, FATES generally captured the magnitude and trend of the observed average soil water content at the San Lorenzo site (Fig. 6). FATES also simulated the soil water content well for different layers of the soil column (Fig. S3). By April 2016, at the peak of the dry season in a dry year, the simulated soil moisture stress factor (averaged over all the soil layers) reached an annual minimum (0.7), corresponding to the observed soil moisture drying trend (Fig. 6). FATES also underestimated g sw and A net by the largest margin in April when compared to our field measurements ( Fig. 5g and k). To explore this further, we conducted additional experiments focused on evaluating the use of the β factor to modify g 0 , g 1 , and V cmax . For the month of April in 2016, we compared a range of different model simulation experiments wherein the β factor was applied in different combinations to g 0 , g 1 , and Figure 5. (a-d) Diurnal change in climate forcing, (e-h) model-data comparison of net photosynthesis rate (A net ), and (i-l) model-data comparison of stomatal conductance (g sw ) for four field campaign dates. In panels (a)-(d), lines and filled points represent climate forcing data used in FATES and in situ measurements, respectively. Different colors are for different types: black for T a , red for PAR, blue for VPD a , green for atmospheric CO 2 concentration, and gold for H a . In panels (e)-(l) shaded areas represent the range of FATES model ensemble results with different measured g 1 values for different species, while lines represent the averages of these ensemble results. Blue shaded areas and lines are for results from the MED model, and red is for the BWB model. Gray filled circles for the measured data represent averages across species. Black error bars for the measured data represent the 95 % confidence interval (CI) across species. Columns correspond to days of measurements and are presented in chronological order for 17 February, 10 March, 21 April, and 24 May 2016.
V cmax (Table 2, Fig. 7). The results from Exp. 1 and Exp. 4 showed high overlap, indicating that considering the β effect on g 0 does not influence modeled carbon and water fluxes. However, when applied to V cmax the β factor reduced g sw and A net by 15 %-20 % (Exp. 2 vs. Exp. 3, Fig. 7). Applying the β factor to g 1 also reduced g sw and A net by 10 %-50 % (Exp. 1 vs. Exp. 3, Fig. 7). Unsurprisingly, comparing model results with β applied to all or no parameters showed the largest differences (30 %-80 %) (Exp. 2 vs. Exp. 4, Fig. 7). Default simulations with the β factor on g 0 and V cmax underestimated A net by 29 % and g sw by 26 % for the MED model. However, the results from simulations with no β effects or with β only applied to g 0 (Exps. 1 and 4) corresponded best to the observations, in which A net was only underestimated by 15 % and g sw by 9 % for the MED model ( Fig. 7a and  c). There was also a significant improvement of performance when the β effects were removed from the equation in the BWB model ( Fig. 7b and d).   (Table 2).

Advances in understanding model difference
We implemented the MED stomatal model in FATES and compared model projections of CO 2 and water vapor exchange to the existing BWB formulation. The two models diverged considerably in the responses of both leaf-level (g sw and A net , Figs. 1 and 2) and canopy-level (ET and GPP, Figs. 3 and 4) fluxes to a wide range of radiation, air temperature, VPD a , and CO 2 concentrations with the default stomatal slope parameter (g 1 ). When parameterization of g 1 was harmonized between the MED and BWB formulations, the difference was much smaller in responses to varying radiation, temperature, and CO 2 conditions but were markedly apparent at VPD a above 1.5 kPa.
Our analysis of the general model responses to synthetic climate forcing presents some advantages over previous evaluations. First, some studies found that different stomatal conductance models varied considerably in water-limited regions (Knauer et al., 2015;Morales et al., 2005) but were unable to attribute the difference to specific climate forcing as all factors, such as temperature and humidity, are closely related (Galbraith et al., 2010;Rowland et al., 2015). In their recent experimental study of a tropical forest, Smith et al. (2020) found that stomatal response to VPD a , rather than to T a , is the primary mechanism for high-temperature photosynthetic declines in tropical forests by separating the temperature effect and VPD a effect. This observation, along with our findings, highlights the need to improve the representation of stomatal conductance response to VPD a in models. Second, most previous modeling studies relied on evaluating model performance against benchmarks such as eddy covariance data and remote sensing products (De Kauwe et al., 2015), which were limited to the current climate conditions and ecoregions. To test model behaviors under all possible climate change scenarios, our studies designed simulations driven by a wide range of climate forcing data. Third, understanding model response to synthetic climate forcing (Figs. 1-4) is a powerful diagnostic tool because the model outputs can be evaluated in comparison to known and measurable physiological responses to environmental variation, such as radiation and CO 2 . The model outputs of GPP and ET also provide insight into how leaf-level responses influence the emergent ecosystem-scale responses, which is relevant for forecasting the responses of ecosystems and biomes to climate change. Fourth, by using the calibrated and default parameters to run the models, we were also able to separate effects of model structure (i.e., stomatal model choice) and parameterization (i.e., g 0 and g 1 ) on model differences.
As highlighted previously by Franks et al. (2017), the influence of parameterization dominated potential differences of g sw and ET due to model choice, further emphasizing the need to develop robust approaches to estimate g 1 and understand covariance with environmental drivers, such as soil moisture availability, and other leaf traits that may facilitate the use of trait-trait or trait-environment relationships to enable model parameterization (De Kauwe et al., 2015;Héroult et al., 2013;Lin et al., 2015;Wu et al., 2020). However, different g 1 values did not markedly change the magnitudes of A net and GPP, suggesting that the difference of g 1 propagates to the simulation of intercellular CO 2 first and finally to A net with attenuated effects. The structural difference is attributable to the different representation of humidity in the BWB and MED models (i.e., H s vs. VPD s ) and is consistent with previous studies (Rogers et al., 2017a;Knauer et al., 2015;Franks et al., 2017). g sw simulated by the power function of the MED model decreases hyperbolically, while that simulated by the linear function of the BWB model drops steeply. The nonlinear response of g sw to VPD a when using the MED model is supported by some observations (Marchin et al., 2016;Héroult et al., 2013;Wang et al., 2009;Domingues et al., 2014), but more measurements of leaf-level VPD a responses would be valuable. Our results suggest that when implemented in a dynamic vegetation demography model (FATES) the choice of stomatal model only has a small effect on projections of leaf-and canopy-level CO 2 and water vapor fluxes under conditions of VPD a below 1.5 kPa. Under higher VPD a , higher g sw values were simulated using the MED model compared with the BWB model and led to higher A net , ET, and GPP. This suggests that the MED formulation would predict tropical evergreen broadleaf forests to be more resistant to extreme atmospheric drought than with the BWB formulation. As the global surface temperature is projected to increase, the VPD a is also expected to increase (Ficklin and Novick, 2017;Yuan et al., 2019;Kolby Smith et al., 2016). Therefore, the difference between the two models under high VPD a conditions will lead to radically different ecosystem carbon and water dynamics under future climate change scenarios.

Model responses under simulated water stress
Our field campaign, which occurred during the 2016 ENSO event, enables us to evaluate model performance under various climate conditions, including extreme drought. Overall, simulations made with FATES with both g sw models captured the dynamics of measured upper canopy leaf-level fluxes well, confirming the utility of the current stomatal conductance models in LSMs for non-stressed conditions. However, at the peak of the dry season this underestimation of g sw and A net was notable and resulted in part from application of a soil water stress (β) factor used to modify leaf physiology in response to reduced soil moisture content. In FATES, the β factor affects g 0 and V cmax through an empirical modification. Experimental evidence about how physiological parameters change in response to soil water conditions is diverse. Some previous studies found that g 1 was relatively stable under water stress (Gimeno et al., 2016;De La Motte et al., 2020;Xu and Baldocchi, 2003). Other studies found a range of responses of g 1 to drought across different plant species (Miner and Bauerle, 2017;Zhou et al., 2013). For g 0 , it was reported to decrease under water stress (Miner and Bauerle, 2017;Misson et al., 2004) but also to show no response to drought (Barnard and Bauerle, 2013). Drought nearly universally lowered V cmax in plants (Zhou et al., 2013). However, some argued that effects of mild and moderate droughts on V cmax were negligible (Aranda et al., 2012;Bota et al., 2004;Cano et al., 2013), and others showed a range of responses resulting in a 10 %-25 % reduction in V cmax (Galmés et al., 2007;Grassi and Magnani, 2005;Keenan et al., 2010;Limousin et al., 2010;Wilson et al., 2000;Zait and Schwartz, 2018).
Despite previous extensive experimental studies of the β effects on plant physiological parameters, understanding of the results of applying β effects in models is still inadequate. The uncertainty of the β calculation is a major challenge. Based on the equations, the β factor is a function of soil water content, modified by parameters related to plant response, root distribution, and soil properties. Due to the lack of in situ measurements, we only used general parameters for the β factor in the simulations. Although soil moisture content was relatively well simulated (Fig. 6), root fraction and other soil properties were difficult to constrain due to scarce observations. In our study, by toggling on and off the β effects on stomatal and photosynthetic parameters, we were able to learn more about how the calculation of β influences model outputs. Overall, we found that the predictions of g sw and A net were closer to the measurements when the β factor was treated as 1 (i.e., no stress). Similar studies also found that the implementation of the β factor in the CLM overestimated the drought-related productivity loss compared with the observations, biased the transpiration rate, or lacked diurnal variability (Powell et al., 2013;Kennedy et al., 2019;Bonan et al., 2014). To improve models, further systematic evaluation of the β effects on photosynthetic capacity, stomatal conductance, and mesophyll conductance in LSMs is highly recommended (Egea et al., 2011;Vidale et al., 2021). More mechanistic approaches such as representation of hydraulic limitations and chemical signaling through abscisic acid (ABA) are emerging as promising ways to represent the plant response to drought in LSMs but come with significant added complexity (Verhoef and Egea, 2014;Sperry and Love, 2015;Kennedy et al., 2019).

Implications for evaluating model performance
Our FATES sensitivity analysis used synthetic meteorological forcing to enable us to isolate the impacts of individual abiotic drivers on model behaviors. By adjusting the specific humidity concurrently with air temperature, we were able to isolate the model response to changing air temperature from typically concurrent change in VPD a . We believe that our sensitivity analysis should be included as a routine approach for evaluating changes in model behavior during model development activities. Using similar simulations regularly during development would provide a powerful check on unexpected or unintended changes related to any changes in structure or parameters.
When the models are driven by synthetic climate forcing, special attention should be paid to the change in the environment conditions at the leaf surface, to which plants directly respond. In most LSMs, leaf surface temperature (T l ) is the balance of environmental drivers and leaf biophysical activities, and it is one of the most important variables regulating leaf biochemical responses such as photosynthesis and respiration (Kumarathunge et al., 2019;Leuning, 2002). But as T l is an emergent variable in FATES we could only control T a rather than T l in our sensitivity analysis simulations. In scenarios with changing radiation, although we have fixed T a , T l was also increasing (Fig. S4a), which resulted in slight decreasing trends of g sw and A net in response to radiation as T l exceeded the temperature optimum of g sw and A net (Figs. 1d and 2d). But the influence of T l change was limited for other response curves (Fig. S4).
Parameterization of g 0 has been shown to be critical for predicting ecosystem fluxes (De Kauwe et al., 2015;Barnard and Bauerle, 2013). However, there is little agreement on how to parameterize g 0 due to different definitions and measurement approaches for this parameter. Whether g 0 should be an intercept from data fitting, a minimum threshold when A net approaches zero, a nighttime g sw , or the cuticular conductance is still an active research topic Duursma et al., 2019;Lamour et al., 2022;Davidson et al., 2022). The slope parameter g 1 we used in the model was from Lin et al. (2015), estimated with the assumption that g 0 was zero. In our implementation, we not only included a nonzero g 0 in the numerical calculation of g sw , but also set a small positive value for g 0 to prevent g sw becoming zero or negative when A net approaches zero. In this way, the leaf stomatal resistance (i.e., the reverse of g sw ) will not become infinitive during the simulations. To understand how different g 0 values influence g sw and A net , we tested the sensitivity of g sw , A net , ET, and GPP to different g 0 values with our synthetic climate forcing listed in Table 1. In addition to the simulations with our default value (1000 µmol m −2 s −1 ), the g 0 was set zero or the commonly adopted value of 10 000 µmol m −2 s −1 (Sellers et al., 1996). A comparison of g 0 = 0 and g 0 = 1000 µmol mol −1 showed a very minor effect on the model response of g sw . Using the 10-fold larger estimate for g 0 (10 000 µmol mol −1 ) only resulted in a small effect on the magnitudes of g sw , A net , ET, and GPP (Figs. S5-S8).
For the model evaluation against site-level measurements, we found it is necessary to check the consistency of climate forcing used in models and that measured by the instruments (Héroult et al., 2013). In our study, the in situ measured T a and VPD a were higher than those recorded by a nearby meteorological station (Fig. 5a-d). The mismatch was partially related to the challenge of matching leaf chamber conditions with ambient conditions, avoiding condensation in the leaf chamber, or the use of a pump to move air across the leaf sur-face during gas exchange measurements. The slight deviation of modeled g sw and A net against measurements when soil was relatively wet (measured in February, March, and May) can be partly attributed to the mismatch of T a and VPD a used in the model compared with the in situ measurements. In this study we evaluated the inclusion of the MED model in FATES in a tropical forest. However, future efforts could include evaluations at sites of different ecosystem types and at regional and global scales for carbon and water cycles, particularly at sites where VPD a routinely rises above 1.5 kPa.

Conclusions
Implementing new plant physiological theories, such as the optimal stomatal conductance model, into dynamic vegetation models is crucial to keep the models up to date and to enable the exploration of new behaviors and capacities to understand potential ecosystem responses to global change. In this study, we added the optimality-based Medlyn model into the state-of-the-art dynamic vegetation model FATES as an alternative to the default Ball-Woodrow-Berry model and then tested model behaviors in response to key independent climate forcing. Our model evaluation demonstrated that the major difference between the two models was caused by the parameterization of the stomatal slope parameter (g 1 ). When parameters were harmonized, the potential for markedly different projections of water vapor and CO 2 fluxes between stomatal conductance models only occurred as VPD a rose above 1.5 kPa. We also compared model performance with gas exchange measurements from an evergreen tropical forest. Modeled CO 2 and water vapor fluxes in the dry season of a drought year were similar between models and closely matched observations, except at the peak of the dry season when a soil moisture correction factor was used to adjust physiological parameters. After removing this adjustment, projections for both models improved. Our study showed that the parameterization of g 1 and the application of the correction factor associated with decreasing soil moisture content are the key targets for improving model representation of CO 2 and water fluxes in tropical forests.