Does dynamically modeled leaf area improve predictions of land surface water and carbon fluxes? Insights into dynamic vegetation modules

Land surface models represent exchange processes between soil and the atmosphere via the land surface by coupling water, energy and carbon fluxes. As a strong mediator between these cycles, vegetation is an important component of land surface models. Some land surface models include modules for vegetation dynamics, which allow for the adjustment of vegetation biomass, especially leaf area index, to environmental conditions. Here, we conducted a model–data comparison to investigate whether and how vegetation dynamics in the models improve the representation of vegetation processes and related surface fluxes in two specific models, ECLand and Noah-MP, in contrast to using prescribed values from lookup tables or satellite-based products. We compared model results with observations across a range of climate and vegetation types from the FLUXNET2015 dataset and the MODIS leaf area product and used on-site-measured leaf area from an additional site. Yet, switching on the dynamic vegetation did not enhance representativeness of leaf area index and net ecosystem exchange in ECLand, while it improved performance in Noah-MP only for some sites. The representation of energy fluxes and soil moisture was almost unaffected for both models. Interestingly, the performance regarding variables of the carbon and water cycles was unrelated for both models such that the weak performance of, e.g., leaf area index did not deteriorate the performance of, e.g., latent heat flux. We show that one potential reason for this could be that the implemented ecosystem processes diverge from the observations in their seasonal patterns and variability. Noah-MP includes a seasonal hysteresis in the relationship between leaf area index and gross primary production that is not found in observations. The same relationship is represented by a strong linear response in ECLand, which substantially underestimates the observed variability. For both water and carbon fluxes, the currently implemented dynamic vegetation modules in these two models did not result in better model performance compared to runs with static vegetation and prescribed leaf area climatology.

1 Introduction

Land surface models (LSMs) represent the energy, water and biogeochemical cycles at the land surface. Traditionally, their main purpose has been to provide a surface component in coupled atmosphere–land models. LSMs are applied in meteorological models, reanalysis products or the Coupled Model Intercomparison Project (CMIP). However, their scope is widening, and new fields of application, like historical land cover change simulations (Lawrence et al., 2018) or flood alert services (Harrigan et al., 2023), are arising. There is active development within the land surface modeling community, with more and more features being added to existing models to make them more realistic (Blyth et al., 2021). Given the wide use of these models and the implications of their results, extensive model validation has been done already. Model validation covers a wide range of water, energy and carbon fluxes at global, regional and site scales (e.g., Niu et al., 2011; Haverd et al., 2018; Lawrence et al., 2019; Boussetta et al., 2021). Such works that introduce individual evaluation schemes are often accompanied by studies that perform comparisons between models (e.g., Best et al., 2015; Krinner et al., 2018). Comparisons like these are conducted for different reasons. For example, one aim is to create a ranking between models that allows for assessment against alternative schemes. Using this method, Best et al. (2015) reported that simple statistical methods achieve a higher performance in energy partitioning at eddy covariance sites than any single LSM tested. One limitation of their study is that they did not report metrics of individual model performance but only normalized ones. This procedure does not allow model users to judge whether the investigated methods have achieved a (dis-)satisfactory performance, since all methods might have a poor individual model performance. Other challenges in these activities are to maintain a standard protocol for model comparison while not creating a superficial performance contest among them and to minimize human errors (Menard et al., 2021). Haughton et al. (2016) more closely explored the cause of poor model performance of LSMs shown in the PLUMBER study by Best et al. (2015), which they presented as the bias for the evaporative fraction (EF) derived from various tower sites exemplarily. From all investigated aspects they concluded that mismatches between modeled and observed heat fluxes are most likely caused by calculations within the models and not related to errors in the observations. Yet, specific reasons for this mismatch, for example, over-parameterization, missing processes or calibration issues, cannot be identified by benchmarking studies or model rankings alone but require further investigation of individual model performance. At the same time, the causes of poor model performance can be multifaceted, rendering their identification challenging (Haughton et al., 2018b). Nonetheless, further LSM development needs understanding of how individual process implementation and parameterization affect model performances. A wealth of studies evaluated different LSMs with respect to radiation, heat fluxes or surface temperature, and carbon fluxes. Carbon fluxes like gross primary production (GPP) are often validated by using global gridded fluxes like FLUXCOM (Ma et al., 2017; Jung et al., 2019; Lawrence et al., 2019). The correct implementation of ecosystem processes and related variables is crucial for using LSMs in assessing impacts due to climate change, for example in drought evaluation (Ukkola et al., 2016; Dirmeyer et al., 2021), because plant transpiration directly links the terrestrial carbon and water cycle. For example, a substantial underestimation of evapotranspiration by eight LSMs during drought conditions was shown across different plant communities (Ukkola et al., 2016). De Kauwe et al. (2015) concluded from their simulations of drought responses for the European FLUXNET sites with the Community Atmosphere Biosphere Land Exchange (CABLE) model that accounting for differing drought sensitivity of plant communities in LSMs may be required to correctly capture drought impacts. Currently, most LSMs are not able to represent direct vegetation control on surface exchange, in part because they underrepresent biophysical responses to changing water availability and oversimplify vegetation dynamics, in particular leaf area index (LAI) (Forzieri et al., 2020). LSMs typically work with climatological LAI, e.g., seasonality read from lookup table files, or calculate LAI as a prognostic variable internally. At the same time, LAI has a large impact on both water and carbon fluxes (e.g., Fisher et al., 2014), and an understanding of how its parameterization impacts flux estimates by LSMs would help to shed light on the known discrepancies in representing vegetation. Here, we investigate model performance for water and carbon fluxes with a focus on vegetation processes. We additionally check the reasons for model–data mismatch by analysis of the underlying computer source code of the models (as stated by Dirmeyer et al., 2018), which can only be done for a limited set of models due to the large effort that is needed. For this scope, we chose ECLand and Noah-MP as frequently used and continuously developing LSMs with available modules for vegetation dynamics. In this paper, we aim to answer the following research questions. (1) Does the representation of net ecosystem exchange (NEE) and LAI improve if ECLand or Noah-MP represents vegetation dynamically? (2) How does dynamic vegetation in ECLand or Noah-MP impact other variables like heat fluxes and soil moisture? Do improvements in model performance for one variable compromise performance for other variables? (3) What are the mechanics behind modeled temporal patterns in vegetation dynamics and occurring misfits to the observations?

2 Methods

2.1 Data basis

2.1.1 Site selection

The FLUXNET2015 dataset (Pastorello et al., 2020) provides measurements from globally distributed eddy covariance sites. We selected a subset from all available FLUXNET sites, focusing on sites with long observation periods, covering different vegetation types and a gradient in aridity within each vegetation type. Vegetation types within FLUXNET rely on the International Geosphere–Biosphere Programme (IGBP) land classification (National Center for Atmospheric Research, 2022). The aridity index of all sites was retrieved from the CGIAR-CSI Global-Aridity and Global-PET Database (Trabucco and Zomer, 2018) and inverted afterwards, bringing it back to the definition as the ratio of the long-term mean annual potential evapotranspiration to the long-term mean annual precipitation by Budyko (1974). We excluded sites with observation periods of less than 6 years because they might not represent the local climate (Haughton et al., 2018a), and extreme years could create a systematic bias. Due to the small number of sites per vegetation type with long observation periods, the vegetation types savanna (SAV), woody savanna (WSA) and open shrubland (OSH) were merged into one savanna group before continuing with the selection procedure. For each vegetation type or group, first, we chose the site with the longest observation record. Next, other sites with similar aridity (±0.1 logarithmic aridity index) were dropped to avoid an overrepresentation of some vegetation type–aridity combinations due to heterogeneous site distribution within FLUXNET. We used logarithmic values to create a linear scale of the aridity index so that a selection of too many dry sites was prevented. Afterwards, we repeated these steps for the remaining sites and continued until no more sites were available for selection in this vegetation type or group. For the selected sites, we double-checked data availability and quality and replaced a site with an alternative site if necessary. The most common reasons for discarding sites were missing or poor-quality soil moisture data or low-quality gap filling, which reduced the length of the observation record below the threshold of 6 years. By doing so, only two sites with mixed forests (MFs) were left, which is critically low. Thus, we included all MF sites into the deciduous broadleaf forest (DBF) vegetation type and repeated the selection for this group. We were left with 24 sites, covering a wide range of site characteristics as recommended by Haughton et al. (2018a), including aridity, vegetation types and observation periods (Fig. 1). Additionally, we also used data of the eddy covariance site Hohes Holz (Rebmann and Pohl, 2023), which is part of the TERENO Harz/Central German Lowland Observatory (Wollschläger et al., 2016) and has been included in the ICOS network since 2019, because on-site-measured LAI data were available for that DBF site.

Figure 1 Selected FLUXNET sites grouped by their vegetation type. For each group, sites were chosen to cover a gradient in aridity (y axis) if available. The vegetation types are grassland (GRA), savanna (SAV), woody savanna (WSA), evergreen broadleaf forest (EBF), cropland (CRO), mixed forest (MF), deciduous broadleaf forest (DBF) and evergreen needleleaf forest (ENF). The color scale represents the duration of the available time series in years.

2.1.2 Variables used and data pre-processing

