Modeling the impacts of diffuse light fraction on photosynthesis in ORCHIDEE (v5453) land surface model

Aerosol- and cloud-induced changes in diffuse light have important impacts on the global land carbon cycle, as they alter light distribution and photosynthesis in vegetation canopies. However, this effect remains poorly represented or evaluated in current land surface models. Here, we add a light partitioning module and a new canopy light transmission module to the ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) land surface model (trunk version, v5453) and use the revised model, ORCHIDEE_DF, to estimate the fraction of diffuse light and its effect on gross primary production (GPP) in a multilayer canopy. We evaluate the new parameterizations using flux observations from 159 eddy covariance sites over the globe. Our results show that, compared with the original model, ORCHIDEE_DF improves the GPP simulation under sunny conditions and captures the observed higher photosynthesis under cloudier conditions in most plant functional types (PFTs). Our results also indicate that the larger GPP under cloudy conditions compared with sunny conditions is mainly driven by increased diffuse light in the morning and in the afternoon as well as by a decreased vapor pressure deficit (VPD) and decreased air temperature at midday. The observations show that the strongest positive effects of diffuse light on photosynthesis are found in the range from 5 to 20 ^∘C and at a VPD < 1 kPa. This effect is found to decrease when the VPD becomes too large or the temperature falls outside of the abovementioned range, which is likely due to the increasing stomatal resistance to leaf CO₂ uptake. ORCHIDEE_DF underestimates the diffuse light effect at low temperature in all PFTs and overestimates this effect at high temperature and at a high VPD in grasslands and croplands. The new model has the potential to better investigate the impact of large-scale aerosol changes and long-term changes in cloudiness on the terrestrial carbon budget, both in the historical period and in the context of future air quality policies and/or climate engineering.

1 Introduction

Process-based land surface models (LSMs), which simulate the water and energy balance as well as biogeochemical processes on land, have been widely used to attribute past changes in carbon (C) fluxes (Piao et al., 2009; Sitch et al., 2015) and to project the future land C budget (Ciais et al., 2013). Despite being useful and widely applied tools, large uncertainties are a limitation of LSMs (Sitch et al., 2008). One of the sources of the uncertainties is the omission or oversimplification of important processes that affect primary production. For instance, the impacts of light quality on photosynthesis is not currently represented in most LSMs, limiting the possibility to predict the variability of the carbon budget driven by changes in the atmospheric aerosol load which may be triggered by volcanic eruptions or variation in air pollution levels.

In situ observations have found that, under the same light level, the increase in the diffuse light fraction can enhance light use efficiency and ultimately photosynthesis or gross primary production (GPP; Gu et al., 2003; Niyogi et al., 2004; Misson et al., 2005; Alton, 2007a; Knohl and Baldocchi, 2008; Mercado et al. 2009; Oliphant et al., 2011; Kanniah et al., 2013; Williams et al., 2014; Cheng et al., 2015; Wang et al., 2018). Several mechanisms explaining this GPP enhancement have been proposed and tested. First, the more isotropic nature of diffuse light means that it penetrates deeper into the canopy to become available for photosynthesis of the lower canopy leaves, which would otherwise be shaded and light limited (Roderick et al., 2001; Urban et al., 2012). Second, the multi-directionality of diffuse light produces a more homogeneous distribution of radiation between sunlit and shaded leaves, enhancing the photosynthesis of upper canopy shaded leaves and limiting the waste of energy in light-saturated sunlit leaves (Li et al., 2014; Williams et al., 2014). Third, a higher diffuse light fraction is often accompanied by a less stressing temperature and vapor pressure deficit (VPD) for photosynthesis. The covariance of these environmental factors may also cause the GPP to increase under cloudier conditions, although not as a direct effect of diffuse light (Gu et al., 2002; Cheng et al., 2015; Li et al., 2014). Lastly, the plant LAI (leaf area index – the area of leaves per unit land area) maximum may become acclimated to the cloudier seasons, which also contributes to higher GPP (Williams et al., 2016).

Currently, most process-based LSMs simulate leaf photosynthesis using equations and parameterizations derived from Farquhar et al. (1980) with different formulations of stomatal conductance, usually with stomatal closure under a high VPD or low relative humidity (Ball et al., 1987; Yin et al., 2009; Medlyn et al., 2011). These parameterizations calculate photosynthesis per unit LAI considering the stress from temperature, VPD and soil water, and then integrate it over the entire canopy volume. Therefore, the effects of temperature and VPD change under cloudier conditions have usually been implicitly considered in current LSMs (e.g., Zhang et al., 2019). However, for the sake of simplicity and computational efficiency and due to the lack of diffuse light fraction data, most global LSMs have assumed a single extinction coefficient for both direct and diffuse light (Sellers et al., 1997; Sitch et al., 2008). Therefore, these LSMs are incapable of investigating the effect of the diffuse light fraction changes on photosynthesis. This limit of LSMs is thought to cause considerable underestimation of the land C sink after the eruption of Mount Pinatubo (Le Quéré et al., 2018).

To better simulate the diffuse light impacts, several earlier works have developed photosynthesis models that consider different light transmission of diffuse and direct radiation (Spitters, 1986; Leuning et al., 1995; de Pury and Farquhar, 1997). Based on these models, a few studies have tried to address the influence of light quality on GPP in LSMs. Dai et al. (2004) introduced a two-big-leaf canopy model to simulate the effects of diffuse and direct radiation in the Common Land Model (CLM 2L). This two-big-leaf scheme was further used in iTem LSM (Chen and Zhuang, 2014) and was partly inherited in later CLM models (Oleson et al., 2013). However, this light transmission scheme assumes a single-layer canopy and can, therefore, not simulate the vertical profile of leaf traits. A multilayer canopy model is more suitable to represent the vertical heterogeneity of leaf traits and radiation transfer (Alton et al., 2007b; Bonan et al., 2012). Differentiating sunlit and shaded leaves in a multilayer canopy LSM was firstly considered in the Joint UK Land Environment Simulator (JULES) LSM (Alton et al., 2007a; Mercado et al., 2009). Using this version of JULES, Mercado et al. (2009) investigated the diffuse light effect and suggested that diffuse light fraction change enhanced the global land C sink during the 1960–1999 period by about a quarter. However, the model of Mercado et al. (2009) was only tested at two forest sites, which cannot sufficiently represent global terrestrial ecosystems. Thus, the need remains to obtain well-evaluated LSMs that distinguish diffuse and direct light in order to test the results of Mercado et al. (2009) and to further investigate the diffuse radiation effect of aerosols. Apart from JULES, the Yale Interactive terrestrial Biosphere model (YIBs) also included a two-stream multilayer canopy light transmission scheme, but few efforts have been made to evaluate the ability of the YIBs model to capture the observed diffuse light fertilization effect, especially at sub-daily timescales (Yue and Unger, 2015).

Here, we introduce a modified version of the LSM ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems; Krinner et al., 2005), referred to as ORCHIDEE_DF, which uses a semiempirical method to calculate the fraction of diffuse light (Weiss and Norman, 1985) as well as a process-based multilayer canopy light transmission model to simulate the effects of diffuse light fraction on photosynthesis (Spitters, 1986). We evaluated the GPP simulated by ORCHIDEE_DF and the same version of the ORCHIDEE code without diffuse light (trunk version, v5453) using observations collected from 159 eddy covariance flux sites over 11 plant functional types (PFTs; Baldocchi et al., 2001). Using both model simulations and observations at the flux sites, we also investigated the interactions between diffuse light fraction and biotic and abiotic factors on GPP, with the objective of understanding when and how much the light quality affects photosynthesis. Because diffuse light is expected to enhance the photosynthesis of shaded leaves in deep canopy, we tested whether the enhancement of GPP due to diffuse radiation is larger in canopies with a larger LAI. We also tested whether environmental factors such as temperature or VPD affect this enhancement from diffuse light.

2 Data and methods

2.1 Model description

2.1.1 Canopy light transmission and photosynthesis in the ORCHIDEE trunk

