Evaluating the vegetation–atmosphere coupling strength of ORCHIDEE land surface model (v7266)

Plant transpiration dominates terrestrial latent heat fluxes (LE) and plays a central role in regulating the water cycle and land surface energy budget. However, Earth system models (ESMs) currently disagree strongly on the amount of transpiration, and thus LE, leading to large uncertainties in simulating future climate. Therefore, it is crucial to correctly represent the mechanisms controlling the transpiration in models. At the leaf scale, transpiration is controlled by stomatal regulation, and at the canopy scale, through turbulence, which is a function of canopy structure and wind. The coupling of vegetation to the atmosphere can be characterized by the coefficient Ω. A value of Ω→0 implies a strong coupling of vegetation and the atmosphere, leaving a dominant role to stomatal conductance in regulating water (H₂O) and carbon dioxide (CO₂) fluxes, while Ω→1 implies a complete decoupling of leaves from the atmosphere, i.e., the transfer of H₂O and CO₂ is limited by aerodynamic transport. In this study, we investigated how well the land surface model (LSM) Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) (v7266) simulates the coupling of vegetation to the atmosphere by using empirical daily estimates of Ω derived from flux measurements from 90 FLUXNET sites. Our results show that ORCHIDEE generally captures the Ω in forest vegetation types (0.27 ± 0.12) compared with observation (0.26 ± 0.09) but underestimates Ω in grasslands (GRA) and croplands (CRO) (0.25 ± 0.15 for model, 0.33 ± 0.17 for observation). The good model performance in forests is due to compensation of biases in surface conductance (Gs) and aerodynamic conductance (Ga). Calibration of key parameters controlling the dependence of the stomatal conductance to the water vapor deficit (VPD) improves the simulated Gs and Ω estimates in grasslands and croplands (0.28 ± 0.20). To assess the underlying controls of Ω, we applied random forest (RF) models to both simulated and observation-based Ω. We found that large observed Ω are associated with periods of low wind speed, high temperature and low VPD; it is also related to sites with large leaf area index (LAI) and/or short vegetation. The RF models applied to ORCHIDEE output generally agree with this pattern. However, we found that the ORCHIDEE underestimated the sensitivity of Ω to VPD when the VPD is high, overestimated the impact of the LAI on Ω, and did not correctly simulate the temperature dependence of Ω when temperature is high. Our results highlight the importance of observational constraints on simulating the vegetation–atmosphere coupling strength, which can help to improve predictive accuracy of water fluxes in Earth system models.

1 Introduction

Accurately representing the land–atmosphere interactions in Earth system models (ESMs) is crucial for analyzing climate variability and projecting climate change (Claussen, 1998; Goldberg and Bernhofer, 2001; Zhu et al., 2017). Among the key interactions, the exchange of latent heat (LE) between the land surface and the atmosphere is one of the most important processes (Trenberth et al., 2009; IPCC, 2014). LE is contributed by several sources, including evaporation from bare soil and canopy interception, vegetation transpiration, snow and ice sublimation (Chapin et al., 2011). In these sources, transpiration has the largest contribution (Jasechko et al., 2013; Wei et al., 2017; Li et al., 2019), but is massively uncertain across models (Stoy et al., 2019), leading to considerable uncertainty in LE simulation in current ESMs (Wild, 2020). The large uncertainties in current transpiration and LE simulations can further result in difficulties in constraining soil moisture and the carbon cycle (Humphrey et al., 2021). Therefore, there is a need to evaluate and improve the simulation of transpiration and LE in ESMs.

The LE parameterization in ESMs is based on Fick's law, using the conductance, or 1$/$resistance of water vapor between vegetation and atmosphere (Bonan, 2019). This conductance is the result of several processes such as stomatal opening, boundary layer turbulence, soil-to-air evaporative resistance; it is thus affected by multiple factors including plant physiology, vegetation structure, vapor pressure deficit (VPD), temperature, net radiation, soil moisture, etc (Igarashi et al., 2015; F. Zhang et al., 2018; Veste et al., 2020). Currently, we can observe total LE at the site scale (i.e., FLUXNET), but we are unable to disentangle the relative contribution of different processes. The complexity of conductance and the lack of process-level observations lead to difficulties in detailed evaluation of the vegetation–atmosphere water exchanges in ESMs based on the underlying processes. As a result, accurately capturing the regulation of LE by biotic and abiotic factors remains a key challenge for the land surface modeling community (Mueller et al., 2013; De Kauwe et al., 2017; Stoy et al., 2019).

An early attempt to quantify the contribution of different conductance processes was made by Jarvis and McNaughton (1986), who developed a metric commonly referred to as the decoupling coefficient, Ω, to describe whether vegetation transpiration is mainly controlled by stomatal or aerodynamic processes. The calculation of Ω is based on the ratio between aerodynamic and stomatal conductance (see “Data and methods” section). At the limit, Ω=0 denotes perfect coupling between vegetation and atmosphere, i.e., the transpiration is entirely regulated by stomata, while Ω=1 denotes complete decoupling, i.e., transpiration is driven entirely by boundary layer turbulence. The concept of Ω can be used at scales from leaf to regional level, and for different fluxes from transpiration only to the total evapotranspiration (ET; e.g., Peng et al., 2019). Because evapotranspiration includes water fluxes from not only leaf but also other surfaces, the stomatal conductance needs to be replaced by a surface conductance (Gs) which integrates all conductances at different surfaces in the evapotranspiration Ω calculation.