From the FLUXNET (Pastorello et al., 2020) and Hohes Holz (Rebmann and Pohl, 2023) datasets, air temperature, downward short- and longwave radiation, wind speed, relative humidity, air pressure, and precipitation were used for model forcing. Turbulent fluxes, i.e., latent heat flux (LE) and sensible heat flux (H), as well as net ecosystem exchange (NEE), gross primary production (GPP) and volumetric soil water content in 10 cm depth, were used for model evaluation. All data were provided and used at half-hourly resolution. FLUXNET data were retrieved from their website (fluxnet.org, 2020). LE and H in FLUXNET2015 are available in two different variables: one is a product that corrects the turbulent fluxes for energy balance closure, while the other one provides a continuous time series filled by marginal distribution sampling. We decided to use the first one as long as it was available in the dataset since LSMs also consider energy balance. Missing data in the Hohes Holz meteorological dataset were filled using a Kalman filter (Sayed, 2003) for short gaps up to 3 h, except for precipitation, which was set to zero. For longer gaps, the Kalman procedure tends to overestimate the observations, which resulted in offsets at the end of the gap-filling periods. Thus, gap-filling data for these gaps were retrieved from the ERA5 (Hersbach et al., 2020) data product (Muñoz Sabater, 2019) with 0.1° spatial resolution and 1 h temporal resolution. For calculation of the evaporative fraction $\frac{\mathrm{LE}}{\mathrm{LE}+H}$, all time steps with H≤0 were excluded. The same time steps were left out for LE to focus the comparison of turbulent fluxes on periods with evaporative demand. For estimation of model performance, we excluded gap-filled periods that were longer than 1 month.

2.2 Model description

We investigated how dynamic vegetation affects model outputs in two land surface models capable of representing both static and dynamic vegetation: ECLand (Balsamo et al., 2009; Dutra et al., 2010; Boussetta et al., 2021) and Noah-MP (Chen and Dudhia, 2001; Ek et al., 2003; Niu et al., 2007, 2011).

2.2.1 ECLand

The European Centre for Medium-Range Weather Forecasts (ECMWF) developed a Carbon-Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (CHTESSEL) (Balsamo et al., 2009; Dutra et al., 2010; Boussetta et al., 2013), which represents the land component of the Integrated Forecasting System (IFS). As part of the IFS, CHTESSEL has evolved into a more flexible system ECLand (Boussetta et al., 2021), which also allows for several modular extensions. Among these, a dynamic vegetation module simulates the temporal evolution of vegetation. Therein, LAI, vegetation biomass and vegetation coverage are calculated from the daily carbon budget instead of taking them from the climatological LAI. However, LAI climatology can still be used for fully static or partly dynamic simulations. In ECLand (IFS version CY46R1), each of the 19 vegetation types receives its own parameter values (e.g., for roughness lengths, stomata resistance, root distribution) from lookup tables (Boussetta et al., 2012, 2021). These vegetation types are categorized into high or low vegetation. Each grid cell has one dominant high-vegetation type and one dominant low-vegetation type, together forming the vegetation of a grid cell (Balsamo et al., 2009). Surface fluxes are computed for the high- and low-vegetation tiles separately and are then merged for the whole grid cell according to their fractional cover. The vegetation coverage is calculated from a prescribed climatological vegetation fraction (part of the input) and a vegetation-type-dependent density (from a lookup table) and corrected by current LAI (Boussetta et al., 2021). Net assimilation results from carbon uptake of atmospheric CO₂ by the current leaf area (defines absorbed radiation) and is restricted by environmental factors such as soil moisture and nitrogen availability (important equations can be found in Appendix A). Together with the dark respiration and after scaling with a quantum use efficiency factor, potential gross assimilation is calculated. This value is then linearly linked to LAI and the humidity-corrected air density, resulting in gross primary productivity (GPP). With activated vegetation dynamics, a potential net assimilation, together with LAI, forms a damping factor for biomass senescence. Biomass senescence is determined from the current biomass, linearly linked to the current LAI, and the damping factor. The change in biomass results from this updated biomass and the net assimilation. Then, biomass is updated again and linearly transferred into updated LAI by using a specific leaf area from a lookup table (Boussetta et al., 2021). For static ECLand, the prescribed climatological LAI is used. LAI in ECLand determines the canopy resistance for water vapor transport and thus the evapotranspiration and the interception (Boussetta et al., 2012, 2013, 2021).

2.2.2 Noah-MP

Noah-MP is the widely used community Noah land surface model (Chen and Dudhia, 2001; Ek et al., 2003) with multi-parameterization options (Niu et al., 2007, 2011). Predicted LAI in Noah-MP is calculated based on leaf carbon allocation and specific leaf area per vegetation type (Ma et al., 2017). In contrast to ECLand, Noah-MP can either use prescribed LAI values per vegetation type or depend solely on dynamic LAI estimates, without the option to mix between the two. In Noah-MP (version HRLDAS 3.9), parameter values (e.g., value range of stomatal resistance, number of rooted soil layers, specific leaf area) of the 27 vegetation types are taken from lookup tables. The vegetated sub-grid area of each grid cell is dominated by one vegetation type forming a one-layer canopy. For the calculation of canopy interception and transpiration, Noah-MP considers aerodynamic and stomatal resistances for the water vapor and carbon fluxes within the canopy and between the canopy and the atmosphere (Ma et al., 2017). Among others, stomatal resistance is predominantly controlled by photosynthesis (Niu et al., 2011), which depends on leaf area, and is limited by light and root zone soil moisture (important equations can be found in Appendix A1). Assimilation depends on LAI and is constrained by physiology and light availability. Assimilated carbon is allocated to different plant tissues (leaf, stem, wood, root), forming GPP, and reduced by respiration, dying and turnover processes, such as drought stress and senescence representing leaf dynamics (Dickinson et al., 1998). The respiration rate is determined by LAI, GPP, temperature and soil moisture stress. Carbon that is allocated to leaves together with biomass losses forms an updated leaf biomass, which converts into the LAI by using a specific leaf area (Ma et al., 2017). Carbon assimilation and allocation, and thus also GPP and NEE estimation, are usually deactivated for the static Noah-MP since a prescribed LAI is given. We adapted the model code in a way that GPP and NEE for the static simulations are calculated anyways but resetting all variables that would be dynamically predicted within the same function to their prior values. This assured that the model still ran in a static configuration.

2.3 Model setup and simulations

Simulations with activated modules that predict LAI time series are referred to as activated vegetation dynamics or dynamic ECLand and dynamic Noah-MP hereafter. For both models, the reference height (level of the forcing input) was set to the flux tower height of the sites, which depends on the vegetation type. The models were set up as closely as possible to the available site information, but there are some technical differences in the structure of the model input, i.e., in the initial files. Forcing and model calculation were done in 30 min resolution if available; otherwise, hourly resolution was applied. We used four-layered soil representation, with the uppermost layer used for the evaluation of soil moisture, which is 7 and 10 cm deep for ECLand and Noah-MP, respectively. Every simulation started with a 10-year spin-up phase by recalculating the first year.

2.3.1 ECLand

We used ERA5-based (Hersbach et al., 2020) global initial data for ECLand and selected the grid cells where the flux towers are located. These initial files contain information on albedo, orography, soil type, surface roughness and monthly LAI, which is not available in the FLUXNET metadata. For the simulations that use alternative LAI forcing, monthly LAI in the initial files was replaced by the scenario-specific alternative values (see Sect. 2.3.3). We defined the vegetation on this grid cell to be either high or low vegetation (and not a mixture) depending on the site information. Forests and savannas were treated as high-vegetation types, while grasslands and croplands were allocated to low-vegetation types. The vegetation type that fits best with the FLUXNET characterization was selected (see Table 1). The coverage of that vegetation type was set to 100 %. Meteorological forcing was taken from the FLUXNET and TERENO datasets mentioned above (Sect. 2.1.2). The ECLand simulations were done with van Genuchten's soil hydrologic parameters (van Genuchten, 1980), activated sub-grid surface runoff and activated snow parameterization (see Table A1).

Table 1 Assignment of vegetation types used in ECLand and Noah-MP and referred initial LAI. The values in brackets for Noah-MP initial LAI refer to sites in the Southern Hemisphere due to shifted seasons.

2.3.2 Noah-MP

The soil type for Noah-MP was taken from a global soil grid (Hengl et al., 2014; Poggio et al., 2021; ISRIC, 2024) by selecting the grid cell including the flux tower location. Initial values for temperatures and soil moisture were taken as the FLUXNET/TERENO observations on 1 January at 00:00 local time in the first year of the simulation period. Vegetation types were chosen to match the USGS vegetation types (University Corporation for Atmospheric Research, 2023) as closely as possible, and the initial LAI values were set according to the defaults in the parameter file (see Table 1). The vegetation cover fraction was set to 100 % so that the entire grid cell represents the vegetation type of the observation site. The minimum green vegetation fraction was set to 1 % to ensure that the whole vegetation cover does not die during winter, which would hinder temperate short vegetation from growing in spring. For the simulations with alternative LAI forcing, the monthly LAI in the lookup table was replaced by the scenario-specific alternative values (see Sect. 2.3.3). The Noah-MP simulations were done with soil parameterization from lookup tables and the Ball–Berry stomatal resistance approach (Ball et al., 1987; Bonan, 1996) using a matric potential limitation. All other selected options can be found in Table A2.

2.3.3 Leaf area index data and scenarios