The ORCHIDEE_DF model is based on ORCHIDEE trunk version 5453 (updated in September 2018). A general description of the physical processes related to the energy and water balance, vegetation dynamics and biogeochemical processes in ORCHIDEE can be found in Krinner et al. (2005). The ORCHIDEE trunk version 5453 brings a number of improvements, and photosynthesis parameters were recently recalibrated against FLUXNET data (Baldocchi et al., 2001) and atmospheric CO₂ observations for the IPSL Earth system model (IPSL-CM6) and the CMIP6 simulations.

The leaf-scale photosynthesis calculation in ORCHIDEE trunk is based on the scheme of Yin and Struik (2009). This scheme is an adaptation of the biophysical model of Farquhar et al. (1980) with a specific parameterization of stomatal conductance. The Farquhar et al. (1980) model calculates assimilation (A) as the minimum of the Rubisco-limited rate of CO₂ assimilation (A_c) and the electron-transport-limited rate of CO₂ assimilation (A_j):

$$\begin{array}{}\text{(1)}& A=\min \left\{A_{\mathrm{c}},A_j\right\}\end{array}$$

Here, A_c is mainly affected by the maximum carboxylation capacity of Rubisco (Vcmax), which is temperature dependent (Yin and Struik, 2009), and the CO₂ concentration at the carboxylation site (C_c):

$$\begin{array}{}\text{(2)}& A_{\mathrm{c}}={\displaystyle \frac{(C_{\mathrm{c}}-{\mathrm{\Gamma }}^{\ast })\mathrm{Vcmax}}{C_{\mathrm{c}}+K{m}_{\mathrm{C}}(\mathrm{1}+O/{\mathrm{Km}}_{\mathrm{O}})}}-R_{\mathrm{d}},\end{array}$$

where Γ^* is the CO₂ compensation point in the absence of dark respiration (R_d); Km_C and Km_O are the Michaelis–Menten constants for CO₂ and O₂, respectively; and O is the O₂ concentration at the carboxylation site.

A_j is calculated as a function of C_c and the electron transport rate (J):

$$\begin{array}{}\text{(3)}& A_j={\displaystyle \frac{J(C_{\mathrm{c}}-{\mathrm{\Gamma }}^{*})}{\mathrm{4.5}{C}_{\mathrm{c}}+\mathrm{10.5}{\mathrm{\Gamma }}^{*}}}-R_{\mathrm{d}}\end{array}$$

Here, J is determined by a temperature-dependent maximum electron transport rate (Jmax) and the photosynthetic photons absorbed by leaves, calculated following Yin and Struik (2009). Due to the attenuation of photosynthetically active radiation (PAR) with depth in the canopy, J also varies vertically. In addition, to account for the distribution of light and maximize the assimilation, plants tend to allocate nitrogen unevenly in the canopy profile (Niinemets et al., 1998; Meir et al., 2002), resulting in a vertical gradient in the enzyme concentration and, consequently, in Vcmax and Jmax. The vertical heterogeneity of canopy photosynthetic properties highlights the need to represent the canopy in a multilayer way.

In order to simulate the vertical transmission and absorption of light within the canopy, ORCHIDEE trunk uses a multilayer canopy with a big-leaf approximation in each layer. The canopy is geometrically divided into up to a maximum number of 20 layers depending on the leaf area index (LAI). The discretization is represented in Fig. 1a, and the LAI at the interface of the layers are given by

$$\begin{array}{}\text{(4)}& \mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i=\mathrm{12}\times {\displaystyle \frac{e^{\mathrm{0.15}\times (i-\mathrm{1})}-\mathrm{1}}{e^{\mathrm{0.15}\times \mathrm{20}}-\mathrm{1}}},\end{array}$$

where LAI_c_i is the cumulative LAI above layer i ($\mathrm{1}\le i\le \mathrm{20}$), and the layers are numbered from top to bottom (Fig. 3). It should be noted that 20 layers are only for canopies with a total LAI larger than 12. The number of layers decreases with the total LAI. For instance, if the LAI is 2, only the first 10 layers are used to calculate the light distribution and photosynthesis (Fig. 1a).

Figure 1 The distribution of light and leaves in the canopy: (a) light distribution in ORCHIDEE trunk; (b) distribution of sunlit and shaded leaves in canopy in ORCHIDEE_DF; and (c) light absorbed by shaded leaves in each canopy layer under different solar zenith angles (SZAs) and fractions of diffuse light (Fdf) in ORCHIDEE_DF. Panel (d) is the same as panel (c) but for sunlit leaves. I denotes downward PPFD at the top of the canopy; I_abs denotes PPFD absorption per leaf area in ORCHIDEE trunk; I_shd,abs denotes PPFD absorption per leaf area in shaded leaves; and I_sun,abs denotes PPFD absorption per leaf area in sunlit leaves.

Light transmission in the multilayer canopy is calculated using the Beer–Lambert law (Monsi and Saeki, 2005) without distinguishing direct and diffuse light. The downward shortwave radiation arriving at the top of canopy (TOC) layer i (I_i) is

$$\begin{array}{}\text{(5)}& I_i=I_{\mathrm{0}}{e}^{-k\times \mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i},\end{array}$$

where k is the light extinction coefficient, taken equal to 0.5; and I₀ is the TOC downward shortwave radiation (W m⁻²).

Because the radiation attenuation between one layer and the one just below is assumed to be due to leaf absorption, the absorbed radiation per leaf area at the top of layer i (Iabs_i) can be estimated as in Saeki (1960):

$$\begin{array}{}\text{(6)}& {\mathrm{Iabs}}_i={\displaystyle \frac{-\mathrm{d}I}{\mathrm{dLAI}\mathrm{\_}\mathrm{c}}}\mathrm{|}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i=k{I}_{\mathrm{0}}{e}^{-k\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i}\end{array}$$

Here, we assume that all canopy layers are thin enough to neglect the difference in light absorption within each canopy layer, i.e., the absorbed radiation does not attenuate within each canopy layer and Iabs_i is used for all leaves in layer i.

It should be noted that the radiation considered to calculate the J term in Eq. (3) is not shortwave radiation in watts per square meter (W m⁻²) but photosynthetic photon flux density (PPFD) in micromoles per square meter per second (µmol m⁻² s⁻¹). Thus, a translation from Iabs_i to the absorbed PPFD per leaf area in canopy layer i (PPFDabs_i) is needed. Currently, there is no standard definition of the wavelength range for shortwave radiation (e.g., Howell et al., 1983; Zhang et al., 2004; Chen et al., 2012). In ORCHIDEE trunk, shortwave radiation (in W m⁻²) is multiplied by a factor of 0.5 to calculate photosynthetically active radiation (PAR; in W m⁻²), and then a quanta-to-energy ratio of 4.6 mmol J⁻¹ is used to convert PAR into PPFD (in µmol m⁻² s⁻¹).

ORCHIDEE accounts for a vertical gradient in the enzyme concentration in the canopy. Vcmax and Jmax are assumed to be linearly related to photosynthetically active leaf nitrogen concentration (per leaf area) in the model (Kattge and Knorr, 2007). Meir et al. (2002) found a decreasing leaf nitrogen concentration as well as decreasing Vcmax and Jmax with increasing canopy depth in different ecosystems, suggesting an acclimation of plants to maximize photosynthesis in a canopy with unevenly distributed radiation. ORCHIDEE trunk lacks an explicit model of dynamic nitrogen allocation to leaves in the canopy; instead, it uses an empirical relationship to represent the impact of the leaf nitrogen concentration on Vcmax and Jmax using the vertical profile of radiation:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathrm{Vcmax}}_i={\mathrm{Vcmax}}_{\mathrm{0}}\left(\mathrm{1}-\mathrm{0.7}\times \left(\mathrm{1}-{\displaystyle \frac{I_i}{I_{\mathrm{0}}}}\right)\right)\text{(8)}& {\displaystyle }{\mathrm{Jmax}}_i={\mathrm{Jmax}}_{\mathrm{0}}\left(\mathrm{1}-\mathrm{0.7}\times \left(\mathrm{1}-{\displaystyle \frac{I_i}{I_{\mathrm{0}}}}\right)\right)\end{array}$$