During the last decades, the number of eddy covariance flux measurements has grown rapidly. Quantification of Ω at site level from eddy covariance flux measurements offers insights into how different vegetation types control turbulent fluxes as a function of their phenology and stomatal physiology during the growing and non-growing seasons (De Kauwe et al., 2017; Goldberg and Bernhofer, 2001). These observation-based Ω provide valuable information to evaluate ESMs on how well they capture the controls of LE. Using this estimate, De Kauwe et al. (2013) found that one of the principal reasons for disagreement among simulated transpiration responses to elevated CO₂ is the differences in the degree of coupling between vegetation and the atmosphere.

The Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) land surface model (LSM) is one of the widely used models for simulating carbon, energy and water budget of terrestrial ecosystems (e.g., Zhang et al., 2021; Schrapffer et al., 2020). ORCHIDEE and the ESM Institute Pierre-Simon Laplace climate model (IPSLCM), which has ORCHIDEE as the land surface module, have participated in various model intercomparison projects including TRENDY, Coupled Model Intercomparison Project (CMIP), etc. In spite of its wide usage, the LE of the ORCHIDEE LSM remains simply calibrated and evaluated against the total evapotranspiration observations (Bastrikov et al., 2018), without considering the detailed processes. A recent study showed that the ORCHIDEE version used in CMIP6 still has biases in LE, especially in tropical regions (Tafasca et al., 2020). However, it remains unclear how the biases happened and which processes need to be improved to better simulate the fluxes. To solve this problem, in this study, we used the Ω dataset derived from eddy covariance data from 90 sites (De Kauwe et al., 2017) to evaluate the vegetation–atmosphere coupling strength of the land surface model ORCHIDEE 2.2 (v7266). We tested whether the calibration of the stomatal response to atmospheric dryness, or using observed canopy height, can improve the simulation of coupling strength. Further, we used random forest (RF) models to investigating the biotic and abiotic factors affecting the coupling strength. The methodology presented here is generic enough to be applied for the benchmarking of other LSMs. The objectives of this study are as follows: (1) benchmark ORCHIDEE using Ω estimated from FLUXNET observations; (2)investigate how different factors affect Ω in the observations and (3) investigate whether ORCHIDEE correctly captured the driving factors.

2 Data and methods

2.1 ORCHIDEE model

We use the ORCHIDEE 2.2 (v7266) land surface model in this study. This model version is the latest version participating in the CMIP6 project under coupled configuration to the atmospheric circulation model in the IPSL–CM6A–LR ESM (Boucher et al., 2020). The ORCHIDEE model consists of three interactive sub-modules (Krinner et al., 2005). The SECHIBA module parameterizes the land surface energy and water balance (Ducoudré et al., 1993). The STOMATE module deals with phenology (Botta et al., 2000) and carbon fluxes of terrestrial ecosystems (Viovy, 1996). The LPJ dynamic vegetation module simulates the dynamics of vegetation (Sitch et al., 2003). In this study, the dynamic vegetation module is turned off because the vegetation types are prescribed at each site.

ORCHIDEE simulates LE by considering plant transpiration, bare soil evaporation, sublimation, floodplain evaporation, and evaporation from canopy water interception. Since this study focuses on the vegetation–atmosphere coupling strength for transpiration and because the data to evaluate this model have been filtered to represent the transpiration (De Kauwe et al., 2017), we only introduce the parameterization of conductance relating to transpiration in ORCHIDEE here.

The stomatal conductance (gs, $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{bar}}^{-\mathrm{1}}$) is calculated in the photosynthesis module which couples the leaf-level photosynthesis and stomatal conductance based on (Yin and Struik, 2009):

$$\begin{array}{}\text{(1)}& \text{gs}=\text{g0}+{\displaystyle \frac{A+R_{\mathrm{d}}}{C_{\mathrm{i}}-C_{\mathrm{i}}^{\ast }}}{f}_{\text{vpd}},\end{array}$$

where g0 is the stomatal conductance when the irradiance is zero ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{bar}}^{-\mathrm{1}}$); A is the rate of CO₂ assimilation ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$); R_d is the dark respiration ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$); C_i is the intercellular CO₂ partial pressure (µbar); $C_{\mathrm{i}}^{\ast }$ is the C_i-based CO₂ compensation point (µbar) in the absence of R_d; and f_vpd is the function for the effect of vapor pressure deficit (VPD, kPa) on stomatal conductance, calculated as

$$\begin{array}{}\text{(2)}& f_{\text{vpd}}={\displaystyle \frac{\mathrm{1}}{\left[\frac{\mathrm{1}}{a_{\mathrm{1}}-b_{\mathrm{1}}\text{VPD}}-\mathrm{1}\right]}},\end{array}$$

where a₁ and b₁ are empirical parameters depending on vegetation type (Fig. S1 in the Supplement). This equation shows that a higher VPD will induce stomatal closure and decrease gs.

The canopy-level stomatal conductance is calculated by integrating gs across all leaves in the canopy.

The aerodynamic conductance (Ga, $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) formulation in ORCHIDEE is

$$\begin{array}{}\text{(3)}& \text{Ga}={\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}{u}_{\mathrm{a}}}{\left[\ln \left(\frac{z_a-d}{z_{\text{0m}}}\right)\ln \left(\frac{z_a-d}{z_{\text{0h}}}\right)\right]}}\text{ps}/\left(RT\right),\end{array}$$

where z_a is the height of the wind measurement, d is the displacement height (i.e., the height at which the wind speed would go to zero), calculated as 0.66 of average canopy height. The wind speed is u_a (m s⁻¹) and κ is the von Karman's constant. Air pressure and temperature are denoted as ps, T, R is the universal gas constant, and z_0m and z_0h are respectively the roughness heights (m) for momentum and heat transfer estimated following Su et al. (2001) and Ershadi et al. (2015) using canopy height (z) and LAI:

$$\begin{array}{}\text{(4)}& z_{\text{0m}}=(z-d)e^{-\frac{\mathit{\kappa }}{\mathit{\eta }}},\end{array}$$