Monthly LAI values are part of the initial input of both models via lookup tables. These tables contain annual cycles of LAI for each vegetation type separately. This default climatology is already based on values from MODIS. For ECLand, the gridded values of LAI were disaggregated to the high- and low-vegetation type of the grid cell for the time span 2000–2008 (Boussetta et al., 2013). For alternative LAI inputs, these values in the lookup tables were replaced manually. LAI values were taken from the MOD15A2H data product from NASA's Earthdata portal (Myneni et al., 2015). One grid cell of 500 m×500 m was selected per eddy covariance tower according to the site coordinates, and LAI values with a temporal resolution of 8 d were extracted for the years 2000 to 2014. To assure reliability of the values, the MODIS15A2H data product comes with numeric quality flags. Although Fang et al. (2012) recommend using all values with quality flags less than 64, we excluded data with quality flags of 8 because many of these LAI values were extremely low during the vegetation period, which is unrealistic. Then again, due to lacking LAI values during winter or wet seasons, values with quality flags of 73 (empirically filled with clouds present), 81 (empirically filled with mixed cloudiness) and 97 (empirically filled for other reasons) were included as a trade-off between excluding as much poorly flagged data as possible and keeping roughly the same amount of data values for each month (see MODIS documentation for more details). Afterwards, we smoothed the remaining values by using a Savitzky–Golay filter (window length: 11, polyorder: 2) (as done by, e.g., Xiao et al., 2011; Huang et al., 2021) from the SciPy package (Savitzky and Golay, 1964; Luo et al., 2005) and prepared a mean annual LAI cycle for all available years with monthly resolution, further named MODIS climatology. For an additional experiment, the monthly LAI from MODIS of each year within the simulation period separately was used as input, called MODIS single year from this point on. Missing LAI values for a month were filled by the average value of the adjacent months. If LAI values for at least 2 consecutive months were not available, the LAI values from the default lookup table were used for those months. For the Hohes Holz site, on-site-measured LAI data were available from digital cover photography (DCP), which was shown to yield comparable results to established methods (Piayda et al., 2015). For each measurement date, we averaged the values from the whole plot area and, afterwards, calculated monthly means over a time span covering 2014–2019. This alternative LAI forcing is called on-site LAI hereafter. The nomenclature of all LAI scenarios can be found in Table 2.

Table 2 Nomenclature of all model scenarios using LAI data sources.

The MODIS LAI was also applied for model evaluation but in a high temporal resolution of 8 d. Due to the usage of single-day values, we solely used data with quality flags 0 (no issues) and 32 (saturated) to lower the uncertainty. Additionally, we refrained from smoothing to avoid an offset of the LAI values and left gaps as they were. For the static runs, comparison with MODIS LAI on a daily basis provides information on how well a LAI climatology represents the local LAI evolution and whether the incorporation of more site-specific climatology can improve the representativeness. For the dynamic simulations, comparing modeled LAI with daily MODIS values is used to examine whether the models are able to capture inter- and intra-annual LAI dynamics.

2.4 Performance evaluation

Model outputs and observational data from the flux towers were averaged/summed to daily values for direct comparison. For LAI, we calculated the 8 d mean of the LAI model output to correspond to the temporal resolution of the MODIS LAI estimates. As performance criteria we used Pearson's correlation coefficient, the normalized standard deviation and a modified relative bias for the model–observation relationship. Pearson's correlation coefficient R describes the fit between model and observation values (Benesty et al., 2009) and is calculated from the NumPy package. The normalized standard deviation s_n is the ratio of the standard deviation of the model predictions to the standard deviation of the observations. It is used to describe the models' ability to reproduce the variability of the observations. The relative bias b applied here was adapted to the domain of the variable to avoid division by zero or by values very close to zero (especially important for NEE). For this purpose, the distribution of the observed values was shifted by their minimum, resulting in only positive values with a minimum of zero:

$$\begin{array}{}\text{(1)}& b={\displaystyle \frac{\stackrel{\mathrm{‾}}{y-x}}{\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{ˇ}}{x}}},\end{array}$$

whereby y represents the model predictions, x the observations, $\stackrel{\mathrm{‾}}{x}$ the mean and $\stackrel{\mathrm{ˇ}}{x}$ the minimum of the observed values. To compare the model performance between simulations with static and dynamic vegetation, we determined the change in relative bias as follows:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }b=\left|b_{\mathrm{static}}\right|-\left|b_{\mathrm{dynamic}}\right|.\end{array}$$

Negative values mean that the relative bias of the dynamic simulation was greater than that of the static simulation and thus that the performance was reduced by activating vegetation dynamics. To investigate the sensitivity of dynamically modeled vegetation to the model performance, we checked how strongly the quality of the model simulation of one target variable (e.g., LE) depends on the model quality of another (e.g., LAI). For this, we used elasticity as a metric. Elasticity is calculated as the ratio of the change in one statistical measure (analogous to Eq. 2) for two different target variables:

$$\begin{array}{}\text{(3)}& E={\displaystyle \frac{\mathrm{\Delta }{m}_i}{\mathrm{\Delta }{m}_j}},\end{array}$$

where m is one of the statistical measures mentioned above, i.e., R, s_n or b, while i and j denote different target variables, e.g., GPP or LE. For variables that are strongly related, like LAI and GPP, we expect elasticity to be positive. Two variables are considered independent if $-\mathrm{0.1}\le E\le \mathrm{0.1}$ because the change in m_j would then need to be larger than 1 order of magnitude to cause a change in m_i. Changes in model performances of the target variables were plotted in Taylor diagrams (Copin, 2021).

3 Results

3.1 Effect of dynamic or prescribed leaf area index on leaf area and carbon uptake prediction

Figure 2 shows the quality metrics for the model performance regarding LAI in a Taylor diagram. The location an optimal model simulation would occupy is indicated with a star. The model performance of the dynamic run is shown with the symbols, while the static runs can be read from the start of each arrow. The direction and length of each arrow highlight the difference in the performance metrics between static and dynamic runs. Shown are simulations started (dynamic) or run (static) with default vs. MODIS climatology. While in the Noah-MP simulations with static vegetation the model performance depended on the LAI forcing applied, the simulation results were unaffected by the type of LAI forcing with vegetation dynamics switched on since the symbols in Fig. 2c and d have the same positions. For ECLand, this was also the case for many sites but not all, e.g., AT-Neu and AU-How (Fig. 2a, b). Initializing ECLand with default climatology (Fig. 2a) and activating vegetation dynamics generally increased the variance of simulated LAI compared to static simulations but also decreased model performance; e.g., the mean Pearson correlation decreased from 0.72 to 0.62. At the same time, whether the predicted LAI fit better to MODIS observations than the default climatology was ambiguous, as can be seen by the shift in relative bias, which ranged between −0.5 and 1.3. On the contrary, the results for Noah-MP showed a different pattern (Fig. 2c) because there was no clear shift to higher variances or worse correlation when activating vegetation dynamics. Short- (GRA+CRO) and sparse (SAV+WSA)-vegetation types in particular had the highest changes towards decreased performance but also enhanced model performance for LAI. For other sites (mostly forests), modeled dynamic LAI correlated well with the observations. For both models, using MODIS climatology instead of default climatology in static simulations resulted in the best performances with regard to the LAI of all simulations (start of the arrows in Fig. 2b, d); e.g., the mean correlation coefficient increased to 0.83 and 0.84, and the mean relative bias (Table S1 in the Supplement) improved to −16 % and −2 % for ECLand and Noah-MP, respectively. This can be expected because MODIS was also used as a reference dataset for LAI evaluation. With activated vegetation dynamics, the performance of both models decreased, as all quality metrics shifted away from the point indicating the best performance on the Taylor diagram (Fig. 2b, d). The same applied to the relative biases of LAI since their shift was predominantly negative. In other words, switching on vegetation dynamics did not result in improved LAI representation compared to just using MODIS climatology.

Figure 2 Change in model quality metrics for LAI when switching on vegetation dynamics for all included sites and when using default climatology (a, c) or MODIS climatology (b, d). The star (“Observ”) marks the location of perfect correlation between the observation and model and perfect agreement between the observed and modeled variance. The model performance of the static runs can be read from the start of each arrow. When no arrow appears, either no correlation could be calculated (e.g., for evergreen forests where default climatological LAI is constant) or values could not be placed on the logarithmic axis. The symbol colors indicate the site aridity (top-right legend) as follows: very humid – aridity index (AI)<0.6, humid – AI<1.25, sub-humid – AI<1.54, dry sub-humid – AI<2, semi-arid – AI<5 and arid – AI≥5 (Ashaolu and Iroye, 2018). Vegetation types are symbolized by different marker types (bottom-right legend).

Forest ecosystems, in general, were better represented by model predictions with vegetation dynamics than short or sparse vegetation. Figure 3 shows the results of the forest site Hohes Holz in more detail. Although the representation of LAI variability deteriorated when simulating dynamic vegetation with Noah-MP, those runs resulted in LAI predictions that closely match MODIS observations (Fig. 3d–f), represented by a relative bias of −18 % and a correlation coefficient of 0.78. ECLand more generally suffered from larger relative biases in LAI, especially when simulating with vegetation dynamics (−30 % on average, Fig. 3c). The only scenario where model performance generally increased for ECLand was through switching on vegetation dynamics compared to static runs with default climatology.

Figure 3 Statistical measures for the variables LAI, NEE and GPP of the model runs for the Hohes Holz site. The categories on the y axis mark the different LAI forcings. The statistical measures of the static and dynamic simulations of the same variable are connected by a horizontal line. The vertical dotted red line marks the optimum of each measure.

In contrast to LAI, the model performance of ecosystem exchange variables in ECLand was less affected by activating vegetation dynamics. A common feature is that the variance predominantly increased when using dynamic vegetation (Fig. 4a, b). Sites with short or sparse vegetation mostly reacted more sensitively to dynamic vegetation modeling in their NEE and GPP representation, especially when forcing with MODIS climatology, which is indicated by the longer arrows in Fig. 4a and b (for GPP, see Fig. S1 in the Supplement). For forest ecosystems in general, the changes in the model performance of NEE and GPP were small, as also shown for the Hohes Holz site (Fig. 3a–c). Nevertheless, the performance of NEE (and GPP) decreased when activating vegetation dynamics, mainly driven by lowered correlation coefficients, on average from 0.41 to 0.37 (0.72 to 0.68). Only three sites showed improvements in NEE representation when predicting with dynamic ECLand and just one site for GPP. Relative bias changed in both directions, towards lower and higher model performance. Dynamic ECLand mainly overestimated NEE by 11 % on average, indicating that ecosystems were predicted to be a smaller carbon sink than observed (Table S2). Instead, dynamic Noah-MP estimated on average 10 % lower NEE compared to the observations for most sites (Figs. 4c, d, 3c, f).