It should be noted that the leaf-scale assimilation variables (e.g., Vcmax) in ORCHIDEE trunk are also affected by the instantaneous leaf temperature and the temperature of the last that which plants have adapted to (Kattge and Knorr, 2007). As there is only one energy budget per grid cell, from which we cannot determine the leaf temperature, in the current ORCHIDEE model, the air temperature is used to represent the leaf temperature in current model. The calculation of C_c depends on the VPD and also on whether the vegetation follows the C₃ or C₄ photosynthesis pathway (Yin and Struik, 2009). For simplicity, the near-surface air temperature and humidity are used for the calculation of assimilation in all canopy layers. Furthermore, there are 13 PFTs in ORCHIDEE (Table S1), and Vcmax and Jmax are PFT dependent.

2.1.2 Canopy light transmission in ORCHIDEE_DF

In ORCHIDEE_DF, we use the same stratification of canopy as in the trunk version (Eq. 4). However, for the light transmission, we use a two-stream radiative transfer model with direct and diffuse radiation treated separately following Spitters (1986). For convenience, we use radiation and I in this section to refer to the PPFD derived from the light partitioning step (see Sect. 2.1.3).

An assumption of the model is that leaves are bi-Lambertian surfaces for radiation, i.e., the reflection and transmission are isotropic. This reflection and transmission are together referred to as leaf scattering. This assumption implies that once direct radiation encounters a leaf, it is either absorbed or scattered as diffuse light. For diffuse radiation, however, the scattered light remains diffuse. The scattering coefficient, σ, is assumed to be equal to 0.2 following Spitters (1986) – meaning that 20 % of the light encountering a leaf is scattered (80 % is absorbed).

Based on this assumption, the radiation penetrating the canopy can be divided into three components (Fig. 3): the direct light which has not been intercepted by leaves (I_dr,dr), the diffuse light generated by leaf scattering of intercepted direct light (I_dr,df) and the diffuse light in the canopy provided by the TOC diffuse radiation (I_df). It should be noted that the diffuse light generated after direct light has been scattered multiple times is grouped into I_dr,df, whereas light from the scattering of TOC diffuse radiation belong to I_df (Fig. 3). The sum of I_dr,dr and I_dr,df, hereafter noted as I_dr, represents the total radiation in each canopy layer derived from the TOC direct radiation, hereafter I_dr,0.

If we also consider direct radiation as parallel beams, only the first leaves in the way of direct light can absorb I_dr,dr. These leaves are referred to as sunlit leaves. The fraction of sunlit leaves in each canopy layer can be calculated by applying the Beer–Lambert law using an extinction coefficient for opaque, nonreflective “black” leaves (Fig. 1b):

$$\begin{array}{}\text{(9)}& {\mathrm{LAIf}}_{\mathrm{sun},i}=e^{-k_{\mathrm{b}}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i}\end{array}$$

Here, LAIf_sun,i is the fraction of sunlit LAI in canopy layer i; LAI_c_i is the cumulative LAI in Eq. (4); k_b is the extinction coefficient if the leaves are assumed to be “black”. A function of θ, the leaf angle distribution index (LA) and leaf clumping index (LC) is used to represent the geometry between the direct radiation and leaves:

$$\begin{array}{}\text{(10)}& k_{\mathrm{b}}={\displaystyle \frac{\mathrm{LA}\cdot \mathrm{LC}}{\cos \mathit{\theta }}}\end{array}$$

For spherically distributed leaves, the LA equals 0.5 (Goudriaan, 1977; Bodin and Franklin, 2012). The LC is defined as in Myneni et al. (1989) and Baldocchi and Wilson (2001) – varying between 0 and 1. Here, we use the value 0.85 instead of 0.84, as recommended by an observationally based study (Baldocchi and Wilson, 2001).

The leaves which cannot absorb I_dr,dr are referred to as shaded leaves. Thus, the fraction of the shaded LAI in canopy layer i (LAIf_shd,i) is the complement of LAIf_sun,i:

$$\begin{array}{}\text{(11)}& {\mathrm{LAIf}}_{\mathrm{shd},i}=\mathrm{1}-{\mathrm{LAIf}}_{\mathrm{sun},i}\end{array}$$

Because I_dr,dr is assumed not to be transmitted as direct radiation through leaves, $I_{\mathrm{dr},\mathrm{dr},i}$, which represents I_dr,dr at layer i, can be calculated in a similar manner as in Eq. (9) using the downward direct radiation at the top of the canopy (I_dr,0):

$$\begin{array}{}\text{(12)}& I_{\mathrm{dr},\mathrm{dr},i}=I_{\mathrm{dr},\mathrm{0}}{e}^{-k_{\mathrm{b}}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i}\end{array}$$

The transmission of I_dr,df is difficult to estimate directly. Here, we calculate it as the difference between I_dr and I_dr,dr in each layer:

$$\begin{array}{}\text{(13)}& I_{\mathrm{dr},\mathrm{df},i}=I_{\mathrm{dr},i}-I_{\mathrm{dr},\mathrm{dr},i},\end{array}$$

where $I_{\mathrm{dr},\mathrm{df},i}$ and I_dr,i represent net (downward minus upward) I_dr,df and net I_dr at layer i, respectively.

The calculation of I_dr,i is based on Goudriaan (1982) and Hikosaka et al. (2016) under the assumptions that there is no difference in optical traits between leaves from different canopy layers and that the canopy is deep enough to neglect the reflection of the soil:

$$\begin{array}{}\text{(14)}& I_{\mathrm{dr},i}=(\mathrm{1}-\mathit{\rho })I_{\mathrm{dr},\mathrm{0}}{e}^{-\sqrt{\mathrm{1}-\mathit{\sigma }}{k}_{\mathrm{b}}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i},\end{array}$$

where ρ indicates the canopy reflection coefficient (i.e., the ratio between the TOC downward and upward radiation), calculated as

$$\begin{array}{}\text{(15)}& \mathit{\rho }=\left({\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-\mathit{\sigma }}}{\mathrm{1}+\sqrt{\mathrm{1}-\mathit{\sigma }}}}\right)\left({\displaystyle \frac{\mathrm{2}}{\mathrm{1}+\mathrm{1.6}\cos \mathit{\theta }}}\right)\end{array}$$

In contrast to the direct light transmission, the diffuse light will not change its directional characteristics when scattered by leaves. Similar to Eq. (5), net I_df at canopy layer i (I_df,i) can be estimated using TOC downward diffuse radiation (I_df,0) in a Beer–Lambert equation:

$$\begin{array}{}\text{(16)}& I_{\mathrm{df},i}=(\mathrm{1}-\mathit{\rho })I_{\mathrm{df},\mathrm{0}}{e}^{-k_{\mathrm{d}}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i},\end{array}$$

where k_d is the light extinction coefficient for diffuse light, calculated following Spitters (1986) as

$$\begin{array}{}\text{(17)}& k_{\mathrm{d}}=\mathrm{0.8}\sqrt{\mathrm{1}-\mathit{\sigma }}\end{array}$$

Similar to Eq. (6), the flux of light that is absorbed per canopy leaf area in layer i from I_df (Iabs_df,i), I_dr (Iabs_dr,i) and I_dr,dr (${\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i}$), respectively, can be written as

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathrm{Iabs}}_{\mathrm{df},i}={\displaystyle \frac{-\mathrm{d}{I}_{\mathrm{df}}}{\mathrm{dLAI}\mathit{\_}\mathrm{c}}}\mathrm{|}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i=k_{\mathrm{d}}{I}_{\mathrm{df},i}\text{(19)}& {\displaystyle }{\mathrm{Iabs}}_{\mathrm{dr},i}={\displaystyle \frac{-\mathrm{d}{I}_{\mathrm{dr}}}{dLAI\mathit{\_}c}}\mathrm{|}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i=\sqrt{\mathrm{1}-\mathit{\sigma }}{k}_{\mathrm{b}}{I}_{\mathrm{dr},i}\text{(20)}& {\displaystyle }{\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i}={\displaystyle \frac{-\mathrm{d}{I}_{\mathrm{dr},\mathrm{dr}}}{\mathrm{dLAI}\mathit{\_}\mathrm{c}}}\mathrm{|}\mathrm{LAI}\mathrm{\_}{\mathrm{c}}_i=k_{\mathrm{b}}{I}_{\mathrm{dr},\mathrm{dr},i}\end{array}$$