where

$$\begin{array}{}\text{(5)}& \mathit{\eta }=\mathrm{0.32}-\mathrm{0.264}{e}^{-\mathrm{3.02}\text{LAI}}.\end{array}$$

z_0h is estimated using z_0m:

$$\begin{array}{}\text{(6)}& z_{\text{0h}}={\displaystyle \frac{z_{\text{0m}}}{e^{\mathit{\kappa }{B}^{-\mathrm{1}}}}},\end{array}$$

where B is the Stanton number; κB⁻¹ is estimated following Su et al. (2001):

$$\begin{array}{}\text{(7)}& \mathit{\kappa }{B}^{-\mathrm{1}}={\displaystyle \frac{\mathit{\kappa }\text{Cd}}{\mathrm{4}\text{Ct}\mathit{\eta }(\mathrm{1}-e^{-\frac{\text{nec}}{\mathrm{2}}})}}{\text{fc}}^{\mathrm{2}}+\mathrm{2}\text{fc}\text{fs}{\displaystyle \frac{\mathit{\kappa }\mathit{\eta }\frac{z_{\text{0m}}}{z}}{C_{\mathrm{t}}^{\ast }}}+\mathit{\kappa }{B}_{\mathrm{s}}^{-\mathrm{1}}{\text{fs}}^{\mathrm{2}},\end{array}$$

where Cd and Ct are drag and heat transfer coefficients of leaves, nec is within the canopy wind profile extinction coefficient, calculated as $\text{nec}=\text{CdLAI}/(\mathrm{2}{\mathit{\eta }}^{\mathrm{2}}$); fc and fs are the fraction of canopy and bare soil, $C_{\mathrm{t}}^{\ast }$ is the heat transfer coefficient of soil; B_s is the Stanton number for bare soil, with $\mathit{\kappa }{B}_{\mathrm{s}}^{-\mathrm{1}}$ estimated following Brutsaert (1999):

$$\begin{array}{}\text{(8)}& \mathit{\kappa }{B}_{\mathrm{s}}^{-\mathrm{1}}=\mathrm{2.46}{\mathit{\text{Re}}}_{\ast }^{\frac{\mathrm{1}}{\mathrm{4}}}-\ln \left(\mathrm{7.4}\right),\end{array}$$

where Re_∗ is the Reynolds number.

2.2 FLUXNET data and empirical calculation of Ω

The empirical Ω reference is derived from the FLUXNET 2015 dataset (Pastorello et al., 2020). This dataset collects eddy covariance measurements of heat and water fluxes, as well as the corresponding meteorological variables above the vegetation canopy in sites over the world and across different plant functional types (PFTs). The detailed information of the flux sites used can be found in Table S1 in the Supplement.

The calculation of Ω was firstly introduced by Jarvis and McNaughton (1986), using the formulation:

$$\begin{array}{}\text{(9)}& \mathrm{\Omega }={\displaystyle \frac{\mathrm{1}+\mathit{ϵ}}{\mathrm{1}+\mathit{ϵ}+\frac{\mathrm{Ga}}{\mathrm{Gs}}}},\end{array}$$

where $\mathit{ϵ}=s/\mathit{\gamma }$; s is the slope of the saturation vapor pressure curve with air temperature (Pa K⁻¹); γ is the psychrometric constant (Pa K⁻¹). It should be noted that the conductance (Ga, Gs) used for Ω calculation depends on the scale of interest. At the scale larger than a leaf, if other water vapor fluxes besides transpiration (e.g., soil evaporation) have significant contribution to LE, Gs must also include such contribution. In such cases, the synthesized Gs was sometimes referred to as surface conductance (Peng et al., 2019). To be accurate, we use the term surface conductance for Gs hereafter to match our scale.

There remains no direct observation of Ga and Gs at flux sites. De Kauwe et al. (2017) developed an empirical method to estimate the two terms. In this method, Ga was estimated as an empirical equation using wind speed and friction velocity (Thom, 1972), and Gs ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) was estimated using the inverted Penman–Monteith equation with measured evapotranspiration (ET, in $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) flux:

$$\begin{array}{}\text{(10)}& \text{Gs}={\displaystyle \frac{\mathrm{Ga}\mathit{\gamma }\mathit{\lambda }\text{ET}}{s(R_{\mathrm{net}}-G)-(s+\mathit{\gamma })\mathit{\lambda }\text{ET}+\mathrm{Ga}{M}_{\mathrm{a}}c\text{VPD}}},\end{array}$$

where λ is the latent heat of vaporization (J mol⁻¹); VPD (Pa) is the vapor pressure deficit; R_net (W m⁻²) is the net radiation flux; G (W m⁻²) is the soil heat flux; M_a (kg mol⁻¹) is molar mass of air, and c is the heat capacity of air ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$).

In this study, Ga, Gs and Ω from the dataset of De Kauwe et al. (2017) are used as the reference to evaluate the ORCHIDEE LSM.

2.3 Simulation setup and modeled Ω calculation

The site simulations with ORCHIDEE are forced with observed meteorology in the FLUXNET 2015 dataset (Pastorello et al., 2020). The variables include half-hourly time series of air temperature (K), surface pressure (Pa), specific humidity (kg kg⁻¹), northerly and easterly wind speeds (m s⁻¹), short-wave downward radiation (W m⁻²), long-wave downward radiation (W m⁻²), rainfall ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and snowfall ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Gaps in the FLUXNET meteorology data are filled following Vuichard and Papale (2015). The PFT classification of FLUXNET is different from the one used in ORCHIDEE. To let ORCHIDEE simulate LE and the conductances without bias, we used a combination of ORCHIDEE PFTs to represent the vegetation type at each site (Table S1).