Figure 4 The same as the Taylor diagram before but with NEE evaluation.

Activating the dynamic vegetation affected the model performance of NEE and GPP for Noah-MP heterogeneously. Some sites showed very small changes (e.g., IT-Lav and IT-Ren, Fig 4c, d), while the model performance of NEE was largely impacted by vegetation dynamics for other sites (e.g., US-Var and AU-DaS). In contrast to ECLand, no evidence could be found that certain vegetation types or aridity classes were more sensitive to activated vegetation dynamics in Noah-MP, and even forests showed larger changes in model performance (Fig. 3e). For GPP on the other hand, the variance predominantly increased by activating vegetation dynamics in Noah-MP (Fig. S1c, d), and the normalized standard deviation of NEE changed in opposing directions, which is another difference compared to ECLand. Despite the higher sensitivity of the NEE model performance to Noah-MP vegetation dynamics, the overall model performance was barely affected since relative bias shifted from −12 % to −11 % and the correlation coefficient from 0.50 to 0.53 on average. Changes in statistical measures can be in opposing directions, as can be seen for the normalized standard deviation and relative bias of the forest site Hohes Holz (Fig. 3e, f), which eliminated trends towards improved or reduced model performance. Only the AU-Stp site clearly improved regarding NEE representation by activating vegetation dynamics in Noah-MP by initializing with either default or MODIS climatology. The GPP performance showed small improvements by activating vegetation dynamics in Noah-MP as the mean correlation coefficient shifted from 0.68 to 0.74, and the range of the relative bias was lowered from between −32 % and +69 % to between −28 % and +42 %. In general, Noah-MP seemed to capture NEE representations better, as the mean deviance from a normalized standard deviation of 1 was 0.33 (ECLand: 0.39) and showed a higher correlation coefficient on average of 0.51 compared to ECLand (Fig. 4c). Remarkably, the four and nine best sites regarding NEE correlation and variance were forests for ECLand and Noah-MP, respectively. At the same time, all evergreen broadleaf forests suffered from poor performance in both models. GPP representation in both models was better than for NEE (Fig. S1 and Table S3). Overall, static and dynamic Noah-MP performed well in representing NEE and GPP for most forest sites apart from the evergreen broadleaf forests. In line with the finding that model performances of dynamic Noah-MP were independent of the prescribed LAI forcing, the availability of on-site LAI data for the Hohes Holz site yielded no improvement in the representation of NEE or GPP compared to other LAI climatologies (Fig. 3). The same was the case for dynamic ECLand. Forcing static ECLand with on-site LAI data resulted in a NEE and GPP correlation and relative bias that are comparable to the forcing with MODIS climatology, with only variability being lower.

3.2 On the sensitivity of heat fluxes and soil moisture to vegetation dynamics in LSMs

For both models, activating vegetation dynamics had a small impact on the representation of turbulent fluxes and soil moisture. The strongest changes occurred for short- or sparse-vegetation types or for drier climates, which had the largest arrows in the Taylor diagrams (Figs. 5, 6). In ECLand, activating vegetation dynamics enhanced the variance of the latent heat flux for most sites (from 0.80 and 0.84 to 0.94 on average for default and MODIS climatology, respectively), but correlation between the simulated and observed values remained unaffected or even diminished (mean change smaller than −0.03). For several sites, LE estimates from dynamic ECLand better represented the observations, as shown by the positive shift in relative bias (reduction from −32 % to −21 %) (Fig. 5a, Table S4), but no relationship regarding vegetation type or site aridity can be seen, and changes are generally small. A major exception appeared for CH-Oe2, which was caused by its default LAI climatology that did not fit the vegetation type. Activating vegetation dynamics in Noah-MP hardly affected model performance of LE (the mean changes in correlation, standard deviation and relative bias were 0.02, 0.00 and 0.02, respectively). Sites that showed some sensitivity predominantly have drier climate (e.g., AU-Stp, US-Var; see Fig. 5c). Several sites showed less bias in LE predictions when using dynamic vegetation predictions in Noah-MP. When using MODIS climatology as the LAI forcing, activating vegetation dynamics could be advantageous for some sites regarding LE representation (AU-Stp, CH-Fru, US-GLE), but it would mostly not lead to higher model performance.

Figure 5 Change in statistical measures for LE modeling when switching on dynamic vegetation for all included sites and by using default climatology (a, c) or MODIS climatology (b, d) as forcing.

Model performance regarding the evaporative fraction (EF) was lower compared to LE as the points are further away from the point of optimal model performance (Fig. 6). Running ECLand with activated vegetation dynamics lowered the representation of the evaporative fraction, which is demonstrated by many points in the Taylor diagram drifting away from the star indicating best performance. Thereby, the mean standard deviation changed from 0.95 to 1.08, and the correlation coefficient was reduced slightly from 0.48 to 0.46 on average (Fig. 6a, b). Exceptions were BE-Lon, US-SRM and US-Ton, where model performance improved slightly regarding correlation and variability. Again, the relative bias of EF changed in both directions without any trend regarding vegetation type or aridity for both models (see also Table S5). For Noah-MP, eight sites showed an improved representation of the evaporative fraction when running the model with vegetation dynamics. This number was reduced to six when the model was initialized with MODIS climatology. But changes were very small on average.

Figure 6 The same as before but for evaporative fraction, which represents the turbulent flux partitioning.

Regarding soil moisture, the model performance was almost insensitive to the vegetation dynamic option used or the type of LAI forcing for both models (Fig. S2). Some sites showed improvement in soil moisture prediction by activating vegetation dynamics for both models, although the improvement was very weak. Interestingly, no humid site was among them. However, the simulation of soil moisture resulted in a broad range of model performances, starting from very well fitting predictions (R>0.9, b≈0 %) up to very poor fitting predictions (R<0.2, $b<-\mathrm{40}\mathit{\%}$ or b>100 %; see Table S6). To investigate the sensitivity of dynamically modeled vegetation to the model performance, we checked how strongly the quality metrics of NEE, GPP, LE and soil moisture change with the quality metrics of LAI. For this, we used elasticity (defined in Eq. 3) as a metric, which is summarized for all sites in Fig. 7. Surprisingly, the quality metrics of these closely related variables were independent of each other; i.e., the elasticity was very low (within the grey band) or randomly distributed around zero. The strongest connection of all pairs tested could be found for the correlation coefficient between LAI and GPP in ECLand when using MODIS climatology but without affecting the normalized standard deviation or relative bias. Here, the mean elasticity of the correlation and normalized standard deviation is positive, meaning that, as expected, an increased model performance in LAI co-occurs with enhanced performance for GPP of the same order of magnitude. In a similar manner, NEE and LAI performances were positively related regarding the correlation coefficient in Noah-MP. Other elasticity values that include LAI were predominantly small. In other words, changes in the model quality for LAI, for most of the sites, do not affect the model performance of LE or soil moisture, and even the model performance of carbon fluxes remains unaffected.

Figure 7 Density plots showing the elasticity of the correlation (left column), normalized standard deviation (middle column) and relative bias (right column) for different variable relationships in both models when activating dynamic vegetation and using default climatology (blue) or MODIS climatology (orange) as forcing. For reasons of practicability, elasticity is used reciprocally. Accordingly, the explanatory variable is the first one of each relationship shown on the y axis. The grey-shaded area marks the range between the thresholds of independence.

3.3 Observed and simulated relationships between ecohydrological variables

One possible explanation for the small contribution of the model quality of LAI to that of the turbulent fluxes could be a weak relation between LAI and carbon exchange in the model. However, this is not the case, as illustrated in Fig. 8. The figure shows the relation between GPP and LAI for four exemplary sites: DE-HoH is a deciduous broadleaf forest in a humid climate, IT-Ren is an evergreen needleleaf forest in a semi-arid climate, GF-Guy is an evergreen broadleaf forest in a tropical climate and US-Var is a grassland in a semi-arid climate. In general, the relationships between GPP and LAI are much more scattered in the observations (top row) compared to the models (other rows), and this is true for both models, across biomes and vegetation types.

The two European sites (left columns, De-HoH and IT-Ren) reach maximum LAI and GPP in JJA and minimum values in DJF, leading to a correlation that is mainly governed by the seasonal cycle. Similarly, at the US site, with an overall closer relation, vegetation productivity and LAI peak together in spring (i.e., MAM). For these three sites, correlation coefficients range between 0.80 and 0.86, indicating a clear but not perfect relation between LAI and GPP. However, the scatter of the observed relation is considerable, with the standard deviation of the residuals (σ_r) being between 58 and $\mathrm{102}\times {\mathrm{10}}^{-\mathrm{6}}$ gCO₂ m⁻² s⁻¹. The variance is the highest for the peak of the growing season, when GPP quickly responds to environmental conditions (e.g., cloudiness, precipitation and soil moisture stress) that LAI responds much slower to. The tropical site in French Guiana (GF-Guy) shows, as expected, no seasonal cycle, leading to an extremely weak relation between LAI and GPP. The latter is comparatively high all year round (GPP between 250 and $\mathrm{600}\times {\mathrm{10}}^{-\mathrm{6}}$ gCO₂ m⁻² s⁻¹), although LAI values from the MODIS dataset surprisingly varied between 1 and 7 m² m⁻². For this tropical site, GPP and LAI dynamics seem decoupled (Fig. 8c).

