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 (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₂ 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₁). After harmonizing g₁ and holding the stomatal intercept parameter (g₀) 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. 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 CO₂ and water fluxes in tropical forests.

1 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₀), which is the conductance when photosynthesis (A) is zero, and the slope parameter (g₁) 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):

$$\begin{array}{}\text{(1)}& g_{\mathrm{sw}}=g_{\mathrm{0}}+g_{\mathrm{1}}{\displaystyle \frac{A_{\mathrm{net}}}{C_{\mathrm{s}}}}{H}_{\mathrm{s}}.\end{array}$$

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:

$$\begin{array}{}\text{(2)}& g_{\mathrm{sw}}=g_{\mathrm{0}}+\mathrm{1.6}\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}}{\sqrt{{\mathrm{VPD}}_{\mathrm{s}}}}}\right){\displaystyle \frac{A_{\mathrm{net}}}{C_{\mathrm{s}}}}.\end{array}$$

Although the functional form of the MED model is similar to the BWB model, g₁ 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₁ parameter has also been found to vary significantly across a wide range of different plant functional types (PFTs) and climate regions (Lin et al., 2015). Better representation of g_sw in LSMs requires efforts to improve the fidelity of g₁ parametrization by PFT. The g₁ parameter can be estimated through field measurement campaigns (Lin et al., 2015; 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₁ 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₂ 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₂-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₀ (e.g., CLM, Lawrence et al., 2019), g₁ (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 leaf-level and canopy-level CO₂ 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?

2 Methods

2.1 Implementation of the Medlyn model into FATES

In FATES, leaf-level photosynthesis (A) in C₃ 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 RuBP-carboxylase-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₂, and energy cycles in forests. It not only controls the water and CO₂ 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₀ 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):

$$\begin{array}{}\text{(3)}& \mathit{\beta }={\sum }_{j=\mathrm{1}}^{nj}{w}_j{r}_j,\end{array}$$

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).

2.2 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 (Wu et al., 2020; 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.

2.3 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₁ 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₁ value of 8 (unitless) for the BET tree PFT in FATES. In our MED-default setup, the MED model was parameterized with g₁ 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₁ value that was harmonized with the BWB model in FATES, which was abbreviated as MED-B. Here, we assumed g₁ 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 BWB-equivalent g₁=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, g₀ was fixed at 1000 µmol m⁻² s⁻¹.

$$\begin{array}{}\text{(4)}& g_{\mathrm{1}\mathrm{b}}={\displaystyle \frac{\mathrm{1.6}}{H_{\mathrm{a}}}}\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}\mathrm{m}}}{\sqrt{{\mathrm{VPD}}_{\mathrm{a}}}}}\right)\end{array}$$

The FATES model is driven by half-hourly longwave radiation, shortwave radiation, air temperature, specific humidity, precipitation, surface pressure, wind speed, and atmospheric CO₂ 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₂ concentration (CO₂), 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⁻¹, surface pressure to 99 626 Pa, wind speed to 4.8 m s⁻¹, and longwave radiation to 407.4 W m⁻² 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.

Table 1 Scenario setting for the sensitivity simulations.

2.4 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₂ concentration, and temperature surrounding the leaves were closely matched to the ambient conditions. Using this dataset, the g₁ 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₁ 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₁ was estimated for BWB and MED models based on the same measurements, g₁ 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₁ 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⁻² s⁻¹ 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₂ 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₂ concentration inside the leaf chamber of the gas exchange equipment.

2.5 Drought effects on physiological parameters

In FATES, a soil water stress factor (β) is used to adjust g₀ 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₁ 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₁ and V_cmax were set as the averages across all the species measured for the BET PFT.

Table 2 Model simulations for studying soil water stress effects on physiological parameters in FATES.

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.

3 Results

3.1 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₁) and MED-default (MED model with default g₁) simulations had similar shapes (Figs. 1 and 2, blue and black lines). The MED-default 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).

Figure 1 The responses of stomatal conductance (g_sw) to scenarios with (a) radiation, (b) CO₂, (c) VPD_a, and (d) T_a for the three model setups: BWB (blue), MED-B (red), and MED-default (black).