Three simulations are performed at each site (Fig. 1). The first simulation named Ctrl uses the default configuration and parameters as used in CMIP6 and TRENDY experiments. The second simulation named Clb_gs uses the same configuration as Ctrl but changes the empirical parameters in Eq. (2). New values for a₁ and b₁ are obtained by constraining the modeled formulation of conductance against a global database of leaf-level observations of stomatal conductance from Lin et al. (2015) for different PFTs (see Table S2 and Fig. S1 in the Supplement). Finally, because the ORCHIDEE model prescribes canopy height for each PFT (Table S3 in the Supplement), which may cause biases in Ga, we performed a last simulation referred to as Clb_ht. The latter simulation also uses the Ctrl configuration but the default canopy height parameters for each PFT are replaced by the canopy height observed at each site. In all the simulations, we kept the distance between measurement height and canopy height consistent with the observations, to ensure unbiased estimates of aerodynamic conductance in the model. Because canopy height and measurement height are required in the last simulation, in this study, we only used 90 sites out of the flux sites in the FLUXNET 2015 dataset where we found both heights.

Figure 1 Workflow of this study.

Although De Kauwe et al. (2017) excluded time steps with precipitation and the subsequent 48 half hours to have the LE mainly contributed by transpiration, referring to Gs as “stomatal conductance” in their paper, we still need to keep in mind that the Gs calculated in this way may also contain contributions from several other processes. It includes the conductance related to bare soil evaporation and the one related to water transport in the leaf boundary layer, in addition to the stomatal conductance integrated over the entire canopy. Therefore, it is more a “surface conductance” than a “stomatal conductance”. To be consistent with the observation-based dataset, we did not use the integrated canopy level stomatal conductance from ORCHIDEE output to calculate Ω. Instead, Gs is diagnosed using the ORCHIDEE output evapotranspiration, R_net and G following Eq. (5).

2.4 Leaf area index data

Because leaf area is an important factor affecting both aerodynamic and surface conductance, it is necessary to take leaf area into consideration when explaining the decoupling coefficient. However, instantaneous leaf area information is not available at most of the flux sites. To match the space and time of observation-based Ω, we extracted the leaf area index (LAI) from the 8 d 500 m dataset, MOD15A2H, derived from the space-borne MODIS observations (Myneni et al., 2015). This LAI dataset shows good consistency with in situ observations (Xu et al., 2018). The LAI for a given date is interpolated by averaging the nearest two high-quality LAI observations from the 8 d time series. For the simulated Ω, we used the LAI from the simulations for analyses to keep consistency between Ω and LAI.

2.5 Analyses

To be comparable with the observation-based Ω dataset, we first used the same criteria to screen the model outputs as De Kauwe et al. (2017), i.e., (1) only the 3 most productive months, to account for the different timing of summer in the Northern (June, July, August) and Southern (December, January, February) hemispheres are included in the study. This is to maximize the role of transpiration in Ω versus bare soil evaporation in the growing season. (2) Only daytime data from 08:00 to 16:00 (local solar time) are used. (3) Time steps during precipitation or within 2 d after precipitation are excluded. Because the 30 min Ω is very noisy, to reduce the noise in data, we used the daytime average of Ω and explanatory variables in all later analyses.

The decoupling coefficient Ω is affected by multiple factors and the relationships between Ω and different factors are often nonlinear. To characterize these relationships, we constructed random forest (RF) models for each of the observation-/simulation-based daily Ω. The goal here is to diagnose the main explanatory variables from the RFs in the observations/simulations, and to gain insights about the model over-/under-representation of their relative importance. The explanatory variables used in the RF models include wind speed, air temperature (Tair), VPD, net radiation (Rnet), LAI, canopy height and PFT. For each model, 90 % of the data are randomly sampled for training and the leftover 10 % are used for testing whether there is overfitting in the RF models (Fig. S2 in the Supplement).

To visualize the role of each factor in the complex RF model, we calculated SHAP (SHapley Additive exPlanations) values. A SHAP value is an index based on the classic Shapley values from the game theory (Lundberg and Lee, 2017). For each daily sample, SHAP calculates the expectation of contribution of each factor to deviate the sample value from the average of all samples. An example explaining the SHAP values can be found in Fig. S3 in the Supplement. Investigating the dependence of the SHAP value to the factor value tells us how this factor affects Ω. Moreover, by averaging the absolute values of the SHAP of 1 factor from all samples, we can get the importance of the factor in the RF model.

The workflow of the simulations and analyses can be found in Fig. 1.

3 Results

3.1 The performance of the ORCHIDEE model

The average growing season daytime Ω estimated from observations and from the ORCHIDEE outputs are shown in Fig. 2. A remarkable difference in the decoupling coefficient is found among plant functional types. According to the observation-based estimation (De Kauwe et al., 2017), the short vegetation types including grasslands (GRA) and croplands (CRO) are generally more decoupled from the atmosphere than forests, with the median values of Ω over sites of 0.31 and 0.38. In forest vegetation types, the evergreen forests (median Ω=0.26–0.35) are more decoupled with the atmosphere than deciduous forests (median Ω=0.16). The wetlands (WET) in observation show a strong decoupling (median Ω=0.42). Considering the large evaporation from open water in this vegetation type, the strong decoupling is not surprising. Besides the difference among vegetation types, we also find large variability in Ω within each type, especially for GRA and CRO (Table S4 in the Supplement).

Figure 2 Box plots of site mean Ω observations (Flux) and different simulations; n indicates the number of sites in each PFT group: EBF – evergreen broadleaf forests; DBF – deciduous broadleaf forests; ENF – evergreen needleleaf forests; MF – mixed forests; SAV – savannahs; WSA – woody savannahs; CSH – closed shrublands; OSH – open shrublands; WET – wetlands; GRA – grasslands; CRO – croplands.