Noah-MP shows a non-linear relationship with a pronounced hysteresis effect at all sites except the tropical one (Fig. 8m–p). Thereby, GPP increases linearly with LAI during biomass buildup to a point where allocation to leaves becomes minimal (vegetation type specific), and it drops considerably without any substantial reductions in LAI towards the end of the growing season (e.g., Fig. 8m). When GPP values reduce below approximately $\mathrm{100}\times {\mathrm{10}}^{-\mathrm{6}}$ gCO₂ m⁻² s⁻¹, then LAI reduces from values of about three to values of zero. This hysteresis is shifted in seasons due to local climate, as for the US-Var site (Fig. 8p). At the tropical site, Noah-MP shows some variability in GPP but almost no change in LAI, which is around a value of five.

ECLand, in general, shows a linear relationship with considerably less variability compared to the observations. The slope and intercept of the linear regression are dependent on the choice of static or dynamic vegetation. Dynamic ECLand shows a very close linear relation between LAI and GPP with much lower scatter compared to the observations (Fig. 8 third row), as R is larger than 0.99, and σ_r is between 10 and $\mathrm{14}\times {\mathrm{10}}^{-\mathrm{6}}$ gCO₂ m⁻² s⁻¹ for all non-tropical sites. With slope values of 104 to $\mathrm{254}\times {\mathrm{10}}^{-\mathrm{6}}$ gCO₂ m⁻² s⁻¹, this relationship is much steeper than in the observations. Even for the tropical site, the relationship between LAI and GPP is clearly and closely linear (Fig. 8k).

Figure 8 Scatter plots of the relationship between LAI on the x axes and GPP on the y axes as 8 d averages for four selected sites (columns). The rows from the top to the bottom show the observations, static ECLand model output, dynamic ECLand model output and dynamic Noah-MP model output. Seasons are represented by different dot colors. The arrows represent the range of GPP and LAI values for the individual seasons. A simple linear regression model was applied (dashed blue line), and its correlation coefficient (R), slope and standard deviation of the residuals (σ_r) are given for each relationship.

4 Discussion

4.1 Using LAI climatology for ECLand and Noah-MP runs is the best way to reproduce leaf area and carbon uptake

Comparison between model output and observational data of LAI, NEE or GPP on a daily basis is rarely done. The ability of the two models to reproduce these observed ecosystem variables was in line with previous results. For Noah-MP, model quality metrics were in the range of other studies (Brunsell et al., 2020; Li et al., 2022; Xu et al., 2021; Liang et al., 2020). However, dynamic LAI modeled by Noah-MP in our assessment with a mean of +70 % was more biased compared to a mean of +20 % for annual LAI values reported by Li et al. (2022). Ma et al. (2017) reported a relative bias in GPP of 40 % on average, which is higher than the relative bias found here. For ECLand, we could not find any comparable study reporting the performance of daily LAI or NEE/GPP specifically, for both dynamic and static simulations. However, for static ECLand, correlation coefficients between the modeled and observed NEE and GPP were in line with those obtained by Boussetta et al. (2013) for 10 d averages at several FLUXNET sites. Moreover, for the mean annual cycles of NEE and GPP, Stevens et al. (2020) found a lower prediction error (RMSD) when using MODIS LAI forcing compared to default prescribed LAI and, like in our study, a substantial bias in LAI.

For both models, using MODIS climatology in static simulations resulted in the best performances concerning LAI. This agreed with expectations. Since all our simulations were validated with MODIS data, using static runs with MODIS climatology itself would likely yield better results than the default values in either model. For Noah-MP, static simulations with MODIS climatology indeed yielded the best performance regarding LAI in some sites, but, interestingly, using the default climatology also performed well for others. LAI deviations with the default climatology occurred specifically in short vegetation, which was also true for the dynamic runs (see below). For ECLand, where the default climatology is already based on MODIS data (Boussetta et al., 2012), the performance of the static run was generally improved compared to the validation dataset, as the higher spatial resolution allowed for a better geographical mapping. Additionally, ECLand default climatology was created by disaggregating the total LAI in the MODIS data to the low- and high-vegetation type on the grid cell. Both points together can explain the better performance for LAI of static ECLand simulations with MODIS climatology compared to default climatology.

Dynamic vegetation did not yield better LAI results compared to using static runs with MODIS climatology for either model. Evergreen broadleaf forests showed the lowest correlation coefficients for dynamic LAI predictions, which was also shown by Yang et al. (2011) for tropical regions simulated by Noah-MP. Additionally, Brunsell et al. (2020) reported overestimation of LAI with dynamic Noah-MP for the eastern Amazonian forest, which we could not find here. ECLand suffered from overall strong relative biases regarding LAI in dynamic simulations. The underestimation of prognostic LAI has already been shown by Boussetta et al. (2021). Substantial biases also occurred in Noah-MP (also shown by  Huang et al., 2022), but particularly for short- or sparse-vegetation types. The latter could be due to LAI overestimation in the early growing season as reported by Cai et al. (2014). Moreover, Liu et al. (2016) found that neither the lookup table LAI nor the predicted dynamic LAI annual cycles seemed to reproduce LAI observations for short vegetation. On the other hand, Pilotto et al. (2015) achieved satisfactory model predictions for crop sites without vegetation dynamics. Thus, for short vegetation such as grasslands, the Noah-MP crop module perhaps better represents LAI dynamics (Liu et al., 2016), which should be tested in the future.

The performance of NEE and GPP in ECLand was not very sensitive to different vegetation dynamics. Generally, using static MODIS climatology yielded the best predictions of GPP and NEE, although the correlation between the modeled and observed NEE was generally low (the mean Pearson correlation coefficient was 0.44). In many sites, even static simulations with default climatology resulted in comparable performances. Interestingly, adding more detailed information by using MODIS single-year LAI forcing did not further improve model performance (not shown), as we would have expected if LAI dynamics contributed substantially to enhancing model performance for the carbon fluxes. However, other authors found improved model performance of turbulent fluxes, GPP and soil moisture for roughly 50 % of their set of sites by updating the LAI forcing using near-real-time data assimilation (Boussetta et al., 2015). In other words, a more frequent reset of LAI to the correct value can improve the ECLand performance in general but did not have an effect on the annual resolution applied here.

Assimilation of LAI during model runs instead of fixed forcing (as in a static case) also improved LAI and GPP model quality in a study by Xu et al. (2021) using dynamic Noah-MP. We therefore expected that LAI dynamics potentially improve model quality regarding carbon fluxes, which was predominantly not the case. However, dynamic Noah-MP is already known to overestimate GPP (Ma et al., 2017; Liang et al., 2020; Brunsell et al., 2020). Short- and sparse-vegetation types in particular suffered from low predictive power, mainly in NEE correlation (Yang et al., 2021) and in GPP relative bias (Li et al., 2022), which could also be observed here. None of the parameter sets tested by Yang et al. (2021) for simulations with dynamic Noah-MP resulted in well-fitting predictions of daily changes in NEE for three of the four sites with short vegetation within ChinaFLUX. Note, however, that LAI, NEE and GPP of short and sparse vegetation were also not well-represented in static runs. Moreover, Kumar et al. (2019) could only achieve marginal improvements in GPP representation with dynamic Noah-MP due to LAI assimilation for crops and grasslands, which suggests that LAI dynamics had only a limited effect on simulated NEE there. One possible reason for the lower predictive power of the models regarding the carbon fluxes of short or sparse vegetation could be that these vegetation types more quickly and dynamically respond to fluctuations in the environment (e.g., soil moisture limitations). Forest ecosystems might be able to compensate for restrictions through a larger intrinsic carbon storage or deeper roots, resulting in less variability of their productivity within and between years. It seems that both models investigated here cannot mimic this differentiation. Nonetheless, it was shown here that correlation coefficients for GPP simulated with dynamic Noah-MP were high (also found by  Liang et al., 2020; Li et al., 2022), and, at the same time, relative bias was small for all forests except the evergreen broadleaf forests (see Sect. 3.1 and Fig. S1). Thus, although some previous studies found substantial uncertainties in modeled GPP for different vegetation types (Ma et al., 2017; Liang et al., 2020; Li et al., 2022), predicting ecosystem variables using dynamic Noah-MP could be useful, at least for forests in studies where LAI climatology cannot be used, such as climate change impact studies.

Considering the opposing biases in NEE (and GPP) indicates that the models differ in their estimates of ecosystem respiration. One important difference is the sequence of the calculation of GPP, NEE and respiration. ECLand first estimates net assimilation and respiration separately, whereby respiration is set to be 11 % of net assimilation, and then both are used to calculate GPP. In Noah-MP, the first estimate is for GPP, which is reduced by respiration to gain a value for NEE, and, additionally, respiration is scaled by GPP and available biomass inclusive of LAI. Including our findings, for dynamic ECLand, the underestimated LAI directly transfers into lower NEE values and thus also to GPP since respiration is a fixed fraction of NEE. Apart from the fraction of GPP that is directly needed for metabolism, the estimation of respiration in dynamic Noah-MP also considers leaf maintenance, which is another difference compared to ECLand. As a result, respiration is slightly overestimated in ECLand and slightly underestimated in Noah-MP.

4.2 Model performance with respect to turbulent fluxes and soil moisture was almost unaffected by vegetation dynamics in both LSMs

The model performance of ECLand and Noah-MP regarding heat fluxes and soil moisture seems almost insensitive to vegetation dynamics.