Figure 2 The responses of net photosynthesis (A_net) to scenarios with (a) radiation, (b) CO₂, (c) VPD_a, and (d) T_a for the three model setups: BWB (blue), MED-B (red), and MED-default (black).

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₂ (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₁. With a harmonized parameterization of g₁ 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₁ for the MED model (MED-default) 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).

Figure 3 The responses of evapotranspiration (ET) to scenarios with (a) radiation, (b) CO₂, (c) VPD_a, and (d) T_a for the three model setups: BWB (blue), MED-B (red), and MED-default (black).

Figure 4 The responses of gross primary productivity (GPP) to scenarios with (a) radiation, (b) CO₂, (c) VPD_a, and (d) T_a for the three model setups: BWB (blue), MED-B (red), and MED-default (black).

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 leaf-level physiology.

3.2 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₂ 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₁ across different species, we ran a series of FATES simulations driven by meteorological forcing data with different g₁ values. These experiments showed that FATES A_net and g_sw were sensitive to different g₁ values for both model formulations (Fig. 5e–l). The MED model ensemble results for A_net and g_sw with different g₁ 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).

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₂ concentration, and gold for H_a. In panels (e)–(l) shaded areas represent the range of FATES model ensemble results with different measured g₁ 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.

3.3 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₀, g₁, 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₀, g₁, and V_cmax (Table 2, Fig. 7). The results from ​​​​​​​Exp. 1 and Exp. 4 showed high overlap, indicating that considering the β effect on g₀ 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₁ 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₀ 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₀ (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).

Figure 6 Annual cycle of the modeled volumetric soil water content (blue line) and corresponding soil water stress (β) factor (green line) from the FATES simulation as well as ERA5 reanalysis soil water content data (orange line) at the San Lorenzo field site in 2016. The soil water content data are means across all soil layers. For the β factor, 1 represents fully saturated soil, while 0 represents very dry soil. Vertical gray lines indicate the four campaign dates in 2016.

Figure 7 Comparison between the model outputs and measurements on 21 April 2016 for (a) the MED model net photosynthesis (A_net), (b) the BWB model A_net, (c) the MED model stomatal conductance (g_sw), and (d) the BWB model g_sw with different soil water stress effects on parameters in FATES (Table 2).

4 Discussion

4.1 Advances in understanding model difference

We implemented the MED stomatal model in FATES and compared model projections of CO₂ 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₂ concentrations with the default stomatal slope parameter (g₁). When parameterization of g₁ was harmonized between the MED and BWB formulations, the difference was much smaller in responses to varying radiation, temperature, and CO₂ 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₂. 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₀ and g₁) 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₁ 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₁ values did not markedly change the magnitudes of A_net and GPP, suggesting that the difference of g₁ propagates to the simulation of intercellular CO₂ 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₂ 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.

4.2 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₀ 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₁ 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₁ to drought across different plant species (Miner and Bauerle, 2017; Zhou et al., 2013). For g₀, 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; Misson 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).

4.3 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₀ 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₀ due to different definitions and measurement approaches for this parameter. Whether g₀ 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 (Lombardozzi et al., 2017; Duursma et al., 2019; Lamour et al., 2022; Davidson et al., 2022). The slope parameter g₁ we used in the model was from Lin et al. (2015), estimated with the assumption that g₀ was zero. In our implementation, we not only included a nonzero g₀ in the numerical calculation of g_sw, but also set a small positive value for g₀ 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₀ values influence g_sw and A_net, we tested the sensitivity of g_sw, A_net, ET, and GPP to different g₀ values with our synthetic climate forcing listed in Table 1. In addition to the simulations with our default value (1000 µmol m⁻² s⁻¹), the g₀ was set zero or the commonly adopted value of 10 000 µmol m⁻² s⁻¹ (Sellers et al., 1996). A comparison of g₀=0 and g₀=1000 µmol mol⁻¹ showed a very minor effect on the model response of g_sw. Using the 10-fold larger estimate for g₀ (10 000 µmol mol⁻¹) 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 surface 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.

5 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₁). When parameters were harmonized, the potential for markedly different projections of water vapor and CO₂ 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₂ 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₁ and the application of the correction factor associated with decreasing soil moisture content are the key targets for improving model representation of CO₂ and water fluxes in tropical forests.