Compared with observations, ORCHIDEE Ctrl simulations show similar median Ω in forests and CRO (Fig. 2 and Table S4). However, in GRA, the Ctrl median Ω (0.15) is much smaller compared to observations (0.31), implying a greater stomatal control in the model than the observations on GRA transpiration. This bias is not contributed by a few outlier sites but by a systematic underestimation of Ω at most of the GRA sites. For WET, ORCHIDEE also shows a significant underestimation of Ω (Fig. 2). This could be due to the lack of wetland PFTs and the corresponding open water in the ORCHIDEE model (Table S3). Despite the biases in GRA and WET, the observed differences in Ω among vegetation types are to a larger degree well reproduced (Fig. 2). The strongest decoupling is found in CRO and deciduous broadleaf forest (DBF), and the evergreen needleleaf forests (ENF) are more coupled than deciduous broadleaf forests.

By calibrating stomatal conductance (VPD dependence parameters leading to the Clb_gs simulation), we obtained Ω estimations closer to observations in short vegetation types (CRO and GRA) than Ctrl (Fig. 2). However, the median Ω estimation for most forest types is degraded after the gs “calibration”, with the Ω being more overestimated in DBF, ENF and mixed forests (MF). In contrast to the large impact from the calibration of stomatal conductance, prescribing realistic canopy height to the model leads to minor changes in Ω (Fig. 2).

In order to understand the reasons for differences in Ω between observations and the ORCHIDEE model, we also look into its components Ga and Gs (Fig. 3). Compared to observations, both Ga and Gs are underestimated in Ctrl. For Ga, the underestimation from the model is ∼1.0 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in forest types and ∼0.4 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in GRA and CRO. Calibrating stomatal conductance (Clb_gs) or prescribing the observed canopy height to the model (Clb_ht) both have a small impact on Ga. For Gs, using the new parameters for stomatal conductance (Clb_gs) can generally correct the Gs bias in DBF, ENF and MF, and improved Gs in GRA and CRO than Ctrl. Although Clb_gs has improved the Gs simulation compared with Ctrl, it does not result in an improvement of Ω and LE simulation, implying a compensation of biases in Ga and Gs in the current ORCHIDEE model.

Figure 3 Same as Fig. 2 but for (a) aerodynamic conductance, (b) surface conductance and (c) latent heat.

3.2 Factors controlling the decoupling coefficient

To better understand the underlying drivers of the variability in decoupling, we separated the importance of hypothesized drivers of decoupling coefficient in RF models using SHAP values (Fig. 4a). Among all the factors, the observation-based RF results show that the variation of Ω is mainly contributed by the variation of VPD, followed by PFT, with each of them having a SHAP value of ∼0.06, i.e., the variation of the factor contributes on average 0.06 of the deviation of Ω (absolute value) from the average of all samples. The other factors show relatively small importance to Ω, with SHAP values smaller than 0.03. Compared to observations, the ORCHIDEE Ω variation is also strongly contributed by VPD. However, contrary to the strong PFT impact found in observations, the modeled Ω is strongly affected by LAI. In Ctrl, the SHAP value of LAI is 0.09, which is much higher than the observation. The calibration of gs increased this value to 0.14. In contrast to the strong impact of LAI, all the modeled Ω show a much smaller contribution from PFT than in observations. It is also notable that the impact of air temperature on Ω is also much smaller in ORCHIDEE simulations than in observations.

Figure 4 Importance of different factors based on absolute SHAP values of Ω (the expectance of factor-induced deviation of Ω from the averages of all samples). (a) Random forest model built by data from all PFTs. (b) Model using only forest data. (c) Model using only grassland and cropland data.

To further understand the differences between tall and short vegetation, we trained random forest models using only forests (evergreen broadleaf forests (EBF), DBF, ENF and MF) and only short vegetation (GRA and CRO) observations/simulations. In forests, the SHAP value of VPD is comparable in the observations and ORCHIDEE simulations, while the LAI SHAP value is strongly overestimated and the canopy height SHAP value is slightly underestimated by the model. For short vegetation, a strong overestimation of the SHAP of LAI is also confirmed in ORCHIDEE. However, for the other factors (Tair, Rnet, VPD and height), the SHAP values are underestimated. It is notable that the SHAP values for VPD in ORCHIDEE is only 60 % of the estimation in observation, probably indicating a strong underestimation of water stress on Ω in short vegetation.

Figure 5 summarizes how different factors affect Ω in each of the observations/simulations in random forest models. The responses of Ω to most factors are generally consistent in observations and simulations. According to all of the random forest models, the vegetation is more decoupled, or having a larger Ω, under conditions with low wind speed, low VPD and large LAI. Also, both observations and simulations agree that GRA and CRO are more decoupled from the atmosphere than the other PFTs. However, for Tair and Rnet, ORCHIDEE does not capture the observed dependence correctly. In observations, a remarkable positive Tair dependence is found, with higher temperature tending to result in higher Ω. While in simulations, temperature shows a very small impact on Ω. The dependence of Ω on Rnet is similar to that of Tair in observation, but only the Clb_gs simulation captured this dependence correctly. Finally, to our surprise, we did not find Ω to strongly depend on canopy height in both observations and simulations. Although the highest canopy tends to have positive SHAP values, the range of SHAP values for smaller height levels is very large with both positive and negative.

Figure 5 Beeswarm plots showing the dependence of Ω SHAP values to different factors. For each data point, the percentile of the factor's value in all samples is shown in color. The SHAP value, or contribution of this factor value to deviate the daytime Ω from the average Ω of all samples, is shown in the x axis. In each subplot, data points at a certain SHAP value level are sorted by the factor percentile (i.e., vertical gradient indicates the distribution of factor values in the data). (a) Based on the observation dataset, (b), (c) and (d) are for Ctrl, Clb_gs and Clb_ht simulations, respectively.