Correlation, variability and bias of turbulent fluxes in this study were comparable to other studies. While evaluating static ECLand with FLUXNET data, Stevens et al. (2020) found correlation coefficients of 0.79 and 0.77 for the annual cycle of latent and sensible heat, respectively, and Boussetta et al. (2013) showed a mean correlation coefficient of 0.81 for 10 d averages of latent heat. For Noah-MP, statistical measures for turbulent fluxes and soil moisture were mostly in line with other studies (Niu et al., 2011; Ma et al., 2017; Yang et al., 2018; Xu et al., 2021), although Pilotto et al. (2015) presented lower correlation coefficients between 0.20 and 0.43. Interestingly, Ma et al. (2017) showed opposing relative bias for evapotranspiration on an annual timescale over the continental US of 4 % and 22 % for static and dynamic simulations, respectively.

For ECLand, whether vegetation was simulated dynamically instead of statically had little impact on turbulent fluxes. Model performance for LE and EF changed only for some sites and towards lower performance (see Sect. 3.2). The predominant underestimation of LE agrees with the findings of Stevens et al. (2020). For dynamic ECLand, the underestimation of GPP and LAI (also in Boussetta et al., 2021) could also be the reason for the poor correlation of EF between modeled and observed values because the energy fraction that is used for transpiration is underestimated. Boussetta et al. (2021) found that dynamic vegetation in ECLand improved numerical weather predictions. Their main improvements in model performance were achieved through updating land cover maps and the LAI in the lookup table or by including LAI seasonality, which are both comparable with our experiment using MODIS climatology in static ECLand simulations. Here, we could not confirm that these findings are related to improved performance in heat fluxes since model performance of LE and EF is almost unaffected by the LAI forcing used, which has already been experienced by others (Stevens et al., 2020; Nogueira et al., 2021). The reason might be that parameters are adapted to the prior vegetation information (Ruiz-Vásquez et al., 2023), and thus the model needs a re-calibration.

For Noah-MP, activating vegetation dynamics had predominantly little impact on LE and EF predictions as well. A slight improvement in model performance was found for some sites with short-vegetation types or semi-arid climates. Ma et al. (2017) found that using LAI climatology resulted in better model performances for LE than simulations with activated vegetation dynamics for Noah-MP using monthly FLUXNET multi-tree ensemble data over the US. However, here we did not find enhanced biases in LE predictions with dynamic Noah-MP compared to the static simulations as they did, which could be due to the differing timescales for model evaluation. Both overestimation and underestimation of LE predicted by dynamic Noah-MP are reported in the literature (Brunsell et al., 2020; Ma et al., 2017; Cai et al., 2014). Brunsell et al. (2020) showed a positive bias of monthly evapotranspiration in the eastern Amazonian forest simulated with dynamic Noah-MP, while we found a negative bias of LE for the FLUXNET site GF-Guy. For short-vegetation types, using the Noah-MP crop module with activated vegetation dynamics might be sufficient in predicting surface fluxes (Liu et al., 2016). Achieved improvement for LE might not be as large as for sensible heat flux (Liu et al., 2016), which could be a reason for the poor performances in EF presented here.

Although vegetation and soil moisture state variables are directly coupled within land surface models, we found almost no impact of different vegetation modeling on soil moisture predictions for both models. Activating vegetation dynamics or changing LAI forcing did not improve soil moisture representation on average. The reason might be due to the implemented interaction of carbon and water processes. First, the potential photosynthetic activity that depends on leaf area and radiative conditions is calculated. Then, the limitation factor of extractable water is estimated according to available soil water and roots. Lastly, the photosynthetic activity is adapted to this restriction, and the transpiration rate is adapted to conductivity and atmospheric conditions. As a result, the only included sequence is that soil moisture impacts photosynthetic activity and biomass buildup. However, there is no feedback indicating that more biomass requires and loses more water, which would be taken from the soil because photosynthetic activity in the models relates to the carbon fluxes only, not the water fluxes. Additionally, modeled soil moisture suffers from substantial biases in both directions, which was also found by Liang et al. (2020) for Noah-MP and by Garrigues et al. (2021) for ECLand, although correlation between the observed and modeled soil moisture can be satisfactory (Beck et al., 2021; Xu et al., 2021; Pilotto et al., 2015; Liang et al., 2019). The reason might be underlying default values for soil characteristics, such as the field capacity and permanent wilting point, which possibly deviate from on-site soil conditions, and optimal values for soil parameters are still uncertain (Li et al., 2020). Alternatively, it could be an effect of differing scales since the observation from FLUXNET refers to point measurements. The multiscale parameter regionalization (MPR) might provide an improved way to estimate soil parameters by applying a pedotransfer function to local soil characteristics and has recently been applied to Noah-MP as a proof of concept (Schweppe et al., 2022).

Overall, model performance with respect to soil moisture and heat fluxes was barely affected by vegetation dynamics or the applied LAI forcing. However, sensitivity to LAI might be present, as van den Hurk et al. (2003) found effects of changed LAI values in TESSEL, a predecessor of ECLand, as did Ma et al. (2017) and Zhang et al. (2016) for Noah-MP. Xu et al. (2021) showed improved LE and soil moisture simulations with more realistic LAI, although the effect not only was site dependent but also differed with season and year. These authors also highlighted that transpiration is only partly determined by LAI, and other factors controlling the canopy conductance to water vapor might play a larger role. Therefore, other compensating mechanisms may explain the low elasticity between LAI and LE or soil moisture (see Sect. 3.2). Yang et al. (2011) demonstrated that the applied runoff scheme more strongly determined model performance with respect to soil moisture and evapotranspiration than the schemes for dynamic vegetation, stomatal resistance and soil moisture stress. Still, optimizing parameters can be effective in improving model predictions, which could be shown by several studies (Bohm et al., 2020; Li et al., 2021, 2020). Moreover, the sensitivity of soil moisture to vegetation parameters was shown to increase with dynamic vegetation representation (Arsenault et al., 2018). Yet, uncertainty remains regarding the optimal values for soil and vegetation parameters in particular (Li et al., 2020). Overall, the impact of vegetation dynamics and LAI on turbulent heat fluxes and soil moisture in this investigation was minor across sites and seasons for both models. Thus, modelers who are mainly interested in the performance of carbon processes should be careful when using performance metrics for hydrological variables as a proxy (e.g., LE) because the model formulation for the latter might have controlling processes other than LAI or NEE, which could dominate the results. Whether applying vegetation dynamics in model simulations is advantageous might depend on the target variables. While using MODIS climatology might be sufficient for heat flux predictions, activating vegetation dynamics could play a role in improving carbon flux predictions at seasonal or annual timescales (Jarlan et al., 2008).

4.3 Discrepancy between observed and simulated GPP–LAI relationship is caused by model structure

The substantial scatter in the observed relation between GPP and LAI is in close agreement with previous work, showing that GPP also depends on the short-term availability of resources (e.g., light, soil water) (Hu et al., 2022). Additionally, Zhang et al. (2021) found that in LSMs the relation between LAI and GPP was too close. We therefore checked the underlying relations in the models causing this. The GPP–LAI relationship in Noah-MP showed a clear exponential hysteresis (see Sect. 3.3). This hysteresis is related to the partitioning of GPP into the carbon pools in the plants and to LAI-reducing processes such as leaf turnover and leaf dieback. Noah-MP uses a non-linear function for allocation of GPP to the leaves that limits the maximum LAI the model can grow, resulting in the LAI saturation in summer that can be seen in Fig. 8m–p. On the other hand, leaf turnover due to leaf aging is implemented as a linear function of leaf mass, while leaf dieback due to environmental limitations follows exponential functions. Taken together, leaf dieback dominates in the later growing season, which results in the hysteresis. The reduction in LAI (i.e., leaf dieback) is implemented to be dependent on both water and temperature stress, but temperature stress is the main driver. In the specific implementation used here, water stress occurs only at a very low soil saturation exemplarily of 0.11 cm³ cm⁻³ for silt loam, which is even below the permanent wilting point of this soil texture type according to the lookup table value. These values are rarely reached, and thus water stress is negligible most of the time. In contrast, temperature stress is implemented as an exponential function, causing the late growing season non-linear decline in GPP observed throughout the non-tropical sites. Temperature stress is at maximum at 5 °C for forest ecosystems, resulting in no active biomass below this threshold. For this reason, LAI values are almost constant at the tropical forest site because temperature is never limiting there.

ECLand with static vegetation shows a similar pattern of seasonal dynamics to Noah-MP with vegetation dynamics but with a less pronounced exponential relationship. In contrast, dynamic ECLand simulates LAI that is strongly coupled to daily meteorological conditions, leading to higher daily fluctuations in LAI than expected, including strong drops in LAI in summer. Three processes govern these daily LAI dynamics: GPP, respiration and senescence. GPP is linearly related to LAI and varies with environmental and meteorological conditions, causing the variability in static runs. In dynamic runs, losses in biomass due to high or low daily GPP linearly affect LAI. In other words, unfavorable GPP can reduce LAI almost immediately. The second process affecting LAI is senescence. ECLand distinguishes between growing and senescence phases by comparing active biomass due to assimilation with the biomass from the previous time step. If active, then senescence is a linear function of active biomass and a folding factor. The folding factor reduces part of the senescent biomass, depending on photosynthesis (reduced in the case of high assimilation) and LAI. Overall, the folding factor changes only slightly with LAI. Additionally, a reduction in LAI and thus active biomass due to reduced GPP (as explained before) causes the model to trigger senescence because the active biomass of the previous time step was higher. The third process is respiration. About 11 % of physiologically possible assimilation is used for dark respiration without considering actual light conditions. This might cause high values of dark respiration compared with possible assimilation based on meteorological conditions and can thus reduce net primary production, even producing negative values. Notably, no aboveground biomass storage is built up, and there is no turnover. Most locations show a linear relationship comparable to ECLand but with a higher variability (Fig. 8 first and third rows). This might be due to the fact that leaf growth and leaf fall, in particular for trees, happen on longer timescales than the daily one as implemented in ECLand, which inhibits immediate effects of GPP on LAI. Overall, the current implementations of leaf dynamics in both models use very different approaches to represent LAI dynamics. In Noah-MP it is mainly temperature driven, and GPP depends minimally on LAI once the canopy is fully developed. In contrast, in ECLand, LAI and GPP are coupled very closely, and thus the LAI dynamics have almost the same sensitivities to water limitation and radiation as turbulent fluxes, which is unrealistic. Realistic LAI is less dynamic and less sensitive to environmental conditions, as also indicated by the observations. Hence for very different reasons, in both models the performance regarding LAI and turbulent and carbon fluxes is disconnected.