The I_dr,df absorbed per canopy leaf area by layer i (${\mathrm{Iabs}}_{\mathrm{dr},\mathrm{df},i})$ is

$$\begin{array}{}\text{(21)}& {\mathrm{Iabs}}_{\mathrm{dr},\mathrm{df},i}={\mathrm{Iabs}}_{\mathrm{d}r,i}-{\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i}\end{array}$$

It should be noted that all leaves can absorb diffuse radiation. Therefore, Eqs. (18) and (21) also represent the absorption of I_df and I_dr,df at the leaf scale. However, I_dr,dr is only absorbed by sunlit leaves; thus, the absorption of I_dr,dr per sunlit leaf area does not equal ${\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i}$, which is the average at the canopy scale. Instead, because I_dr,dr does not change its intensity, the absorption of I_dr,dr per sunlit leaf area can be written as

$$\begin{array}{}\text{(22)}& {\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i,\mathrm{sun}}=(\mathrm{1}-\mathit{\sigma })k_{\mathrm{b}}{I}_{\mathrm{dr},\mathrm{0}}\end{array}$$

We have assumed that shaded leaves can only absorb diffuse light. Thus, the radiation absorbed (per leaf area) by the shaded leaves layer i ($I_{\mathrm{shd},i,\mathrm{abs}})$ is

$$\begin{array}{}\text{(23)}& {\mathrm{Iabs}}_{\mathrm{shd},i}={\mathrm{Iabs}}_{\mathrm{df},i}+{\mathrm{Iabs}}_{\mathrm{dr},\mathrm{df},i}\end{array}$$

Sunlit leaves absorb both direct light and diffuse light; thus, the radiation received by sunlit leaves can be calculated as

$$\begin{array}{}\text{(24)}& {\mathrm{Iabs}}_{\mathrm{sun},i}={\mathrm{Iabs}}_{\mathrm{shd},i}+{\mathrm{Iabs}}_{\mathrm{dr},\mathrm{dr},i,\mathrm{sun}}\end{array}$$

Apart from light transmission, all other parameters (e.g., Vcmax, Jmax) in ORCHIDEE_DF are kept the same as in ORCHIDEE trunk.

2.1.3 Light partitioning in ORCHIDEE_DF

The lack of light quality (diffuse light fraction) information in most datasets to drive LSMs is one of the main difficulties when simulating the diffuse light effect. This field can be calculated in atmospheric light transmission models when aerosol and cloud information is available (Yue and Unger, 2017; Malavelle et al., 2019). However, the aerosol and cloud information is not always available. Here, we use the empirical equations following Weiss and Norman (1985) to partition the 30 min downward PAR, which can be derived from the shortwave radiation, into diffuse and direct components. Compared with another empirical method (Spitters et al., 1986), we found that this method more successfully reproduces the observed diffuse light fraction at the flux sites used in this study (Figs. 2, S1). The diffuse PAR fraction (Fdf_PAR) above the canopy is estimated as follows:

$$\begin{array}{}\text{(25)}& {\mathrm{Fdf}}_{\mathrm{PAR}}=\mathrm{1}-{\displaystyle \frac{{\mathrm{PAR}}_{\mathrm{p},\mathrm{dr}}}{{\mathrm{PAR}}_{\mathrm{p}}}}(\mathrm{1}-{\left({\displaystyle \frac{a-R}{b}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}),\end{array}$$

where PAR_p and PAR_p,dr are the potential total and direct PAR, i.e., the total and direct PAR which would arrive at land surface under clear-sky conditions. a and b are parameters, which take values of 0.9 and 0.7, and R is the ratio of observed to potential total downward shortwave radiation (SW_obs and SW_p) reaching the top of the canopy:

$$\begin{array}{}\text{(26)}& R={\displaystyle \frac{{\mathrm{SW}}_{\mathrm{obs}}}{{\mathrm{SW}}_{\mathrm{p}}}}\end{array}$$

The potential downward shortwave radiation consists of potential downward PAR (visible, range 0.4–0.7 µm) and potential downward near-infrared radiation (NIR, range 0.7–5 µm). Also, the potential PAR and NIR are the sum of direct (PAR_p,dr, NIR_p,dr) and diffuse (PAR_p,df, NIR_p,df) components, given by

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\mathrm{SW}}_{\mathrm{p}}& ={\mathrm{PAR}}_{\mathrm{p}}+{\mathrm{NIR}}_{\mathrm{p}}={\mathrm{PAR}}_{\mathrm{p},\mathrm{dr}}+{\mathrm{PAR}}_{\mathrm{p},\mathrm{df}}+{\mathrm{NIR}}_{\mathrm{p},\mathrm{dr}}\\ & +{\mathrm{NIR}}_{\mathrm{p},\mathrm{df}}\end{array}\end{array}$$

A simple atmospheric light transfer model modified from Weiss and Norman (1985) is used to estimate potential radiation. The potential direct PAR, PAR_p,dr is calculated as follows:

$$\begin{array}{}\text{(28)}& {\mathrm{PAR}}_{\mathrm{p},\mathrm{dr}}={\mathrm{PAR}}_{\mathrm{TOA}}{e}^{-\mathrm{0.185}(p/p_{\mathrm{0}})m}\cos \mathit{\theta },\end{array}$$

where PAR_TOA is the PAR at top of atmosphere (TOA); p and p₀ indicate the local and standard sea level air pressure, respectively; and m is the optical air mass, calculated using the solar zenith angle θ:

$$\begin{array}{}\text{(29)}& m=(\cos \mathit{\theta }{)}^{-\mathrm{1}}\end{array}$$

The potential diffuse TOC PAR, PAR_p,df is assessed as

$$\begin{array}{}\text{(30)}& {\mathrm{PAR}}_{\mathrm{p},\mathrm{df}}=\mathrm{0.4}({\mathrm{PAR}}_{\mathrm{TOA}}\cos \mathit{\theta }-{\mathrm{PAR}}_{\mathrm{p},\mathrm{dr}})\end{array}$$

and expresses that 40 % of the PAR flux that is extinguished in the atmosphere through scattering and absorption is available as diffuse PAR at the surface. Similarly, the potential direct and diffuse NIR at the top of the canopy (NIR_p,dr and NIR_p,df respectively), can be estimated as follows:

$$\begin{array}{}\text{(31)}& {\displaystyle }{\mathrm{NIR}}_{\mathrm{p},\mathrm{dr}}=({\mathrm{NIR}}_{\mathrm{TOA}}{e}^{-\mathrm{0.06}(p/p_{\mathrm{0}})m}-\mathit{\omega })\cos \mathit{\theta }\text{(32)}& {\displaystyle }{\mathrm{NIR}}_{\mathrm{p},\mathrm{df}}=\mathrm{0.6}({\mathrm{NIR}}_{\mathrm{TOA}}\cos \mathit{\theta }-{\mathrm{NIR}}_{\mathrm{p},\mathrm{dr}}-\mathit{\omega }\cos \mathit{\theta }),\end{array}$$

where ω is a flux term accounting for atmospheric water vapor absorption, calculated as a function of the solar constant (SC, in W m⁻²) and m:

$$\begin{array}{}\text{(33)}& \mathit{\omega }=\mathrm{SC}\times {\mathrm{10}}^{(-\mathrm{1.195}+\mathrm{0.4459}{\log }_{\mathrm{10}}m-\mathrm{0.0345}({\log }_{\mathrm{10}}m{)}^{\mathrm{2}})}\end{array}$$

Using the results from Eqs. (28), (30), (31) and (32), we are able to calculate the SW_p to obtain the value of R in Eq. (26).

It should be noted that the quanta-to-energy ratio (in mmol J⁻¹) is different under different sky conditions, because atmospheric scattering varies spectrally with the air mass and the cloud amount (Dye, 2004). For this consideration, the calculation of PPFD from PAR in ORCHIDEE_DF uses the observation-oriented empirical equations from Dye (2004):

$$\begin{array}{}\text{(34)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{t}}=\mathrm{4.576}-\mathrm{0.03314}{\mathrm{Fdf}}_{\mathrm{PAR}}\text{(35)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{df}}={\displaystyle \frac{\mathrm{4.5886}{\mathrm{Fdf}}_{\mathrm{PAR}}}{\mathrm{0.010773}+{\mathrm{Fdf}}_{\mathrm{PAR}}}},\end{array}$$

where the β_t is the quanta-to-energy ratio for the total PAR (PAR_t) at the top of the canopy, whereas β_df is for its diffuse component (PAR_df):

$$\begin{array}{}\text{(36)}& {\displaystyle }{\mathrm{PPFD}}_{\mathrm{t}}={\mathit{\beta }}_{\mathrm{t}}{\mathrm{PAR}}_{\mathrm{t}}\text{(37)}& {\displaystyle }{\mathrm{PPFD}}_{\mathrm{df}}={\mathit{\beta }}_{\mathrm{df}}{\mathrm{PAR}}_{\mathrm{df}}\end{array}$$

The diffuse PPFD fraction (Fdf_PPFD) can thus be calculated as follows:

$$\begin{array}{}\text{(38)}& {\mathrm{Fdf}}_{\mathrm{PPFD}}={\displaystyle \frac{{\mathrm{PPFD}}_{\mathrm{df}}}{{\mathrm{PPFD}}_{\mathrm{t}}}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{df}}}{{\mathit{\beta }}_{\mathrm{t}}}}{\mathrm{Fdf}}_{\mathrm{PAR}}\end{array}$$

