On the impact of canopy model complexity on simulated carbon, water, and solar-induced chlorophyll fluorescence fluxes

Lack of direct carbon, water, and energy flux observations at global scales makes it difficult to calibrate land surface models (LSMs). The increasing number of remotesensing-based products provide an alternative way to verify or constrain land models given their global coverage and satisfactory spatial and temporal resolutions. However, these products and LSMs often differ in their assumptions and model setups, for example, the canopy model complexity. The disagreements hamper the fusion of global-scale datasets with LSMs. To evaluate how much the canopy complexity affects predicted canopy fluxes, we simulated and compared the carbon, water, and solar-induced chlorophyll fluorescence (SIF) fluxes using five different canopy complexity setups from a one-layered canopy to a multilayered canopy with leaf angular distributions. We modeled the canopy fluxes using the recently developed land model by the Climate Modeling Alliance, CliMA Land. Our model results suggested that (1) when using the same model inputs, model-predicted carbon, water, and SIF fluxes were all higher for simpler canopy setups; (2) when accounting for vertical photosynthetic capacity heterogeneity, differences between canopy complexity levels increased compared to the scenario of a uniform canopy; and (3) SIF fluxes modeled with different canopy complexity levels changed with sun-sensor geometry. Given the different modeled canopy fluxes with different canopy complexities, we recommend (1) not misusing parameters inverted with different canopy complexities or assumptions to avoid biases in model outputs and (2) using a complex canopy model with angular distribution and a hyperspectral radiation transfer scheme when linking land processes to remotely sensed spectra.