4.4 Implications and limitations

For the modeling of LAI and carbon fluxes, using dynamic vegetation modules in their current implementation in either model is not yet efficient because they increase model complexity, encompassing more dynamic processes and parameters without improving the predictive skill. As the dynamic vegetation components in ECLand are still under development, findings from this study can help better understand and represent the processes involved to improve its performance in modeling carbon and energy fluxes. For Noah-MP, we also showed that the dynamic vegetation module has potential for improvement, especially related to the relationship between GPP (and thus also NEE) and LAI. Underlying processes, such as carbon allocation, root dynamics, plant hydraulics, feedbacks on photosynthesis and their parameterization, can still be worked on (Ma et al., 2017; Li et al., 2021). Overall, we recommend using MODIS climatology forcing for static simulations, which yielded the best model performances for carbon and water fluxes. This might be valid for other remote sensing LAI products as well but would need to be tested beforehand. The value of a model evaluation like in this study depends on the reliability of the included datasets. Uncertainty in the forcing data might have a larger impact on the model runs than processes within the models (Zhang et al., 2016), but Haughton et al. (2016) demonstrated that observational errors, in general, are unlikely to cause poor model performance. Nonetheless, model evaluations are also restricted by uncertainty in the reference data (Li et al., 2022), especially when considering flux measurements (Li et al., 2019). We tried to address this by carefully inspecting the time series data from FLUXNET2015 before their usage. However, as in all measurements, there are still uncertainties, e.g., from instrumental errors or incomplete energy balance closure. Moreover, the MODIS dataset harbors uncertainty originating from cloud coverage, especially in the tropics. We tried to minimize this uncertainty by excluding all days from the dataset that were flagged with significant cloudiness. However, saturation also limits the representativeness of the LAI measurements. Even when using only data with the highest possible quality flag, we found suspiciously low LAI values in summer for temperate forests and grasslands and especially for tropical forests throughout the year (Fig. 8c).

Noisy and uncertain LAI data from MODIS for tropical forests have already been reported in the literature (Weiss et al., 2007; Garrigues et al., 2008; Xiao et al., 2016; Zhang et al., 2024). As a result, reference data remain a source of uncertainty, and a deviation in model outputs from these data is expected. In any case, reference data are essential for model verification, calibration and validation but should be treated carefully concerning their reliability and uncertainty.

5 Conclusions

Land surface models often include modules for dynamic vegetation processes. Yet, evaluations of the representativeness of key variables such as leaf area index or net ecosystem exchange are rarely done at high temporal resolution. The impact of different parameterizations of vegetation processes on water and carbon flux estimates by land surface models is still poorly understood. Therefore, we evaluated the change in model performance with respect to ecohydrological target variables when dynamic vegetation processes are included for two land surface models and further gained insight into critical process implementations that lead to the observed patterns. Surprisingly, including modules for dynamic vegetation in the model implementation did not improve the model predictions of ecohydrological variables for both ECLand and Noah-MP. We expected vegetation dynamics in these land surface models to better capture the higher variability in ecosystem exchange, especially that of highly dynamic short- or sparse-vegetation types, but this was predominantly not the case. Using alternative input for leaf area index to default climatology also had a negligible effect on the model performance, but this needs to be evaluated in more detail since our data sources were limited. Moreover, model performances in modeling carbon and hydrological fluxes appeared to be weakly coupled. Therefore, the question arose whether exchange fluxes themselves in these land surface models are sensitive to changes in leaf area index estimates and not only to changing parameter sets. Indeed, different leaf area index estimates lead to different predictions in exchange fluxes but without affecting the overall model performance for these variables. This might be caused by the mismatch in the seasonal patterns between observations and models for the relationship between gross primary productivity and leaf area index. While this relationship in dynamic Noah-MP showed a logarithmic hysteresis, mainly driven by temperature, both variables are closely linearly coupled in dynamic ECLand without allowing for the leaf area index to remain unchanged in suboptimal conditions for photosynthesis. This deeper analysis of the model performance for ecohydrological fluxes that pinpoints the reasons for model behavior was only possible with a reduced number of models. We used specific setups for the two land surface models evaluated here. Adapting or changing parameters and investigating the effect of other processes within the models were beyond the scope of this study. At this point, it remains unclear how representative our model selection is of the performance and process evaluation of other land surface models, since their processes are implemented differently. Nonetheless, we highlighted some crucial relationships in the implementation of vegetation processes that have the potential for further improvement. Additionally, these might be a good starting point for a similar intensive investigation with other land surface models or other alternative LAI climatology.

Appendix A: Dynamic ECLand processes

For more details, see previously published model descriptions (Boussetta et al., 2012, 2013, 2021). The photosynthesis model is based on Calvet et al. (1998). Therein, potential net assimilation A_n is estimated from physiological constrains as

$$\begin{array}{}\text{(A1)}& A_{\mathrm{n}}=A_{\max }\cdot \left(\mathrm{1}-e^{-\frac{g_{\mathrm{meso}}\cdot (c_i-c_{\mathrm{comp}})}{A_{\max }}}\right),\end{array}$$

where A_max is the leaf photosynthetic capacity, g_meso is the mesophyll conductance, c_i is the leaf-internal CO₂ concentration and c_comp is the CO₂ compensation point. Potential gross assimilation A_g is then calculated as

$$\begin{array}{}\text{(A2)}& A_{\mathrm{g}}=(A_{\mathrm{n}}+R_{\mathrm{d}})\cdot \mathit{ϵ},\end{array}$$

where R_d is the dark respiration from

$$\begin{array}{}\text{(A3)}& R_{\mathrm{d}}=A_{\mathrm{n}}\cdot f_{\mathrm{R}},\end{array}$$

where $f_{\mathrm{R}}=\frac{\mathrm{1}}{\mathrm{9}}$ is the dark respiration factor, and ϵ is a quantum use efficiency factor, estimated as

$$\begin{array}{}\text{(A4)}& \mathit{ϵ}=\mathrm{1}-e^{-\frac{{\mathit{ϵ}}_{\mathrm{0}}\cdot E_{\mathrm{PAR}}}{A_{\mathrm{n}}+R_{\mathrm{d}}}}.\end{array}$$

Here, ϵ₀ is the maximum quantum use efficiency, and E_PAR is the absorbed photosynthetic active radiation. Actual gross assimilation GPP results from

$$\begin{array}{}\text{(A5)}& \mathrm{GPP}=A_g\cdot \mathrm{LAI}\cdot {\mathit{\rho }}_{\mathrm{a}},\end{array}$$

where LAI is the leaf area index of the prior time step, and ρ_a is the air density corrected for humidity. A_n is used as the maximum leaf assimilation for the senescence model (Calvet and Soussana, 2001). To avoid immediate leaf dieback, a damping factor for senescence f_s is introduced as

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{f}_{\mathrm{s}}& =max\left({\displaystyle \frac{{\mathit{\tau }}_{\lim }\cdot t_{\mathrm{s}}}{\mathrm{100}\cdot N_{\mathrm{day}}}},max\left({\mathrm{10}}^{-\mathrm{8}},{\displaystyle \frac{t_{\mathrm{s}}}{N_{\mathrm{day}}}}\cdot \right.\right.\\ & \left.\left.min\left(\mathrm{1},{\displaystyle \frac{A_n}{A_{\max }}}\right)\cdot {\displaystyle \frac{max\left(\right(r_{\mathrm{meso}}\cdot \mathrm{1000}{)}^{\mathrm{0.321}}\cdot \mathrm{LAI},\mathrm{1})}{f_{\mathrm{LAI}}}}\right)\right)\end{array}\end{array}$$

where τ_lim is a limiting factor for immediate biomass loss, t_s is the damping time for senescence (which is basically the number of seconds per year), N_day is the number of seconds per day, A_max is the maximum photosynthesis rate with optimal conditions and f_LAI is a LAI correction parameter that reduces mortality at high LAI values which would occur due to shadowing. The amount of biomass loss B_loss is then

$$\begin{array}{}\text{(A7)}& B_{\mathrm{loss}}=min\left(B-{\mathrm{LAI}}_{\min }\cdot f_{\mathrm{LAI}-B},B\cdot \left(\mathrm{1}-e^{-\frac{\mathrm{1}}{f_s}}\right)\right),\end{array}$$

where B is the biomass of the prior time step, and f_LAI−B is a conversion factor between LAI and B. Biomass B is then updated by subtracting B_loss. The change in biomass due to assimilation B_gain results from

$$\begin{array}{}\text{(A8)}& B_{\mathrm{gain}}=max\left({\mathrm{LAI}}_{\min }\cdot f_{\mathrm{LAI}-B}-B,A_{\mathrm{n}}\cdot f_{\mathrm{Cbiom}}\right),\end{array}$$

where f_Cbiom≈0.68 is a factor converting the amount of CO₂ uptake from assimilation to carbon in dry biomass. Biomass B is updated again by adding B_gain. In the end, this updated biomass is transferred to an updated LAI value by