2.2 Flux data and site-level simulations

To evaluate ORCHIDEE_DF, we collected flux site measurements from the La Thuile dataset, which includes 965 site-year observations from 252 sites in total (https://fluxnet.fluxdata.org/data/la-thuile-dataset/, last access 1 November 2020). Because our ORCHIDEE simulations assume that the ecosystems are in equilibrium and do not experience disturbances (e.g., logging, fire), we selected flux sites without strong disturbances over the last 10 years. For sites that also provided growing season LAI information, we also removed forests site with LAI < 2, which may be considered as sparse forests with understory vegetation. In the end, observations of 655 site-years from 159 sites were retained (Table S2). The annual mean temperature of the sites spans from −9 to 27 ^∘C, and the annual precipitation spans from 67 mm yr⁻¹ to over 3000 mm yr⁻¹ (Fig. S2), which is representative of most of the climate conditions over the globe. The dataset provides in situ meteorology, net ecosystem exchange (NEE), gross primary productivity (GPP) and data quality information at 30 min time steps. The GPP provided by this dataset is partitioned from NEE and gap filled using the method of Reichstein et al. (2005). Specifically, 64 of the 159 sites provided measurements of both total and diffuse PPFD, which allows us to evaluate the light partitioning parametrization (Eqs. 25–38). The gaps and missing variables in meteorology are filled using the approach from Vuichard and Papale (2015) to meet the model input requirements.

Because ORCHIDEE has different photosynthesis parameters for different PFTs, we classified the vegetation at each site into the 13 ORCHIDEE PFTs (Table S1) according to the International Geosphere-Biosphere Programme (IGBP) land cover types specified on the FLUXNET website (https://fluxnet.org/, last access: 1 November 2020). If the IGBP land cover type is not specified or may match more than one ORCHIDEE PFT (e.g., shrublands, savannas and wetlands), the PFT is determined according to the dominant plant species described in related references. Specifically, the mixed forests (MF) type exists in the IGBP classification but not in the ORCHIDEE PFTs. Because MF sites are mostly located in temperate regions, we assume that they are composed of 50 % temperate broadleaf deciduous forests and 50 % temperature evergreen needleleaf forests. Detailed information of flux sites is found in Table S2.

To evaluate the model, 30-years spin-up simulations are firstly conducted on ORCHIDEE_DF at each site to equilibrate the leaf area index with site conditions. The simulations with 30 min output are then conducted with ORCHIDEE trunk and ORCHIDEE_DF using the full span of the FLUXNET La Thuile series at each site, respectively. It should be noted that we use the same spin-up for ORCHIDEE trunk and ORCHIDEE_DF to ensure the same initial states for the two simulations. A test has shown that different spin-up simulations do not affect the simulation of GPP in the following years (not shown).

2.3 Analyses

When evaluating the modeled GPP response to diffuse light, we did not use all the 30 min data points due to several concerns. First, all nighttime data points were excluded from the analyses given that GPP is zero at night. Second, all data points flagged with poor quality in the FLUXNET archive were removed. Third, ORCHIDEE might not be perfect at capturing the seasonality of leaf flushing and shedding. In order to minimize the uncertainty from phenology, we used only data from the growing season at each site, which is defined as months when

$$\begin{array}{}\text{(39)}& {\mathrm{GPP}}_{\mathrm{m}}>{\mathrm{GPP}}_{\mathrm{m},\min }+{\displaystyle \frac{{\mathrm{GPP}}_{\mathrm{m},\max }-{\mathrm{GPP}}_{\mathrm{m},\min }}{\mathrm{4}}},\end{array}$$

where GPP_mis the observed monthly GPP, and GPP_m,min and GPP_m,max are the observed minimum and maximum monthly GPP at the corresponding sites, respectively.

To assess the effect of variable diffuse light fraction on both GPP and light use efficiency (LUE – the ratio of GPP to incoming shortwave radiation), we look at the difference in GPP and the LUE under sunny and cloudy conditions. We define sunny and cloudy conditions as those when the fraction of diffuse PPFD at the top of the canopy (Fdf_PPFD) is smaller than 0.4 and greater than 0.8, respectively, and calculate the average sunny and cloudy GPP and LUE at each site. To ensure that the comparison between sunny and cloudy conditions are at the same PPFD level, the sunny time steps with PPFD larger than the maximum PPFD under cloudy conditions are removed from the average, and vice versa. In addition, to make sure that the difference in GPP between sunny and cloudy conditions is not an artifact of different LAI, sites with average modeled LAI under cloudy and sunny conditions that differ by more than 0.3 are excluded from this analysis.

3 Results

3.1 Diffuse light modeling

Figure 2 shows the relationship between 30 min modeled and measured Fdf_PPFD at flux sites (64 sites). The data points are generally distributed along the 1:1 line, indicating an unbiased estimation of our diffuse light model. In total, our simple model explains over 51 % of the variance in observed diffuse PPFD fraction. Although this model is imperfect, we currently have no better way of reproducing the diffuse PPFD at the flux site scale.

Figure 2 Modeled and observed diffuse PPFD fraction (a): the dark area in the scatter plot indicates high data density, whereas the light area indicates low data density. (b) Density distribution of the observed diffuse PPFD fraction. (c) Density distribution of the modeled diffuse PPFD fraction.

Figure 3 Diagram of the canopy light transmission in ORCHIDEE_DF. I_dr,0 denotes the downward direct PPFD at the top of the canopy; I_df,0 denotes the downward diffuse PPFD at the top of the canopy; LAI_c_i denotes the cumulative LAI above canopy layer i; $I_{\mathrm{dr},\mathrm{dr},i}$ denotes downward direct PPFD at the top of canopy layer i; $I_{\mathrm{dr},\mathrm{df},i}$ denotes the net diffuse PPFD derived from the scattering of I_dr,0 at the top of canopy layer i, which is equal to the difference of its downward ($I_{\mathrm{dr},\mathrm{df},\mathrm{down},i})$ and upward ($I_{\mathrm{dr},\mathrm{df},\mathrm{up},i})$ components; I_df,i denotes the net diffuse PPFD derived from I_df,0 at the top of canopy layer i, which is equal to the difference of its downward ($I_{\mathrm{df},\mathrm{down},i}$) and upward ($I_{\mathrm{df},\mathrm{up},i}$) components.

3.2 General model performance

The performance of both ORCHIDEE trunk and ORCHIDEE_DF for the 30 min GPP from each PFT (all sites) is presented in Fig. 4. Generally, ORCHIDEE trunk underestimated the standard deviation (SD) of GPP at a 30 min time step compared with observations as well as across all PFTs except boreal evergreen needleleaf forests and C₄ croplands (Fig. 4a). The correlation coefficients between ORCHIDEE trunk GPP and observations are generally between 0.5 and 0.7 among PFTs (Fig. 4b). In tropical broadleaf forests, this correlation coefficient is about 0.2, which is much smaller than in other PFTs and likely due to the limited seasonality of primary production in the tropics. The GPP simulated by ORCHIDEE_DF shows comparable performance to ORCHIDEE trunk, although with a slightly smaller SD (Fig. 4a).

Figure 4 Performance of ORCHIDEE trunk and ORCHIDEE_DF at different PFTs. (a) Taylor plot of GPP, where all valid 30 min observations are used as a reference: the filled circles indicate ORCHIDEE trunk, and open circles indicate ORCHIDEE_DF. (b) Comparison of the correlation coefficients between the two models against observations. Panels (c) and (e) are the same as panel (a) but for cloudy (diffuse light fraction > 0.8) and sunny (diffuse light fraction < 0.4) conditions only. Panels (d) and (f) are the same as panel (b) but for cloudy and sunny conditions only.

Similar evaluations on the GPP from the two models are performed under cloudy and sunny conditions, respectively (Fig. 4c, d, e, f). Under cloudy conditions, ORCHIDEE trunk and ORCHIDEE_DF both underestimated the GPP SD. The correlation coefficients to observations are generally between 0.5 and 0.8 (Fig. 4d). Compared with ORCHIDEE trunk, ORCHIDEE_DF shows slightly worse correlation coefficients but an improved SD for most of the PFTs except tropical broad-leaved evergreen forests (TrEBF) and temperate needleleaf evergreen forests (TeENF) (Fig. 4c).

Compared with cloudy conditions, the GPP simulated by the two models under sunny conditions shows a weaker correlation with observations. The correlation coefficients generally vary between 0.3 and 0.6 among PFTs. However, it should be noted that ORCHIDEE_DF more successfully reproduced GPP variation under sunny conditions compared with ORCHIDEE trunk in most PFTs except temperate deciduous broad-leaved forests (TeDBF) and C₄ croplands (C4Cro) (Fig. 4f). The GPP SD derived from ORCHIDEE trunk simulations under sunny conditions show larger variability among PFTs than under cloudy conditions, whereas the GPP SD under sunny and cloudy conditions show similar bias compared with observations for ORCHIDEE_DF (Fig. 4e).

3.3 Effects of diffuse light on GPP and LUE

Because the modification of ORCHIDEE_DF was limited to light transmission, the pertinent process-oriented evaluation of the two models should focus on their ability to capture the observed GPP differences between cloudy and sunny conditions (hereafter ΔGPP), rather than on correlations or RMSE with observations, that may result from different structural and parametric errors of the model, which are not related to diffuse light.

Figure 5 Observed GPP and GPP modeled by ORCHIDEE trunk and ORCHIDEE_DF under cloudy (diffuse light fraction > 0.8) and sunny (diffuse light fraction <0.4) conditions at selected sites (with relatively long time series) from each PFT.

Figure 5 shows the observed and modeled GPP under sunny and cloudy conditions at different PPFD levels at flux sites with relatively long time series of observations from each PFT. For all the sites selected, the observed GPP under cloudy conditions is larger than under sunny conditions. However, the GPP simulated by ORCHIDEE trunk shows no or a small difference between cloudy and sunny conditions at most sites. In contrast, ORCHIDEE_DF reproduces this GPP difference in most PFTs except tropical deciduous broad-leaved forests (TrDBF), Boreal deciduous broad-leaved forests (BoDBF) and C₄ grasslands (C4Gra). However, there is only one TrDBF site and very few C4Gra sites in our dataset. Furthermore, at most C4Gra sites, we are not able to find PPFD levels where sunny and cloudy conditions coexist. Therefore, we are not able to carry out further evaluation of cloudy-minus-sunny GPP differences for TrDBF and C4Gra. At three of the four BoDBF sites, the modeled GPP difference under cloudy and sunny conditions is relatively small (not shown). This might be because the model overestimated the deleterious effect of low temperature on photosynthesis at the BoDBF sites (mean annual temperature < 3 ^∘C). In total, observations from about 70 % of the sites show a remarkably higher GPP under cloudy conditions compared with sunny conditions. This percentage is only 30 % in ORCHIDEE trunk simulations but amounts to 60 % in ORCHIDEE_DF simulations.

To summarize the site-level results, we investigated the distribution of the GPP difference between cloudy and sunny conditions (hereafter refer to as ΔGPP; Fig. 6a). Observations and ORCHIDEE_DF show a positive bias in ΔGPP, with ΔGPP values between 0 and $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ gC m⁻² s⁻¹ at most sites. However, for ORCHIDEE trunk, ΔGPP is near zero at most sites. This result confirms that ORCHIDEE_DF performs much better than ORCHIDEE trunk in simulating differences in GPP under different light conditions.

Figure 6 Site distribution of (a) the GPP difference between cloudy (diffuse light fraction >0.8) and sunny (diffuse light fraction < 0.4) conditions. Panel (b) is the same as panel (a) but for the LUE. The sunny and cloudy time steps are sampled at equal light levels.

It should be noted that ΔGPP can be affected by PPFD. At sites where sunny and cloudy conditions only coexist at a relatively low PPFD level, the ΔGPP should also be small. To remove the effect of the PPFD level on ΔGPP, we analyzed the difference in the LUE, i.e., ΔLUE, between the two conditions (Fig. 6b). Compared with ΔGPP, positive ΔLUE values are more evenly distributed around 0–$\mathrm{15}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol (of photon)⁻¹ for observations and the ORCHIDEE_DF simulation. For ORCHIDEE trunk, the ΔLUE has a range of 0–$\mathrm{8}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹, with the upper range smaller than in the observations and ORCHIDEE_DF.

Figure 7 Same as in Fig. 6 but differentiated for three times of the day.

We further refined this analysis to investigate if the effects of diffuse light differ at different times of the day (Fig. 7). Results for three different periods of the day show that cloudy conditions result in a higher GPP of 0–$\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ gC m⁻² s⁻¹ than sunny conditions at most sites in the morning and afternoon, which is generally captured by ORCHIDEE_DF but missed by ORCHIDEE trunk in the morning (Fig. 6a, c). At midday, due to the dependence of Fdf on the PPFD (Eqs. 25, 26), we fail at many sites to find PPFD levels where sunny and cloudy conditions coexist. Nevertheless, the result generally indicates larger midday ΔGPP than those in the morning and afternoon, although the modeled ΔGPP is slightly smaller than the observation. It should be noted that this large difference is captured by both ORCHIDEE_DF and ORCHIDEE trunk (Fig. 7b). Because direct and diffuse light are not distinguished in ORCHIDEE trunk, this midday ΔGPP should be mainly contributed by environmental factors other than the diffuse light fraction. The underestimation of midday ΔGPP could be a result of error in the current ORCHIDEE parameterizations. The ΔLUE derived by ORCHIDEE_DF also shows a largely similar distribution to that in the observations, but ORCHIDEE trunk underestimates the morning and afternoon ΔLUE (Fig. 7d, e, f).

3.4 Interactions between diffuse light and environmental factors

As implied by Fig. 7, the diffuse light fraction is not the only factor causing ΔGPP; other possible factors include temperature and the VPD (Gu et al., 2002; Cheng et al., 2015; Li et al., 2014). Thus, we investigate the diffuse light effect along temperature and VPD gradients in Fig. 8. To remove the effect of PPFD level, we only show ΔLUE.

Figure 8 The dependence of the LUE difference between cloudy and sunny conditions on climate factors. In the observations (blue), ORCHIDEE trunk (red) and ORCHIDEE_DF (green), the average and error bars indicate statistics of site-level means. (a) The dependence of the LUE difference on temperature. (b) The dependence of the LUE difference on VPD. Panels (c) and (e) are the same as panel (a) but for forest sites and short vegetation (grasslands and croplands) sites only. Panels (d) and (f) are the same as panel (b) but for forest sites and short vegetation sites only.

ΔLUE shows a unimodal curve along the temperature gradient for the observations and the two models (Fig. 8a). At low temperature, both models indicate a very low ΔLUE of 1 gC  µmol⁻¹, which is about one-third of the ΔLUE derived from observations. With increasing temperature, the observed ΔLUE shows a maximum at 10–20 ^∘C, with a magnitude of $\sim \mathrm{8}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹ and declines slightly at higher temperatures. The peak of ΔLUE simulated by ORCHIDEE_DF has a magnitude comparable to that of observations but at a higher temperature (20–25 ^∘C) than for the observations. The ΔLUE simulated by ORCHIDEE trunk is much smaller, with a peak of $\sim \mathrm{4}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹ at 10–15 ^∘C.

The effect of the VPD on ΔLUE is shown in Fig. 8b. For observations and both model simulations, a monotonic decreasing trend of ΔLUE along the VPD gradient is found. The ΔLUE from observations and ORCHIDEE_DF show a comparable magnitude, from $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹ at VPD < 0.5 kPa to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹ at a 2–4 kPa VPD level. The ΔLUE simulated by ORCHIDEE trunk is smaller than observations.

Apart from environmental factors, the effects of diffuse light may also differ among PFTs because different PFTs have different canopy structures and photosynthetic parameters (e.g., Vcmax). Here, we analyzed the ΔLUE in forests and short vegetation (grasslands and croplands) separately (Fig. 8c, d, ,e, f). In forests (Fig. 8c, d), ORCHIDEE_DF underestimates ΔLUE at temperatures lower than 20 ^∘C; it also largely captures the observed ΔLUE trend with VPD. ORCHIDEE trunk, in comparison, underestimates ΔLUE in all cases. Compared with forests, observations in short vegetation (Fig. 8e, f) show a stronger decline of ΔLUE at high temperatures (> 25 ^∘C) and high VPD conditions (> 0.5 kPa). However, for ORCHIDEE_DF, the short vegetation ΔLUE remains as high as for forests.

Figure 9 shows the distribution of ΔLUE in the temperature–VPD dimensions. Observations indicate that the largest ΔLUE is reached under conditions when temperature is in the range from 5 to 20 ^∘C and the VPD < 1 kPa (Fig. 9a). This temperature is thought to be more favorable for photosynthesis, as it is generally consistent with the photosynthesis optimum temperature detected by Huang et al. (2019) at latitudes where most of the sites are located. Under these conditions, the ΔLUE is usually over $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹. When the temperature is lower than 5 ^∘C or higher than 20 ^∘C or when the VPD becomes larger than 1 kPa, ΔLUE tends to decline. Compared with observations, the ΔLUE simulated by ORCHIDEE_DF shows a similar decreasing trend with VPD at all temperature levels (Fig. 9c); however, no obvious decline in ΔLUE is found at high temperatures. The ΔLUE simulated by ORCHIDEE trunk is much smaller than that for the observations (Fig. 9b).

Figure 9 The distribution of the LUE difference between cloudy and sunny conditions (ΔLUE) in the temperature–VPD field. The upper numbers in each grid square indicate the average of the site-level ΔLUE, while the numbers in parentheses indicate the number of sites with valid data. Panel (a) shows the ΔLUE based on observations, panel (b) shows the ΔLUE based on ORCHIDEE trunk, and panel (c) shows the ΔLUE based on ORCHIDEE_DF.

The ΔLUE from forests and short vegetation are shown separately in Fig. 10. Based on site-level observations (Fig. 10a), both vegetation types show a larger ΔLUE at lower VPD values between 5 and 20 ^∘C. In forests, there is also large ΔLUE under high temperature conditions, which mainly occurs in tropical forests (Fig. S3). Nevertheless, ORCHIDEE_DF still overestimates the ΔLUE at high temperatures (Fig. 10e), which is mainly due to the overestimation of ΔLUE at high temperatures for temperate forests (Fig. S3).

Figure 10 Same as Fig. 9 but for forests (a, c, e) and short vegetation (b, d, f).

Compared with forests, short vegetation shows a much stronger decline in ΔLUE at higher VPD levels (Fig. 10b); however, it is not well captured by ORCHIDEE_DF (Fig. 10f). In most cases, ORCHIDEE trunk tends to strongly underestimate ΔLUE unless the observed ΔLUE is small or negative (e.g., VPD > 2 kPa for short vegetation).

4 Discussion

4.1 Improvement of ORCHIDEE_DF

The role of diffuse light on photosynthesis has been found and modeled in different vegetation types (Gu et al., 2003; Niyogi et al., 2004; Misson et al., 2005; Alton et al., 2007b; Knohl and Baldocchi, 2008; Mercado et al. 2009; Oliphant et al., 2011; Kanniah et al., 2013; Williams et al., 2014; Cheng et al., 2015; Wang et al., 2018). However, very few studies have attempted to account for the diffuse light effect in a global land surface model, and fewer studies have used large FLUXNET datasets for evaluation. Here, by using flux observations from 159 sites over the globe, we show that – by separating the direct and diffuse light – ORCHIDEE_DF improves the simulation of GPP under sunny conditions and, more importantly, reproduces the observed impacts of diffuse light on GPP and LUE for most of the PFTs (Figs. 4, ,5, 6). Under cloudy conditions, ORCHIDEE_DF seems to perform slightly worse than ORCHIDEE trunk (Fig. 4). However, it should be noted that ORCHIDEE_DF has not been recalibrated, and all parameters are those from ORCHIDEE trunk despite the substantial changes in the code with respect to light partitioning and canopy light transmission. Furthermore, the GPP simulated by ORCHIDEE trunk shows different GPP SD biases under sunny and cloudy conditions, whereas ORCHIDEE_DF gives a more systematically underestimated GPP SD, which should be more easily corrected in a future calibration. The site-level comparison (Fig. 5) also explains how ORCHIDEE_DF reproduces the GPP increase compared with ORCHIDEE trunk. At most sites, the GPP simulated by the two models shows a similar magnitude under cloudy conditions, but the GPP simulated by ORCHIDEE_DF is significantly smaller under sunny conditions. This is because all light is considered to be diffuse light and is evenly distributed in each leaf layer in the one-stream canopy light transmission model in ORCHIDEE trunk. This simplified approach to the modeling of the light distribution leads to larger GPP under sunny conditions because the effect of light saturation on sunlit leaves is ignored. As ORCHIDEE trunk was calibrated using both sunny and cloudy data, but ORCHIDEE_DF corrected the overestimation under sunny conditions, ORCHIDEE_DF may give an overall underestimation using current parameters.

4.2 Factors affecting the response of GPP to diffuse light

Although diffuse light can increase the photosynthesis of shaded leaves, the GPP increase under cloudy conditions is not contributed solely by this effect. A recent field study suggested that photosynthesis from part of the canopy (especially sunlit leaves) benefits from the lower VPD rather than the higher diffuse light fraction under cloudier conditions (Wang et al., 2018). Our results show that higher diffuse PAR fraction is the main factor causing larger GPP under cloudy conditions compared with sunny conditions during the morning and the afternoon, as only ORCHIDEE_DF reproduced the observed positive ΔGPP during the two periods (Fig. 7). At midday, in comparison, the larger GPP under cloudy conditions should be mainly due to lower temperature or VPD rather than due to diffuse light, because ORCHIDEE trunk, which does not simulate the diffuse light effect, also reproduces this effect (Fig. 7). A similar effect is also reported by Cheng et al. (2015), who found that the midday GPP increase under cloudier conditions in croplands is mainly caused by lower temperature and lower VPD rather than by diffuse light. Photosynthesis is often considered to be limited by either carboxylation or electron transportation (Farquhar et al., 1980). It is when the shaded leaf photosynthesis is limited by light that diffuse light can increase GPP. At midday, large VPD values may cause stomatal closure, leading to a carboxylation-limited photosynthesis. Our results imply that it might be important to consider the diurnal cycle of environmental factors to better understand the effect of diffuse light.

It should be noted that the covariation of environmental factors with more diffuse light under cloudier conditions does not always benefit photosynthesis. For instance, if vegetation is cold stressed under cooler conditions, the decrease in temperature under cloudier conditions may strengthen this stress and offset the effect of diffuse light. Our analyses indicate that the effect of diffuse light on photosynthesis is weakened under most stressed conditions (Figs. 9, 10).

Another important factor is the light itself. When there is no light saturation of shaded leaves, under the same diffuse light fraction, stronger light levels are likely to benefit the shaded leaves more, resulting in higher ΔGPP (Fig. 5; ΔGPP tends to be larger at higher PPFD level at most sites). Nevertheless, apart from GPP, we also investigated the LUE (the photosynthesis per unit PPFD), which has removed this effect.

Besides environmental factors, canopy structure is also very important. Theoretically, thicker canopies with a large LAI tend to be more sensitive to diffuse light because a larger fraction of leaves are light limited due to shading (Fig. 1). As expected, ORCHIDEE_DF has shown an increasing ΔLUE with LAI (Fig. 11). However, the analyses based on LAI observations suggested a very weak positive effect of LAI on ΔLUE (Fig. 11). This insensitive response of ΔLUE to the LAI detected here should be treated with caution because the LAI observations are not well defined (maximum or average) and remain very limited in the current FLUXNET dataset (less than 10 in each LAI interval). Using more detailed LAI and CO₂ flux observations, Wohlfahrt et al. (2008) clearly exhibited the influence of LAI on diffuse light-induced photosynthesis changes in a grassland.

Figure 11 Same as Fig. 8a but for the LAI.

4.3 Uncertainty and limitations

Many empirical methods have been proposed to partition solar radiation into diffuse and direct light (e.g., Spitters et al., 1986; Weiss and Norman, 1985; Erbs et al., 1982). However, biases remain in the predicted diffuse light fraction under all aerosol and cloud conditions, which inevitably introduce some uncertainties to our analyses. Nevertheless, such methods are currently the most feasible approach at the flux site level. More continuous measurements of direct and diffuse surface radiation at more sites are desirable.

Table 1 Variables in this study.

Another source of uncertainty is from the light transmission model. In ORCHIDEE_DF, we used a two-stream radiative transfer approximation. In this model, the canopy trait parameters such as leaf scattering, leaf orientation and leaf clumping factors are assumed to be the same for all PFTs; however, real canopies are very diverse (Smith et al., 2004). In situ observations are required to obtain better parameters. Furthermore, the validity of the light transmission model in ORCHIDEE_DF depends on several assumptions described in Sect. 2.1. These assumptions are not always valid; for example, because direct solar beams are not exact parallels, leaves in canopies are not always sunlit or shaded, and they may also fall in penumbra regions, (i.e., regions where only part of the incoming direct solar beams are blocked; Smith et al., 1989; Cescatti and Niinemets, 2004). These more complex processes should be considered in future model development. Nevertheless, our simplified light transmission already succeeds at reproducing the observed diffuse light impact.

There are other sources of uncertainties in complex land surface models. Although ORCHIDEE_DF reproduces the magnitude of the diffuse light effects, it fails to reproduce the response of ΔLUE to temperature. For all PFTs, ORCHIDEE_DF underestimates the ΔLUE at low temperatures, and it overestimates ΔLUE at high temperatures (Fig. 8). The low temperature underestimation is also found in ORCHIDEE trunk, indicating that the models may have underestimated the tolerance of plants to low temperatures, whereas ORCHIDEE_DF tends to underestimate the impact of heat stress at high temperatures. This bias might be due to the parameterization of temperature acclimation which is based on observations mainly from a narrow temperature range (11–29 ^∘C) (Kattge and Knorr, 2007). For short vegetation, the introduction of diffuse light into the model results in an increase in ΔLUE at high temperatures and high VPD (Figs. 8, 10), indicating that the vegetation simulated by ORCHIDEE trunk remains light limited under such conditions. However, the strong decreasing trend of the observed ΔLUE along temperature and VPD gradients indicates heat and VPD stress. This implies that parameters in the current ORCHIDEE version may have underestimated the response of grassland and cropland photosynthesis to heat and VPD stress.

Besides the possible bias in parameters, both ORCHIDEE trunk and ORCHIDEE_DF lack a representation of the response of leaf temperature to radiation. Instead, the air temperature is used directly to represent the leaf temperature throughout the canopy for simulating gas exchange processes in the current model. As shown by Chen and Zhuang (2014), the changes in the radiation regime due to aerosols can significantly affect leaf temperature, which could potentially affect GPP. For now, ORCHIDEE_DF remains incapable of dealing with this response of leaf temperature. Further developments are needed to disentangle the role of leaf temperature and diffuse light on GPP.

5 Conclusion

In this study, we added a module to partition the downward surface solar radiation into diffuse and direct components as well as a new canopy radiative transfer model, which separates the existing multilayer canopy into sunlit and shaded leaves, to the ORCHIDEE trunk model. The resulting new land surface model, ORCHIDEE_DF, is evaluated using the La Thuile flux dataset over 159 sites comprising 11 PFTs. Compared with ORCHIDEE trunk, ORCHIDEE_DF improves the GPP simulation under sunny conditions. This improvement successfully reproduces the observed enhancement of GPP under cloudier conditions at most of the sites.

Using observed and modeled GPP, we found an increase in GPP under cloudier conditions at all times of the day; however, the mechanisms causing this effect are different at midday from the morning and afternoon. During morning and afternoon, the increase in GPP is mainly caused by increased diffuse light fraction, whereas the GPP increase is mainly due to weaker stress from temperature and VPD at midday.

Observations indicate that the maximum LUE difference can be over $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{8}}$ gC µmol⁻¹ under cloudy and sunny conditions for the same light level. The maximum LUE is found at temperature and VPD conditions more favorable for photosynthesis (5–20 ^∘C for temperature and < 1 kPa for VPD). With increasing VPD or under lower or higher temperatures, the LUE may decrease. Compared with observations, ORCHIDEE_DF underestimates the diffuse light effect at low temperature and overestimates it at high temperatures, possibly due to imperfect temperature acclimation parameterization in the current ORCHIDEE model. In grasslands and croplands, ORCHIDEE_DF overestimates the diffuse light effect on LUE, which might be due to an overestimation of their tolerance to dry conditions.

As ORCHIDEE_DF is a land surface model which is able to capture the effect of diffuse light for a large number of sites over the globe, we are confident that, with this improved model framework and proper calibration, we can investigate the effect of aerosols on global biogeochemical cycles and assess the impact of aerosol emission policies and aerosol-related climate engineering on such cycles.

Appendix A: List of acronyms

Fdf	Fraction of diffuse radiation
GPP	Gross primary production
LAI	Leaf area index
LSM	Land surface model
LUE	Light use efficiency
NIR	Near-infrared radiation
PAR	Photosynthetically active radiation
PFT	Plant functional type
PPFD	Photosynthetic photon flux density
SW	Downward shortwave radiation at the top of canopy
TOA	Top of atmosphere
TOC	Top of canopy
VPD	Vapor pressure deficit