Introduction
Land surface models (LSMs) simulate the carbon, water, and energy fluxes at the land-atmosphere interface at regional and global scales and are a key component for Earth system models (ESMs). The ability of LSMs to accurately model the carbon, water, and energy fluxes within vegetation canopy largely determines the predictive skills of the ESMs. Modeling canopy carbon, water, and energy fluxes dates back to the early 20th century, and various canopy models have different complexities from a single layer to multiple layers (see Bonan et al., 2021, for an overview). To date, the most widely used canopy models in the LSM community are the "big-leaf model family".
It should be noted that a big-leaf model may refer to different models within the last decades given their interchangeable uses (Luo et al., 2018). According to Luo et al. (2018), the big-leaf model can be categorized at least as the following types given the purposes for which they were developed. (1) The one-big-leaf canopy model considers a canopy to be a single big leaf and was typically used with the Penman-Monteith equation (Penman, 1948;Monteith, 1965) to compute land surface evaporation in early LSMs. Sellers et al. (1992) updated the one-big-leaf model by adding an exponentially diminishing photosynthetic rate within the canopy depth to upscale photosynthesis for the carbon-water coupled LSMs. Yet, this scheme often underestimated the canopy assimilation rate, as the exponential function cannot properly represent the vertical light and photosynthesis profiles.
(2) A two-leaf radiation scheme separates the canopy into a group of sunlit leaves and a group of shaded leaves (Norman, 1982;De Pury and Farquhar, 1997;Campbell and Y. Wang and C. Frankenberg: Canopy model complexity Norman, 1998;Chen et al., 1999) and was used to account for the horizontal and vertical light heterogeneity in the canopy.
(3) A two-big-leaf canopy model combines the one-big-leaf canopy model and two-leaf radiation scheme to upscale carbon and water fluxes and treats each of the sunlit and shaded fractions as a single big leaf, where leaf biochemical parameters and radiation are upscaled to the canopy level (De Pury and Farquhar, 1997;Wang and Leuning, 1998). (4) A two-leaf canopy model uses a two-leaf radiation scheme and treats each of the sunlit and shaded fractions as a leaf with average traits for its representation (not an integrated value as in a big leaf) (Chen et al., 1999Sprintsin et al., 2012). Therefore, one needs to be cautious when using the term bigleaf model, as it may refer to (i) a two-leaf radiation scheme which is a canopy radiative transfer model or (ii) upscaling schemes which differ in the way leaf biochemical parameters are integrated (such as the one-big-leaf and two-big-leaf models) or averaged (such as the two-leaf canopy model).
The biggest advantage of the big-leaf model family is computational efficiency given the simple mathematical formulation. The potential disadvantages of the big-leaf model family are also obvious; e.g., the model is too simplified and thus not able to resolve vertically varying profiles and microclimates in the canopy, such as air temperature, humidity, and wind speed (Bonan et al., 2021). Thus, there is an increasing demand for LSMs to move from a simple one-layered canopy to a multi-layered one.
One of the most important functions of canopy models is to predict carbon, water, and energy fluxes globally in the future to determine whether the land will remain a carbon sink. Though canopy models with different complexity levels have been extended to a global scale in different LSMs, researchers are facing a key problem: a lack of direct globalscale carbon, water, and energy flux observations. The lack of data makes it difficult to calibrate the LSMs at global scales, particularly those using more complex canopy setups given the more parameters required. As a result, though it is shown that a multi-layered canopy model better resolves energy fluxes in the canopy (Bonan et al., 2021), little is known about whether the multi-layered canopy models show improved predictive skills (particularly in terms of carbon and water fluxes) compared to the big-leaf models which are widely used in existing LSMs.
To better constrain LSMs with data, people realized the promise of remote sensing data given their global coverage and satisfactory spatial and temporal resolutions. Regarding carbon, research has shown that solar-induced chlorophyll fluorescence (SIF) and the near-infrared reflection of vegetation (NIRv) are correlated with plant gross primary productivity (GPP; Frankenberg et al., 2011;Zhang et al., 2016;Sun et al., 2018;Badgley et al., 2019). Regarding water, researchers also found SIF to be useful for inverting the transpiration rate by prescribing stomatal responses to the environment (e.g., Shan et al., 2021) and the vegetation optical depth to be useful in sensing aboveground biomass and canopy water stress (Momen et al., 2017;Zhang et al., 2019). Regarding energy, various models and algorithms have been used to detect the surface energy balance using optical light and microwaves (e.g., Roerink et al., 2000;Norman et al., 2003). Further, methods and applications have been developed to invert plant traits from remote sensing data, such as leaf area index (e.g., Colombo et al., 2003;Deng et al., 2006) and chlorophyll content (Croft et al., 2020).
Despite the increasing number of inverted fluxes and plant trait datasets, there is limited research into testing the capability of these data in improving LSM predictions. Among the various reasons that hamper the fusion of large-scale datasets into LSMs, incompatibility between model and data assumptions seems to be the major reason. For example, the disagreement in canopy complexity may introduce errors into modeling if one uses the data inverted from a canopy complexity level (e.g., one-layered canopy) in a model with a different canopy complexity level (e.g., multilayered canopy). Further, the flux and trait maps inverted from remote sensing data often use simplified plant physiological representations, which are, however, key processes in land modeling. For example, studies that derive GPP from SIF or NIRv often assume linear correlation between them, whereas vegetation models must account for light saturation (Zhang et al., 2016).
Ideally, LSMs can be constrained using raw reflection and fluorescence spectra. This, nevertheless, requires the LSMs to move from broadband canopy radiation to a hyperspectral representation and from sunlit and shaded fractions to leaf angular distributions (such as the land model developed by Climate Modeling Alliance, CliMA Land; Wang et al., 2021b). This way, the LSM can be directly coupled to remotely sensed canopy spectra (e.g., Shiklomanov et al., 2021) rather than to reprocessed datasets using often incompatible assumptions. The increasing canopy complexity, however, comes with high costs: (a) many more computational resources are required by the increasing number of leaves (e.g., CliMA Land canopy has a default of 6500 leaves per tree in the canopy, whereas a two-leaf canopy has two "leaves"); (b) canopy radiation and fraction (e.g., the CliMA Land model calculates the radiation and fraction based on leaf angular distribution for a default of 6500 leaves) are more complicated; and (c) most importantly, there is increasing difficulty for research communities when understanding or using the model.
To resolve the problems of a complicated canopy, we examined by how much carbon, water, and SIF fluxes may differ when using different canopy complexity representations in the CliMA Land model, spanning from a one-layered canopy to a multi-layered canopy with hyperspectral radiation and leaf angular distributions (Fig. 1). With the model simulations, we were able to answer the following questions: (1) how does canopy complexity impact modeled canopy fluxes, and (2) could data inverted using different canopy complexity levels be compatible? Regarding the ease of un- Figure 1. Canopy complexity levels. 1X: single-layer canopy without sunlit or shaded fractions. 2X: single-layer canopy with sunlit and shaded fractions. KX: multiple-layer canopy without sunlit or shaded fractions. 2KX: multiple-layer canopy with sunlit and shaded fractions per layer. IJKX: multiple-layer canopy with sunlit and shaded fractions per layer, with the sunlit fraction being further partitioned based on leaf inclination and azimuth angular distributions.
derstanding and using an LSM with various canopy complexities, we presented and suggested the highly modularized CliMA Land model, which can be easily set up to simulate canopy fluxes using different canopy complexity levels.

Materials and methods
We used the CliMA Land model (v0.1) to evaluate how canopy model complexity impacts the simulated carbon, water, and SIF fluxes. The CliMA Land model mechanistically addresses soil-plant-air continuum processes and is able to simulate canopy carbon and water fluxes as well as SIF simultaneously (Wang et al., 2021b). CliMA Land model code and documentation are freely and publicly available at https: //github.com/CliMA/Land (last access: 15 November 2021).