A comparison of all individual controlling factors between the observations and the ORCHIDEE simulations is shown in Fig. 6. The dependence of Ω on wind speed generally has similar patterns in observations and in ORCHIDEE. Similar patterns are also found in Ga and Gs between simulations and observations at wind speeds larger than 1 m s⁻¹. In observations, we found positive SHAP values of wind speed at wind speeds smaller than 1 m s⁻¹. This might be due to coincidence because low wind speed will cause large uncertainty in the eddy covariance measurements and there are very few valid observation-based Ω available at low wind speed.

Figure 6 Dependence of Ω (top), Ga (middle) and Gs (bottom) SHAP values on different factors (in order from left to right: wind speed, air temperature, VPD, net radiation, LAI and canopy height). The colors indicate observation or simulations. Red: observation-based dataset, green: Ctrl, blue: Clb_gs, brown: Clb_ht. The shaded dots show the distribution of SHAP values in sample.

The observed dependence of Ω on Tair is not captured by ORCHIDEE. Observations indicate an increase of Ω when Tair is lower than 30^∘, and a slight decrease at a higher temperature, while ORCHIDEE simulations show a much smaller impact from Tair. This model bias is caused by differences in the relationships of Gs on Tair at a high temperature. A strong decline of the Gs SHAP values is found when the Tair is more than 20^∘ in ORCHIDEE, while the observations show a slight increase of Gs SHAP values at the same temperature. This difference probably indicates an underestimation of optimal temperature for photosynthesis in ORCHIDEE in PFTs that have been acclimated to hot weather.

In terms of the VPD, ORCHIDEE generally captures the negative dependence of Ω to VPD at a VPD smaller than 2 kPa. However, when the VPD is larger, observations show continuous negative dependence of Ω, while ORCHIDEE simulations show no significant changes in Ω with VPD. The decomposition into components of Ω shows that this difference is mainly contributed by different dependence of Gs on VPD (Fig. 6).

Compared with the observations, ORCHIDEE simulations show a different dependence of Ω to Rnet when the net radiation is <100 W m⁻². This difference is also mainly contributed by differences in Gs. In observations, the Gs SHAP values start to decrease rapidly when Rnet is lower than 200 W m⁻², while in ORCHIDEE simulations, the decrease of SHAP values is smaller and happens when Rnet is below 50 W m⁻².

Regarding the dependence of Ω to LAI, ORCHIDEE simulations show a significant increase of Ω with LAI across the entire range of LAI, due to a strong increase of Gs along with LAI, with the Gs SHAP values increasing by 0.2–0.4 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ from LAI=0 to LAI=5. However, the observations show that SHAP values increase only by less than 0.05 $\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the same change in LAI, resulting in a weak dependence of Ω on LAI.

Both observations and ORCHIDEE show weak dependence of Ω on canopy height. However, all of the data agree with a positive impact of canopy height on Ga. A strong increase of Ga is found when the height is below 15 m.

3.3 Interactions among factors

To further understand how the model biases the controls of Ω, we explored the interactions between factors that have significantly different impacts between ORCHIDEE and observations (Figs. 7 and 8).

Figure 7 The interaction between VPD and air temperature in controlling Ω (contribution of temperature) in forests (top) and in grasslands and croplands (bottom). The y axis is the SHAP value of Tair for Ω, colors indicate the VPD of each data point.

Figure 8 Same as Fig. 7 but for interactions between VPD and LAI.

The interactions between VPD and Tair are shown in Fig. 7. The observation data show that when the Ω SHAP value is positive ($\text{Tair}>\mathrm{25}{}^{\circ }$), data with larger VPD have smaller Ω values than those with smaller VPD.

In ORCHIDEE simulations, although Ω SHAP values vary differently along the temperature gradient compared with observations, similar interactions between VPD and Tair are also found, i.e., for a given temperature, when the Ω SHAP value is positive, large VPD values tend to result in smaller Ω. In another words, the dependence of Ω to Tair in hot weather is weakened by a high VPD level. This weakening of Ω dependence on Tair is due to weakened dependence of Gs on Tair under high VPD conditions (Fig. S3).

A similar interaction between VPD and LAI is also found in both the observations and ORCHIDEE simulations (Fig. 8). The data points with VPD>3 kPa show SHAP values close to zero, indicating that higher VPD tends to also weaken the dependence of Ω on LAI. ORCHIDEE underestimated the weakening effect of high VPD to the Ω to LAI dependence as the SHAP values under high VPD conditions remain very positive/negative compared with the observation.

4 Discussion

4.1 How can models correctly simulate the coupling strength

Accurately resolving the land–atmospheric water and energy exchanges is critical in simulating the climate system. To ensure this, LSMs must be carefully calibrated and validated with observations before use. The ORCHIDEE model has been calibrated several times for carbon and water fluxes against flux observations including the use of dedicated data-assimilation systems (e.g., Bastrikov et al., 2018). As a result, the ORCHIDEE model with the most recent set of parameters does not show large biases in LE (Fig. 3c).

Nevertheless, there remains no evaluation of the components and processes of LE, as well as their biotic and abiotic controls, leading to potential biases in LE simulation if climate changes. Disentangling and assessing processes and components of LE are difficult due to the lack of direct observation (Nelson et al., 2020). Although not perfect, evaluating the coupling strength and its components gives a possible way to further constrain the models.