$$\begin{array}{}\text{(A9)}& \mathrm{LAI}={\displaystyle \frac{B}{f_{\mathrm{LAI}-B}}}.\end{array}$$

LAI determines the interception reservoir W by

$$\begin{array}{}\text{(A10)}& W=W_{\max }\cdot (c_{\mathrm{B}}+c_{\mathrm{H}}\cdot {\mathrm{LAI}}_{\mathrm{H}}+c_{\mathrm{L}}\cdot {\mathrm{LAI}}_{\mathrm{L}}),\end{array}$$

where W_max is the maximum thickness of the water layer on leaves or bare ground; c_B, c_H and c_L are the fractions for bare soil, high vegetation and low vegetation on a grid cell, respectively; and LAI_H and LAI_L are the LAI values for high and low vegetation, respectively (Boussetta et al., 2012). Additionally, canopy resistance r_c depends on LAI via

$$\begin{array}{}\text{(A11)}& r_{\mathrm{c}}=f_{\mathrm{1}}{f}_{\mathrm{2}}{f}_{\mathrm{3}}\cdot {\displaystyle \frac{r_{\mathrm{s},\min }}{\mathrm{LAI}}},\end{array}$$

where r_s,min is the minimum stomatal resistance, and f_n is the restriction factor for low input in shortwave radiation, soil moisture stress and saturated atmospheric conditions (Boussetta et al., 2012).

A1 Dynamic Noah-MP processes

For more details, see previously published model descriptions (Niu et al., 2011; Ma et al., 2017; Oleson et al., 2010). The model for leaf dynamics within Noah-MP is based on Dickinson et al. (1998). Leaf biomass C_leaf is balanced over time with

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\mathit{\delta }{C}_{\mathrm{leaf}}}{\mathit{\delta }t}}=f_{\mathrm{leaf}}\cdot A_{\mathrm{tot}}-(d_{\mathrm{stress}}+d_{\mathrm{turnover}}+R_{\mathrm{leaf}})\cdot C_{\mathrm{leaf}},\end{array}$$

where A_tot is the total carbon assimilation rate; f_leaf is the fraction of allocation to the leaves; d_stress is the dying rate caused by cold and drought stress; d_turnover is the turnover rate due to senescence, herbivory loss, or mechanical loss as a vegetation-type-dependent parameter; and R_leaf is the respiration rate of the leaf biomass. f_leaf is determined by LAI via

$$\begin{array}{}\text{(A13)}& f_{\mathrm{leaf}}=e^{\mathrm{0.01}\cdot \mathrm{LAI}(\mathrm{1}-e^{\mathit{\chi }\cdot \mathrm{LAI}})},\end{array}$$

where χ=0.75 is a parameter defining the partitioning of carbon allocation between the leaves and stem. A_tot is split up into photosynthesis rates from sunlit and shaded leaves, respectively:

$$\begin{array}{}\text{(A14)}& A_{\mathrm{tot}}=\mathrm{12}\cdot {\mathrm{10}}^{-\mathrm{6}}\cdot (A_{\mathrm{sunlit}}\cdot {\mathrm{LAI}}_{\mathrm{sunlit}}+A_{\mathrm{shaded}}\cdot {\mathrm{LAI}}_{\mathrm{shaded}}),\end{array}$$

where the first factor is for unit conversion. The partitioning of sunlit and shaded LAI results from a two-stream radiation transfer scheme (Niu et al., 2011). The assimilation rate for sunlit and shaded leaves is estimated with a bottleneck principle as

$$\begin{array}{}\text{(A15)}& A=I_{\mathrm{g}}min(A_L,A_C,A_S),\end{array}$$

where I_g is a growing season index according to leaf temperature, and A_L, A_C and A_S are the photosynthesis rates limited by light, RuBisCO and export, respectively (Bonan, 1996). A_L results from

$$\begin{array}{}\text{(A16)}& A_L={\displaystyle \frac{\mathrm{4.6}\cdot \mathit{ϵ}\cdot E_{\mathrm{PAR}}(c_i-c_{\mathrm{comp}})}{c_i+\mathrm{2}{c}_{\mathrm{comp}}}},\end{array}$$

with c_i being the leaf-internal CO₂ concentration, c_comp being the CO₂ compensation point, ϵ being the quantum use efficiency and E_PAR being the absorbed photosynthetic active radiation. Additionally, $A_S=\mathrm{0.5}\cdot V_{\max }$ and

$$\begin{array}{}\text{(A17)}& A_C={\displaystyle \frac{V_{\max }(c_i-c_{\mathrm{comp}})}{c_i+K_{\mathrm{c}}(\mathrm{1}+\frac{c_{\mathrm{o}}}{K_{\mathrm{o}}})}},\end{array}$$

where c_o is the atmospheric O₂ concentration; K_c and K_o are the Michaelis–Menton constants for CO₂ and O₂ (Collatz et al., 1992), respectively; and V_max is the maximum carboxylation rate, defined by

$$\begin{array}{}\text{(A18)}& V_{\max }=V_{\max ,\mathrm{25}}\cdot {\mathit{\alpha }}_{\max }^{\frac{T_{\mathrm{v}}-\mathrm{25}}{\mathrm{10}}}\cdot f_{\mathrm{N}}{f}_{T_{\mathrm{v}}}\mathit{\beta },\end{array}$$

where V_max,25 is the maximum carboxylation rate at 25 °C, α_max is a temperature conversion factor, T_v is the vegetation temperature, f_N is a factor for nitrogen limitation of the leaves, $f_{T_{\mathrm{v}}}$ is a factor for temperature limitation (Collatz et al., 1992) and β represents the limitation by available soil moisture. d_stress for the leaf mass balance is estimated from

$$\begin{array}{}\text{(A19)}& d_{\mathrm{stress}}=d_{\mathrm{cold}}\cdot e^{-\mathrm{0.3}\cdot max(\mathrm{0},T_{\mathrm{v}}-T_{\min })}{\displaystyle \frac{C_{\mathrm{leaf}}}{\mathrm{120}}}+d_{\mathrm{dry}}\cdot e^{-\mathrm{100}\mathit{\beta }},\end{array}$$

where T_min is a vegetation-type-dependent threshold temperature for leaf survival; β is the soil moisture limitation factor; and d_cold and d_dry are vegetation-type-dependent dying rates (prescribed parameter) for temperature and dryness stress, respectively. Leaf respiration R_leaf is calculated with

$$\begin{array}{}\text{(A20)}& R_{\mathrm{leaf}}=f_{\mathrm{res}}\left(f_{\mathrm{leaf}}-{\displaystyle \frac{\mathrm{LAI}}{\mathit{\chi }\cdot f_{\mathrm{leaf}}}}\right)\cdot A_{\mathrm{tot}}-R_{\mathrm{l}},\end{array}$$

where f_res is a factor defining the fraction of assimilation that is used for respiration, and R_l is the respiration for leaf maintenance from

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}{R}_{\mathrm{l}}& =min\left({\displaystyle \frac{C_{\mathrm{leaf}}-C_{\mathrm{leaf},\min }}{\mathrm{\Delta }t}},\mathrm{0.5}\cdot \mathrm{12}\cdot {\mathrm{10}}^{-\mathrm{6}}\cdot r_{\mathrm{l}}\left(T_{\mathrm{v}}\right)\right.\\ & \left.\cdot \mathrm{LAI}\cdot \mathit{\beta }\cdot {\displaystyle \frac{c_{\mathrm{N}}}{c_{\mathrm{N},\max }}}\right),\end{array}\end{array}$$

where C_leaf,min is the minimum leaf biomass, Δt is the time step duration, 0.5 is a reduction factor for respiration during the non-growing season, r_l(T_v) is the vegetation-type-dependent respiration rate for leaf maintenance at T_v and $\frac{c_{\mathrm{N}}}{c_{\mathrm{N},\max }}$ is the nitrogen saturation within the leaves. Afterwards, net primary production (NPP) is estimated as

$$\begin{array}{}\text{(A22)}& \mathrm{NPP}=\left(f_{\mathrm{leaf}}-{\displaystyle \frac{\mathrm{LAI}}{\mathit{\chi }\cdot f_{\mathrm{leaf}}}}\right)\cdot A_{\mathrm{tot}}-R_{\mathrm{leaf}}-R_{\mathrm{l}}.\end{array}$$

GPP is set to A_tot, and LAI is updated with

$$\begin{array}{}\text{(A23)}& \mathrm{LAI}=C_{\mathrm{leaf}}\cdot f_{\mathrm{LAI}-\mathrm{B}},\end{array}$$

where f_LAI−B is the leaf area per biomass. Assimilation rate A determines the stomatal resistance r_s by

$$\begin{array}{}\text{(A24)}& {\displaystyle \frac{\mathrm{1}}{r_{\mathrm{s}}}}=g_{\min }+{\displaystyle \frac{m\cdot p_{\mathrm{air}}\cdot A}{c_{\mathrm{air}}}}{\displaystyle \frac{e_{\mathrm{air}}}{e_{\mathrm{sat}}\left(T_{\mathrm{v}}\right)}},\end{array}$$

where g_min is the minimum stomatal conductance, m is an empirical parameter for the relationship between transpiration and CO₂ flux, p_air is the surface air pressure, c_air is the CO₂ concentration at the leaf surface, e_air is the vapor pressure at the leaf surface, and e_sat(T_v) is the saturation vapor pressure inside the leaves (Ball et al., 1987; Bonan, 1996). r_s is then used to estimate latent heat flux and thus evapotranspiration.

A2 Model setup options

Table A1 Options chosen for ECLand processes.

(Niu et al., 2007)(Yang and Niu, 2006)Sakaguchi and Zeng (2009)

Table A2 Options chosen for Noah-MP parameterization.