Canopy complexity levels
Leaf physiological responses to light are highly non-linear, such as stomatal conductance to water vapor (g sw ) and net photosynthetic rate (A net ). Typically, when absorbed photosynthetically active radiation (APAR) is low, g sw and A net increase with higher APAR (Fig. 2a); when APAR is high, g sw and A net saturate. If one has a leaf with low APAR and a leaf with high APAR (e.g., closed circles on the solid curves of Fig. 2a), the mean behavior of the two leaves ought to be the average g sw and A net values (closed circles on the colored dashed lines of Fig. 2a). However, if one uses the mean APAR of the two leaves and calculates g sw and A net based on the mean APAR, g sw and A net would be overestimated (open circles on the colored solid curves of Fig. 2a).
Note that averaging APAR values that are beyond the turning point, say 350 µmol m −2 s −1 , may not result in any bias in modeled g sw and A net (such as for sunlit and shaded leaves in the top canopy layer); however, averaging APAR for leaves with high APAR and low APAR, say 300 and 50 µmol m −2 s −1 , would result in overestimated g sw and A net (such as for shaded leaves in an upper and lower canopy as typically done in the two-leaf radiation scheme). Thus, an overly simplified canopy model may overestimate canopylevel carbon and water fluxes, because of the inappropriately averaged APAR, as f (x) = f (x) when averaging non-linear functions (A net (APAR) > A net (APAR) in leaf photosynthesis).
To evaluate how much canopy model complexity matters, we modeled the canopy using five different levels of complexity, and they are denoted as "1X", "2X", "KX", "2KX", and "IJKX" (Fig. 1). 1X represents the scenario in which the canopy is treated as a single average leaf without sunlit or shaded fractions, and leaf radiation is averaged for the entire canopy. 2X complicates 1X by partitioning the average leaf to sunlit and shaded fractions. KX enhances 1X by partitioning the canopy to multiple layers (but no sunlit or shaded fractions per layer). 2KX partitions each canopy layer of KX to sunlit and shaded fractions. IJKX further modifies 2KX by accounting for leaf inclination and azimuth angle distributions per layer (Fig. 1). See Table 1 for the canopy model complexity adopted by other vegetation models (see https://yujie-w.github.io/PAGES/dev/methods/#Vegetationcanopy-model-complexity for a growing list).
For IJKX, we simulated the canopy radiative transfer using the CliMA Land-adapted SCOPE model (Soil Canopy Observation, Photochemistry and Energy fluxes; SCOPE v1.7; Yang et al., 2017). The adaptations included that carotenoid absorption was accounted for as APAR (Wang et al., 2021b) and that canopy clumping was addressed using a clumping index (Braghiere et al., 2021). At layer i, the shaded-leaf fraction (relative to total leaf area in the canopy) is p sh,i , and the sunlit-leaf fraction (relative to total leaf area in the canopy) is p sl,i (incl, azi) ("incl" is the inclination angle, and "azi" is the azimuth angle). The fraction of sunlit leaves relative to total canopy leaf area in a given canopy layer is computed as incl,azi where LAI is total leaf area index, LAI(i) is the leaf area index above the ith canopy layer, LAI(i + 1) is the leaf area index in and above the ith canopy layer, k is the extinction coefficient as a function of leaf inclination angle distribution, and is the clumping index. Then p sl,i (incl,azi) is computed 32 Y. Wang and C. Frankenberg: Canopy model complexity Figure 2. Non-linear leaf responses to the environmental and physiological parameters. (a) Responses of stomatal conductance to water vapor (g sw ; cyan solid curve) and net photosynthetic rate (A net ) to absorbed photosynthetically active radiation (APAR). The black dotted vertical lines indicate two leaves at low-and high-light conditions. Mean behavior of the two leaves ought to be the closed circles on the colored dotted lines. However, using mean APAR for the leaves would result in overestimated g sw and A net (open circles). (b) Non-linear g sw and A net responses to leaf photosynthetic capacity, represented by the maximum carboxylation rate (V cmax ). Environmental and leaf physiological settings for the simulations are the following: air and leaf temperatures at 298.15 K, atmospheric vapor pressure at 1500 Pa (relative humidity at 0.47), atmospheric CO 2 partial pressure at 40 Pa, atmospheric pressure at 101 325 Pa, V cmax (for panel a) at 60 µmol m −2 s −1 , and maximal stomatal conductance at 0.3 mol m −2 s −1 .
where I is the number of inclination angles and J is the number of azimuth angles. The fraction of shaded leaves relative to total canopy leaf area in a given canopy layer is computed as Corresponding APAR values for the shaded and sunlit leaves are APAR sh and APAR sl,i (incl, azi), respectively. We used default values of I = 9 inclination angles, J = 36 azimuth angles, and K = 20 vertical layers for IJKX (K = 20 for 2KX and KX as well).
The 2KX fraction and APAR were derived from IJKX by weighing APAR for sunlit leaves per canopy layer: The KX fraction and APAR were derived from 2KX by weighing APAR for all sunlit and shaded leaves per canopy layer: The 2X fraction and APAR were derived from 2KX by weighing APAR for sunlit and shaded leaves for all canopy layers, respectively, the following: 1X APAR was derived from KX by weighing APAR for all layers: We emphasize here that to derive canopy fluorescence spectrum and its sun-sensor geometry, we need to simulate the canopy radiative transfer using hyperspectral reflectance, transmittance, and fluorescence. Due to the high spectral resolution and multiple layers required, radiative transfer and canopy fractions in complex models such as SCOPE are computed numerically. In comparison, radiative transfer and sunlit/shaded fractions are computed analytically in the twoleaf radiation scheme, as the model is single layered and uses broadband reflectance and transmittance (Campbell and Norman, 1998;Bonan et al., 2021). Yet, the two-leaf radiation Figure 3. Comparison of profiles of mean absorbed photosynthetically active radiation (APAR) for four different canopy complexity levels. 1X: single-layer canopy without sunlit or shaded fractions. 2X: single-layer canopy with sunlit and shaded fractions. KX: multiple-layer canopy without sunlit or shaded fractions. 2KX: multiple-layer canopy with sunlit and shaded fractions per layer. IJKX: multiple-layer canopy with sunlit and shade fractions per layer, with the sunlit fraction being further partitioned based on leaf inclination and azimuth angular distributions. The abbreviations "sl" and "sh" stand for sunlit and shaded leaves, respectively. schemes that use broadband radiative transfer are not adequate for accurate fluorescence modeling. Crucially, the difference in the analytic and numerical solutions could result in biases in the simulated APAR and fraction. To avoid such bias, we computed APAR and sunlit/shaded fractions for the simpler canopy setups numerically using the algorithm in IJKX. See Fig. 3 for the APAR profiles for 2KX, KX, 2X, and 1X derived from IJKX. Also, we note here that leaf biochemical parameters and APAR were not integrated within a canopy layer or sunlit/shaded fractions; instead, we used average APAR and leaf traits in our simulations. Thus, our 1X model is a one-leaf model rather than a one-big-leaf model, and our 2X model resembles the two-leaf model rather than the two-big-leaf model.

Vertical canopy profile
Leaf traits in the canopy are not uniform among the canopy layers. Typically, leaf photosynthetic capacity (usually represented by the maximum carboxylation rate at 25 • C, V cmax ) is higher in upper canopy because of the better light environment. Further, leaf physiological responses to V cmax are also highly non-linear, and using average V cmax may also result in overestimated g sw and A net and thus carbon and water fluxes (e.g., shift from solid circles to open circles in Fig. 2b).
To examine how much the vertical V cmax profile impacts modeled canopy flux simulations, we ran the model simula-34 Y. Wang and C. Frankenberg: Canopy model complexity tion in two scenarios, one using uniform V cmax in the canopy and one using decreasing V cmax towards the lower canopy. For the latter scenario, V cmax at layer i was tuned using an exponential function following De Pury and Farquhar (1997) and Chen et al. (2012): where V cmax,top is V cmax at the top of the canopy and 0.15 is the shape factor that describes the decreasing V cmax with canopy depth. Note that as leaves are experiencing dynamically changing light environment throughout the day, it is unrealistic to assume the sunlit and shaded leaves have different traits; thus, we only accounted for the vertical heterogeneity but neglected the horizontal heterogeneity in each canopy layer, namely using the same characteristics for leaves within the same canopy layer. The V cmax profile was applied to IJKX, 2KX, and KX directly, whereas weighed mean V cmax = i V cmax,i · KX p i was used in 1X. The V cmax profile or value stayed constant in these four scenarios throughout the simulation, as sunlit/shaded fractions did not impact them. We note here that mean V cmax changed with sunlit/shaded fractions in 2X  and particular averages of V cmax for sunlit and shaded fractions ( 2X V cmax,sl and 2X V cmax,sh , respectively) need to be updated with sunlit and shaded fractions: Note that we tuned the maximum electron transport rate and leaf respiration rate in the same manner as V cmax .

Canopy flux simulations
We simulated the canopy carbon and water fluxes using a stomatal optimization model developed in Wang et al. (2020) given the good model performance and scalability (Wang et al., 2021a). The stomatal optimization model posits that stomatal opening is optimized when the difference between carbon gain and risk is maximum: where E is the leaf transpiration rate and E crit is the critical transpiration rate of the leaf beyond which leaf hydraulic conductance drops below 0.1 % of the maximum (see Sperry et al., 2016, andWang et al., 2021b, for more details of E crit ). At each canopy complexity level, for a given environmental condition set, we were able to obtain the steady-state stomatal conductance for each APAR, from which we computed steady-state A net using the classic C3 photosynthesis model (Farquhar et al., 1980) and E as well as leaf fluorescence quantum yield (φ F ) using the model developed in . Stand-level carbon flux, namely net ecosystem exchange (NEE; normalized per ground area), was computed using NEE = LAI · (A net · p) − R remain : where LAI is the leaf area index and R remain is the ecosystem respiration rate per ground area excluding the leaves. The transpiration rate from the canopy is computed and used as a proxy for estimating the difference in model ecosystem evapotranspiration (ET; normalized per ground area) using ET ≈ LAI · (E · p): We remind the reader here that soil evaporation is a function of soil water content, soil surface temperature, and atmospheric-vapor-pressure deficit and that soil evaporation should be the same for all tested canopy complexity models; this is also the case for evaporation from intercepted water on plant surface. Therefore, the modeled ET difference is 100 % caused by canopy transpiration, and using transpiration would not result in any biases in the relative difference of modeled ET.
For IJKX, we used φ F computed for each sunlit and shaded leaf at each layer to compute the canopy-level SIF spectrum.
For 2KX, we plugged the φ F calculated for sunlit fraction into all the sunlit leaves of the corresponding layer of IJKX and the shaded φ F into the shaded leaf of the corresponding layer of IJKX. Then we re-simulated the SIF spectrum at IJKX and used it as that of 2KX. For KX, we plugged the φ F calculated for the whole layer into all the leaves of the corresponding layer of IJKX and recalculated the SIF spectrum. For 2X, we plugged the φ F of the sunlit fraction into all the sunlit leaves in IJKX and shaded φ F into all the shaded leaves in IJKX and recalculated the SIF spectrum. For 1X, we plugged the φ F into all the leaves in IJKX and recalculated the SIF spectrum. We compared SIF at 740 nm (SIF 740 ) among different complexity levels.
Despite the importance of vertical microclimate heterogeneity in modeled canopy energy fluxes (e.g., Bonan et al., 2021), we held environmental conditions constant among vertical canopy layers for all tested canopy complexities. Doing this allowed us to tease apart the impact of APAR distribution in the canopy (due to canopy complexity) on simulated carbon, water, and SIF fluxes.

Sensitivity analysis
We ran a sensitivity analysis to environmental cues for all five complexity levels to examine how much they differ in predicted carbon, water, and SIF fluxes. The tested cues included solar radiation, atmospheric-vapor-pressure deficit (VPD), temperature, soil water potential ( soil ), and atmospheric CO 2 partial pressure (P CO 2 ). When we altered temperature, we changed the air and leaf temperature at the same time and held air relative humidity (RH) constant at 0.47 (fraction; unitless). Saturated water vapor pressure was computed using the Clausius-Clapeyron equation: where P triple is the vapor pressure at the triple point (in Pa), T is the temperature (in K), T triple is the temperature at the triple point (in K), c p is the difference in isobaric specific heat of vapor and liquid (in J Kg −1 K −1 ), R v is the gas constant of water vapor (in J Kg −1 K −1 ), and LH v0 is the latent heat of vaporization at the triple point. Atmospheric vapor pressure was computed using P sat · RH. For each tested environmental cue, we changed only the tested cue while holding all other environmental conditions constant. We ran the sensitivity test in two scenarios: (a) V cmax was uniform throughout the canopy, and (b) V cmax decreased exponentially in the lower canopy. For the two scenarios, we let the entire-canopy mean V cmax be the same (namely mean V cmax at 1X). We compared the modeled site-level NEE, ET, and SIF 740 among canopy complexity levels.

Diurnal cycles
To evaluate how much the canopy complexity models differ in real-world simulations, we ran the model using weather data from a flux tower located at Ozark, Missouri, USA (US-MOz; Gu et al., 2016). We used the weather and soil moisture data from day 177 to 179 in the year 2019 and prescribed leaf temperature and soil water potential to maximally reduce uncertainty among model setups. Briefly, we used outgoing longwave radiation from flux tower measurements to invert canopy temperature and used it as leaf temperature; we also used soil water content to estimate soil water potential and used it as a boundary condition for the soil-plant-air continuum. Prescribing leaf temperature and soil water potential allowed us to tease apart the difference caused by canopy complexity from that caused by environmental and physiological differences. See Wang et al. (2021b) for the model setup details for US-MOz. In addition to the observations that were used to set up the CliMA Land model (Wang et al., 2021b), we further applied vertical V cmax profiles in the simulations (note that V cmax changed in the sunlit and shaded fractions with time for 2X and stayed constant for the other four complexity levels). We tuned V cmax and the whole-plant hydraulic conductance to let IJKX predict reasonable NEE and ET values and used these tuned parameters in all the tested canopy complexity levels. We note here that we were not trying to argue one complexity was better than others but to examine how much the complexity levels differ when we used exactly the same model input parameters. We compared the model-predicted carbon, water, and SIF fluxes. Note here that observed SIF depends on the sunsensor geometry and that SIF retrievals often have different sun-sensor geometries (e.g., the TROPOMI satellite; Köhler et al., 2018). Thus, it is necessary to examine how the sunsensor geometry may impact the SIF flux across canopy complexity levels. We ran the test using the weather data from (a) 12:00-12:30 and (b) 16:00-16:30 of day 177 in the year 2019. At each tested time window, we computed the theoretical SIF at 740 nm for a series of viewing zenith angles from 0 to 85 • and relative azimuth angles (angle between the sensor and sun) from 0 to 360 • . We compared how much 2KX, KX, 2X, and 1X differed from IJKX.
Given that averaging APAR theoretically results in overestimated carbon and water fluxes, we expected that the difference among different canopy complexity levels meets the following trends: (a) 1X > 2X > 2KX and (b) 1X > KX > 2KX. Further, as V cmax also theoretically results in overestimated carbon and water fluxes, we expected that adding a vertical V cmax profile further increases the difference in fluxes across canopy complexity levels. . (e) Responses to atmospheric CO 2 partial pressure (P CO 2 ). This sensitivity analysis was done assuming uniform photosynthetic capacity in the canopy.

Sensitivity analysis
When a uniform V cmax profile was applied, all tested five canopy complexity levels exhibited similar carbon and water flux responses to changing environmental cues (Fig. 4). The responses included increasing canopy photosynthesis and transpiration with higher radiation (Fig. 4a), increasing and then decreasing photosynthesis and increasing transpiration with a higher temperature (Fig. 4b), decreasing photosynthesis and increasing transpiration with a higher VPD (Fig. 4c), decreasing photosynthesis and transpiration with drier soil (Fig. 4d), and increasing photosynthesis and decreasing transpiration with higher atmospheric CO 2 partial pressure (Fig. 4e). Further, as expected, 1X, 2X, KX, and 2KX all overestimated canopy photosynthesis and transpi-ration compared to the IJKX mode; and the overestimation ratios met 1X > 2X > 2KX and 1X > KX > 2KX.
The SIF responses to changing environmental cues in general agreed in trends among tested complexity levels (Fig. 4). However, SIF responses to radiation, temperature, and atmospheric CO 2 differed dramatically among the five canopy complexity levels given the different response magnitudes (Fig. 4b and e). 1X and KX often resulted in different trends compared to IJKX (Fig. 4). 2X and 2KX overall tracked the SIF responses though slightly overestimated SIF of IJKX well. Notably, we found high disagreement between 2X and IJKX at intermediate radiation and increasing disagreement at higher atmospheric CO 2 (Fig. 4a and e).
When an exponential vertical V cmax profile (lower V cmax in the lower canopy) was applied when simulating canopy fluxes, we found similar trends compared to the scenario with constant V cmax (Fig. 5). The differences, however, were that all carbon, water, and SIF fluxes were lower when we ap- . (e) Responses to atmospheric CO 2 partial pressure (P CO 2 ). This sensitivity analysis was done assuming exponentially decreasing photosynthetic capacity in the lower canopy. plied a vertical V cmax profile (Fig. 5). Again, like the scenario of a uniform V cmax , we also found divergent SIF responses to radiation and increasing disagreements among 2X, 2KX, and IJKX for elevated CO 2 (Fig. 5a and e). The divergent flux responses to P CO 2 underlined the importance of adopting a more complex canopy in future land modeling given that (i) CO 2 concentration within the canopy airspace may change dramatically within a diurnal cycle due to plant carbon fixation and (ii) atmospheric mean CO 2 is increasing rapidly due to anthropogenic emissions.
2KX and 2X had a lower difference from IJKX compared to KX and 1X, and 2KX had the lowest error given the better-resolved APAR fractions (Figs. 4 and 5). Combining all response curves together from Fig. 4, we found that when V cmax was evenly distributed in the canopy, relative differences between 2KX and IJKX for carbon, water, and SIF fluxes were 2.4 %, 1.2 %, and 2.8 %, respectively (Fig. 6).
When accounting for a vertically heterogeneous V cmax profile, we still found a lower difference between 2KX and IJKX, and the relative differences were 11.1 %, 3.7 %, and 7.9 % (the differences for 2X were 23.4 %, 8.2 %, and 13.2 %; Fig. 6). Overall, 2KX had a relative error lower than 10 %. Further, the higher error when adopting a vertical V cmax profile agreed with our expectation as the impacts from APAR and V cmax added up (canopy fluorescence was lower for the simpler canopy model at low radiations; Fig. 6).

Diurnal cycle
Our model simulations suggest that all tested canopy complexity levels can qualitatively capture the trends of carbon and water fluxes at the tested flux tower site (Fig. 7). How- ever, the tested complexity levels differed dramatically in the magnitudes of carbon and water fluxes. In general, 1X had the highest fluxes for both carbon and water fluxes (represented by NEE and ET), followed by KX, 2X, 2KX, and IJKX ( Figs. 7 and 8). Though 2KX and 2X, in general, had relatively small differences from IJKX, we were still able to distinguish the difference (Figs. 7 and 8). We note here again that Figs. 7 and 8 were meant to highlight the difference between canopy complexity levels in model simulations but not to say that some models were better than others.

Sun-sensor geometry
Using less complicated canopy complexity (namely 2KX, KX, 2X, and 1X) impacted the observed SIF depending on the sun-sensor geometry (Figs. 9 and 10). For the tested time window at 12:00-12:30, 2KX has the least difference from IJKX, followed by 2X, KX, and 1X. In general, 2KX had a difference lower than 11 % at any viewing zenith angle or relative azimuth angle for the tested time window (Figs. 9). The impact of sun-sensor geometry changed with time because of changes in solar zenith angle and total radiation (e.g., at 16:00-16:30 in the afternoon; Fig. 10). While 2KX still had lower overestimated SIF compared to 2X, KX had better agreement with IJKX, and 1X even underestimated SIF. The dramatic changes in SIF from KX and 1X were due to lower incident radiation from 16:00 to 16:30 (Fig. 5a).

Fluorescence and radiation
While simpler canopy models in general predicted higher carbon, water, and SIF fluxes, there were some scenarios that the simpler models predict contrasting SIF responses compared to IJKX: (a) when total radiation increased, SIF of the simpler canopy models was lower than that of IJKX at low radiation but were higher than that of IJKX at high radiation (Figs. 4 and 5); (b) 1X model SIF increased and then decreased and stayed unchanged with higher atmospheric CO 2 (Figs. 4 and 5); and (c) KX model SIF increased marginally with higher atmospheric CO 2 for a canopy without a vertical V cmax gradient but decreased with higher CO 2 for a canopy with a vertical V cmax profile (Figs. 4 and 5). These contrasting patterns of the simpler models resulted from the different photosynthesis system II (PSII) quantum yield and fluorescence quantum yield (namely φ F ) responses to APAR and CO 2 (Fig. 11a, b). The PSII quantum yield measures the efficiency of converting absorbed photons to electrons by PSII; and φ F measures the efficiency of converting ab- sorbed photons to fluorescence photons. The PSII yield increases and then saturates with higher leaf internal CO 2 and lower APAR. In our model, φ F follows the parameterization of van der Tol et al. (2014) but fitted on leaves measured by Flexas et al. (2002), as first used in Lee et al. (2015). Typically, the PSII-to-φ F relationship depends on the state of non-photochemical quenching (NPQ; Porcar-Castell et al., 2014). φ F has a maximum at intermediate PSII levels (around 0.6) but decreases at lower PSII yields (increased NPQ) as well as higher PSII yields (increased competition with photochemical quenching). This general behavior explains what we see the following: φ F (a) stays unchanged at low radiation with higher leaf internal CO 2 , (b) increases and then decreases and stays unchanged with higher leaf internal CO 2 at intermediate APAR, (c) increases with higher CO 2 at high CO 2 , and (d) increases and then decreases with higher APAR (Fig. 11a, b). Though φ F in general agrees with the PSII yield patterns at high APAR (typical experimental and top-of-canopy scenarios), the disagreements at low APAR could result in problems when APAR is inappropriately averaged. In our case, the turnover from APAR regions in which PSII and φ F are anticorrelated (light-limited) to the region in which they are correlated (increase in NPQ) happens at around 200 µmol m 2 s −1 .
When total radiation was higher, the product of φ F and APAR (leaf-level SIF) increased (Fig. 11c). When φ F stayed unchanged at low APAR, leaf-level SIF increases linearly with higher APAR, and SIF increases faster when φ F starts to increase after a certain threshold (the threshold increased with higher leaf internal CO 2 ; Fig. 11b). Then leaf-level SIF slowed down with higher APAR due to decreasing φ F at higher APAR and was higher when leaf internal CO 2 was higher (Fig. 11b, c). As leaf internal CO 2 was theoretically lowest for 1X, then followed by KX, 2X, and 2KX given the way APAR was averaged, it was expected that 2KX increased earliest with higher APAR and that 1X had the highest SIF value at high radiation (Figs. 4a and 5a). Therefore, in the diurnal-cycle simulations, 1X SIF overestimated SIF at noon when radiation was high ( Fig. 9) but underestimated SIF in the late afternoon as a result of lower radiation (Fig. 10). The inconsistent SIF patterns at low and high radiation from simpler canopy models may potentially result in biases in modeled diurnal and seasonal SIF, and thus we suggest using a complex canopy model when possible to minimize the impact from heterogeneous canopy radiation and leaf physiology.
The 1X model SIF response to atmospheric CO 2 ought to depend on the mean canopy APAR (Fig. 11b): (a) if mean APAR was low, 1X SIF should stay constant with higher CO 2 ; (b) if mean APAR was moderate, 1X SIF ought to increase and then decrease and stay constant with higher CO 2 ; and (c) if mean APAR was high, 1X SIF would increase with higher CO 2 (Fig. 11b, d). For the simulations in Figs. 4 and 5, mean APAR was 156 µmol m −2 s −1 , and thus 1X SIF increased and then decreased with higher CO 2 .
The KX model SIF response to atmospheric CO 2 was impacted by both leaf internal CO 2 and the vertical APAR profile given the heterogeneous APAR. As φ F was higher in the middle layers at lower atmospheric CO 2 when there is no vertical V cmax gradient, modeled SIF showed marginal increase with higher CO 2 (Figs. 4 and 11e). However, when there was a vertical V cmax gradient, φ F was much higher in the lower canopy at lower CO 2 , potentially resulting in higher SIF at low atmospheric CO 2 , which was contrary to the IJKX prediction. The erroneous predicted SIF patterns of 1X and KX highlighted the importance of appropriately averaged leaf APAR, particularly the partitioning of sunlit and shaded leaves.

Dataset compatibility
Our model simulations showed that different canopy complexity levels predicted divergent carbon, water, and SIF fluxes. 1X and KX without partitioning the canopy to sunlit and shaded fractions, in particular, showed very high biases compared to the other three levels of complexity, namely 2X, 2KX, and IJKX. Further, as we expected, IJKX, which has the most complex canopy, had the lowest predicted carbon and water fluxes, followed by 2KX and 2X and then KX and 1X. Moreover, when we accounted for a profile of vertical canopy photosynthetic capacity, the difference among canopy complexity levels increased. Though 2KX and 2X were, in general, close to IJKX in predicted canopy fluxes, the disagreements may range up to > 20 % (maximum) for 2KX and up to > 40 % (maximum) for 2X. Given the differences in predicted fluxes using different canopy complexity levels and the varying difference (not a constant ratio), we do not recommend using photosynthetic parameters inverted from different canopy complexity models; i.e., parameter fitting has to be performed with the same underlying model as for the full forward modeling. Given the higher realism of the enhanced complexity models, however, leaf-level fits of photosynthetic parameters could be employed in models of higher complexity but would result in high biases when used in simple big-leaf models.
The disagreements among canopy complexity levels make it difficult to parameterize a land model using a complex canopy setup and thus hamper the fusion of large-scale remote-sensing-based datasets with land models at a global scale. Thus, it is necessary to revisit the flux and plant trait inversions using more applicable land model setups to make sure the inverted datasets and land models are consistent in their assumptions. This is the only way to ensure that inverted parameters are quantitatively useful in future land surface modeling. Moreover, it is also possible for land models to go without the inverted fluxes or traits if the land model runs using a complex canopy such as IJKX. This way, the model can be directly compared against satellite observations (Shiklomanov et al., 2021) without an intermediate step that performs the inversion from radiation observation canopy properties and thus surface water and carbon fluxes.

Necessity of a complex canopy
As suggested by Bonan et al. (2021), modelers need to move to a multi-layered canopy modeling to account for the vertical profiles and microclimates in the canopy. Further, to better utilize the broadly available remote sensing data, modelers need to move from broadband radiation to hyperspectral radiation and from sun/shade fractions to leaf angular distribution. One may ask whether it is necessary to implement a way more complex and inefficient multi-layered canopy with leaf angular distributions to account for an average 5 %-22 % difference, while the difference can be compensated by tuning plant traits such as photosynthetic capacity and hydraulic conductance. The answer varies depending on what types of data are used in the model. If one uses parameters meant to use with 2X (namely a two-leaf canopy), using a multilayered canopy such as 2KX and IJKX would not improve the model performance but instead could result in higher biases. In this case, we suggest keeping the same canopy complexity as used to derive plant traits. However, if one wants to bridge plant physiology to both leaf-level measurements as well as remotely sensed data such as the reflection and fluorescence spectra, we would suggest using IJKX or using an even more complicated canopy model to be as accurate as possible. We note here that 2KX approximates IJKX well with an average 3 %-12 % difference, and 2KX would be useful to speed the calculations for more qualitatively oriented research as the trends generally agree between 2KX and IJKX.
We recognize that increasing model complexity can make it (a) less user-friendly for researchers to use (e.g., when implementing the model into their research projects) and (b) slower to run the model, particularly using less efficient programming languages such as Python and R (compared Figure 11. Leaf fluorescence responses to radiation, CO 2 partial pressure, and the leaf maximum carboxylation rate. (a) Photosynthesis system II quantum yield responses to leaf absorbed photosynthetically active radiation (APAR) and leaf internal CO 2 partial pressure. (b) Leaf fluorescence quantum yield (φ F ) responses to APAR and leaf internal CO 2 . (c) Product of φ F and APAR vs. APAR for leaves with different internal CO 2 (number labeled next to each curve; unit: Pa). (d) Product of φ F and APAR vs. internal CO 2 at different APAR values (number labeled next to each curve; unit µmol m −2 s −1 ). The simulations of panels (a-d) are done at a leaf temperature of 25 • C and a maximum carboxylation rate (V cmax ) of 60 µmol CO 2 m −2 s −1 . (e) φ F of different canopy layers was at a given atmospheric CO 2 partial pressure. V cmax was the same among canopy layers. The simulation results are from Fig. 4e. (f) φ F of different canopy layers was at a given atmospheric CO 2 . V cmax was lower in the lower canopy. The simulation results are from Fig. 5e. to C). In our highly modularized CliMA Land model, we use Julia, a just-in-time compiled programming language that allows for the versatility of a scripting language like Python but with the speed of fully compiled languages such as C and Fortran (see https://julialang.org/benchmarks/, last accedss: 22 December 2021). The CliMA Land model can simulate canopy radiation using either the mSCOPE-based radiative transfer scheme (Yang et al., 2017) or the traditional sunlitand shaded-fraction scheme (e.g., De Pury and Farquhar, 1997;Campbell and Norman, 1998). Further, CliMA Land supports both stomatal optimization models (including those from Sperry et al., 2017;Anderegg et al., 2018;Eller et al., 2018;Wang et al., 2020) and empirical stomatal models (including those from Ball et al., 1987;Leuning, 1995;Medlyn et al., 2011). For the empirical stomatal models, CliMA Land supports using an ad hoc tuning factor to account for stomatal responses to soil moisture through tuning either the empirical fitting parameter (such as g 1 in the Ball et al., 1987, andMedlyn et al., 2011, models) or the leaf photosynthetic capacity (as done in Kennedy et al., 2019). Users may freely customize the model setup by choosing among the provided alternatives. We believe the practice of making land models more open and modular will benefit the land model and plant physiology communities in future research.

Conclusions
We evaluated how much canopy carbon, water, and SIF fluxes differ when using five different canopy complexity levels in a land model. We found that when using the same model inputs, simpler canopy models predicted higher carbon, water, and SIF fluxes, and when we accounted for a profile of vertically heterogeneous photosynthetic capacity, we found more disagreements among canopy models with varying complexity levels. We also found that the modeled SIF varied with sun-sensor geometry among tested canopy complexity levels. Our model results suggest that misusing parameters inverted from different canopy complexities and assumptions may have resulted in biases in predicted canopy fluxes, and thus we recommend more cautious model parameterization regarding canopy complexity levels. Further, we recommend using complex canopy models with leaf angular distribution and a hyperspectral radiation transfer scheme to compare against remote sensing data in order to accurately mimic observed radiation. However, the use of complex canopy models in land surface modeling may be less efficient and not user-friendly for researchers. We believe more open and modular land models like CliMA Land will help lower the threshold to researchers.
Code and data availability. We coded our model and did the analysis using Julia (version 1.6.2), and the current version of the CliMA Land model is available from the project website at https: //github.com/CliMA/Land under the Apache License 2.0. The exact version of the model used to produce the results employed in this paper is archived on CaltechDATA (https://doi.org/10.22002/D1.2316, Wang and Frankenberh, 2021).
Author contributions. YW and CF designed and conducted the research, performed the general data analysis, and wrote the paper.
Competing interests. The contact author has declared that neither they nor their co-author has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.