In this study, we showed that the current ORCHIDEE model captures the coupling strength at most of the sites but fails to correctly represent the processes. The tuning of current LSMs often adjusts a few uncertain parameters to produce a small number of target variables (C fluxes, LE, sensible heat flux) close to the observation. In a complex model, this kind of calibration may result in overfitting and errors that compensate for each process. In the end, the model may get the correct result for the wrong reasons. Therefore, calibrating the model at the process level is helpful. For instance, the calibration of the a₁ and b₁ parameters in stomatal conductance calculation using independent observation-constrained values from Lin et al. (2015) leaf-scale data synthesis has significantly improved our estimation of f_vpd (Fig. S1), consequently correcting some biases in Gs and resulted in better Ω in short vegetation. In forest sites, Ω seems worse after this calibration, but this is because of the biases in modeled Ga, probably due to a bad assumption in calculating the displacement height.

In spite of the improvement from gs calibration, large biases in Gs remain in short vegetation (grasslands and croplands). Our analyses on the controlling factors shed light on where the problems are and give a direction to improve: we expect the model performance to improve if the dependence of Gs on temperature is corrected and the impact of VPD on stomatal conductance is further constrained. We did not do further calibration here because the responses of gs to VPD are an emergent area of concern for LSMs and more process-level modeling and calibration efforts remain needed (Yang et al., 2019). Also, it is out of the scope of this evaluation study. Nevertheless, the framework we used here would be helpful for models to identify their problematic processes and potentially fix their biases.

4.2 Factors controlling vegetation coupling strength

Due to the complexity of processes, as well as the lack of data, it is difficult to attribute the variation of coupling strength to different factors. Previous studies either focus on one or a few meteorological factors such as VPD, radiation or wind speed (Kumagai et al., 2004; Nicolás et al., 2008; Z. Z. Zhang et al., 2018), or biotic factors like LAI or PFT (Tateishi et al., 2010; Zhang et al., 2016). Our new framework for disentangling the impacts of different factors provides a systematic view to understand the impact of these factors.

Among all the factors, VPD was the most intensively investigated factor due to its strongest impact on stomatal conductance. A previous study showed that vegetation tends to be more decoupled in a wet season with low VPD compared to a dry season with high VPD (Kumagai et al., 2004). In this study, we found that VPD is the most important factor affecting Ω and affected Ω in similar way to the previous study (Fig. 6). This effect is mainly due to the reduction of Gs under dry conditions as plants tend to close the stomata under high VPD conditions to reduce water loss. In addition, high VPD conditions often coincide with low soil moisture, which hampers soil water uptake by plants, also leading to low Gs. It should be noted that this VPD–Ω relationship is obtained using daily data. At a subdaily timescale, this VPD–Ω relationship is not easily observed due to the strong impacts of other factors, such as radiation (Wullschleger et al., 2000; Z. Z. Zhang et al., 2018).

The impact of Tair on Ω is through two possible pathways. First, Tair can directly affect VPD by changing saturated water vapor pressure, leading to changes in Ω. Second, Tair can affect the photosynthesis rate by changing enzyme activities. Because stomatal conductance is strongly coupled with carbon assimilation rate (Cowan and Farquhar, 1977), the changes in photosynthesis rate can thus affect gs, and consequently Ω. In this study, we found that the responses of Ω and Gs to Tair differ from those to VPD, implying that the impacts of Tair through the second pathway is not negligible. The differential Tair impacts on Gs and Ω between observations and model simulations are probably due to an incorrect Tair adaptation of vegetation in the ORCHIDEE model.

Besides VPD and Tair, some studies found significant impacts from net radiation (Nicolás et al., 2008) or photosynthetically active radiation on Ω (which is strongly correlated to net radiation used in our analyses) (Z. Z. Zhang et al., 2018). Similar to Tair, changing radiation can also alter leaf photosynthesis rate. Due to the coupling between stomatal conductance and carbon assimilation, the changes in radiation thus result in Ω changes. Nevertheless, the impact of radiation should be considered with caution because radiation is strongly correlated with other environmental or biotic factors that have diurnal and seasonal cycles (e.g., temperature, LAI). Besides the short-term effect, long-term changes of radiation can affect soil moisture by altering LE, which may potentially change the coupling strength of the vegetation.

In terms of wind speed, we detected a negative dependence of Ω on wind as expected. This is because wind can accelerate the mixing of the boundary layer, increasing Ga. In this study, we did not find wind speed to be as important as VPD or vegetation types in explaining the variation of Ω. However, it needs to be kept in mind that the importance of factors depends on vegetation type. In ecosystems with a small vegetation cover (meaning small Gs), or in ecosystems where Gs has small variability, the importance of wind speed will increase.

Apart from the abiotic factors, the biotic factors or vegetation properties also play important roles in controlling Ω. The PFT is found to be the second most important factor affecting Ω after VPD in observation data (Fig. 4). In ORCHIDEE simulations, the PFT impact on Ω is weaker but still important, especially for different forest types. The pattern of Ω among PFTs found in this study agree well with De Kauwe et al. (2017). The influences from PFTs on Ω may be due to various reasons. Besides leaf area and canopy height (investigated in this study), different PFTs often have different canopy structure and leaf traits, leading to differences in Ga and Gs. Meanwhile, the climate and environmental conditions (e.g., soil types) that different PFTs adapted to are also different. More detailed data are needed to further explain the PFT impacts.

In the two biotic factors, canopy height is thought to be an important factor in affecting Ω because it directly affects the roughness length and the aerodynamic resistance (Ershadi et al., 2015). Higher canopies with larger roughness tend to enhance the turbulence for a given wind speed above the canopy. In this study, we found a positive but weak dependence of Ga on canopy height when the height is under 15 m. This result is consistent with Peng et al. (2019), who found that when controlling leaf area, Ω decreases (corresponding to Ga increase) with canopy height in vegetation with a height of <20 m. In higher canopies, Ga and Ω become less sensitive to canopy height.

Besides canopy height, LAI is also an important control. On the one hand, observations have shown that large LAI can increase the roughness (Alekseychik et al., 2017), which can lead to an increase of Ga. Along with LAI, leaf size might also be important in affecting the roughness and Ga, but is not available at most sites, neither simulated by ORCHIDEE model. On the other hand, LAI affects Gs since a larger LAI means a larger area for transpiration. This effect might be further regulated by environmental factors such as VPD (Fig. 8). Besides the influence from environmental factors, we also expect the impact of LAI on Gs to saturate for high LAI because of increasing self-shading. The shaded leaves in lower canopy tend to have smaller transpiration due to the low interception of radiation (Roberts et al., 1993), resulting in a decrease of average transpiration per leaf area. Also, the Gs at the ecosystem level is a synthesis of different processes including the vapor diffusion within the canopy. A large LAI may slow down the diffusion of water vapor within the canopy, potentially resulting in smaller Gs and smaller Ω.

4.3 Limitations

Although the simulations and analyses we performed in this study clearly showed how and why the ORCHIDEE LSM has biases in its estimation of the coupling strength, there are still some questions that need to be answered before we can calibrate the processes underlying these biases.

First, the coupling strength is the consequence of multiple processes. In this evaluation of Ω, strict criteria have been used to screen the data to have only time steps with LE mainly contributed by transpiration. The effect of other processes (e.g., soil evaporation) can potentially affect the coupling strength under some circumstances. For instance, the wetland Ω is also strongly affected by evaporation from open water. An understanding of these processes is also important, and our evaluation cannot draw conclusions on how well ORCHIDEE simulates these processes.

Second, due to the meteorological requirements of eddy covariance methods, the current selected observations have an incomplete coverage of the real meteorological conditions. We could not obtain valid observations under conditions with very low wind speeds. However, plants still transpire water to the atmosphere under such conditions. New observation methods are needed to fill this gap so that future calibrations can ensure the models to correctly simulate vegetation under the whole range of conditions.

The data used in this study are all daytime values. But for some vegetation types, transpiration also happens at nighttime (Dawson et al., 2007). Although the nighttime transpiration is smaller than the daytime transpiration, it can still affect the water and energy balance at longer timescales. These changes can potentially affect vegetation. However, the processes controlling the nighttime transpiration, as well as how coupled the ecosystems are at night remains poorly understood. Current LSMs also lack representations of such processes. We are not able to consider these processes in our evaluation/simulation.

Besides the missing processes, uncertainty may also come from the method to estimate Ω. In the observation-based estimates, Ga was estimated using an empirical method from Thom (1972), which was derived from a bean crop. The Ga estimates from this method are found to be 81 %–116 % of the estimates of a more physically based method (Knauer et al., 2017) in six forest sites. To test how biased Ga affects our evaluation, we increased/decreased Ga by 30 % and re-estimated Gs and Ω (Fig. S6 in the Supplement). We found that perturbing Ga does not result in large changes in Gs. When Ga is 30 % smaller than current observation-based estimates, we obtained smaller biases in Ga and Ω in the ORCHIDEE Ctrl simulation of forest PFTs. However, decreasing the reference Ga in short PFTs leads to even larger biases in Ω, indicating that the large biases in model vegetation coupling strength in short vegetation is not due to uncertainties in the observation-based estimates.

For Gs, the inverted Penman–Monteith equation may also result in some uncertainties. On the one hand, the energy budget is not always closed in flux observations. De Kauwe et al. (2017) used the value zero when soil heat flux observation is absent for estimating Gs, which could lead to biases in Gs and consequently Ω if the actual soil heat flux is not negligible. When the energy imbalance is corrected by adjusting the Bowen ratio following De Kauwe et al. (2017), we obtained larger Gs estimates (Fig. S6), resulting in even larger modeled Gs bias than in this study. The increased biases in the corrected Gs compensate for the existing biases in Ga, leading to a “good” performance of Ω simulation in forest PFTs. On the other hand, the Penman–Monteith equation is still not perfect for estimating LE. A recent study (McColl, 2020) showed that the linear approximation of the Clausius–Clapeyron relation in the Penman–Monteith equation can contribute ∼5.7 W m² biases to daytime and ∼1.2 W m² to nighttime LE. This bias is remarkable when there is a large difference between ambient air temperature and surface temperature (often with small Ga). A higher surface than ambient air temperature (daytime) tends to overestimate Gs in the inverted Penman–Monteith equation with observed LE, which can further overestimate Ω. However, since ORCHIDEE used the same method to estimate Gs as the observation, the uncertainties from the Penman–Monteith equation should not significantly affect our findings and conclusion.

5 Conclusion

In summary, in this study we evaluated the vegetation–atmosphere coupling strength, Ω, in the ORCHIDEE LSM using an observation-based dataset at 90 flux sites. We found that short vegetation (grassland and cropland) in ORCHIDEE is too tightly coupled to the atmosphere compared to the observation-based estimates, while the coupling strength of forests is generally well estimated by ORCHIDEE. Nevertheless, biases remain in both modeled Ga and Gs. Calibration of parameters controlling the dependence of the stomatal conductance to VPD reduces the biases of Gs in the ORCHIDEE model to a small extent and improves the Ω estimates in short vegetation. Using a set of random forest models and analyses on SHAP values, we found that vegetation tends to be more decoupled to the atmosphere at low wind speed, high temperature, low VPD and large LAI conditions and in short vegetation. ORCHIDEE generally agrees with this pattern but underestimated the VPD impacts when VPD is high, overestimated the contribution of LAI and did not correctly simulate the temperature dependence when temperature is high. Canopy height affects Ga but does not show a strong direct impact on Ω. Our results highlight the importance of observational constraints on simulating the vegetation–atmosphere coupling strength, which can help to improve the predictive accuracy of water fluxes in Earth system models.