Assessing forest biomass and productivity with data-driven vegetation indices: insights from 900 000 simulated forest stands

Vegetation indices (VIs) are widely used in remote sensing, but deriving novel VIs for estimating specific forest attributes remains challenging. Here, “data-driven” VIs, yielding information based on correlations identified in large datasets of forest and reflectance data, could help. In this study, we systematically consider simulated reflectances of temperate forests (400–2400 nm range) and evaluate their correlations to above-ground biomass, leaf area index (LAI), yearly gross primary production (GPP), and yearly net primary production (NPP) production (NPP).

Considering 900 000 forest stands simulated via a classical forest model in combination with a radiative transfer model, we found that data-driven VIs could provide highly accurate estimates for the four analyzed forest attributes. Particularly VIs combining near infrared with shortwave infrared reflectances yielded good estimates. The wavelength combinations best suited for estimating above-ground biomass, LAI, and GPP showed considerable overlap. For less dense and structurally heterogeneous forests, reflectances from the visible band gained importance.

We introduced a new class of “non-parametric” vegetation indices and compared them with linear indices in scenarios of different environmental and physiological variability. Both the functional form of the VIs as well as the variability did not primarily affect the achievable accuracy of the model estimates, rather than the range of wavelengths from which good indices could be constructed. This suggests that data-driven vegetation indices can yield valuable results if the wavelength choice is optimized. These findings open new pathways for utilizing recent hyperspectral satellite missions such as EnMAP or CHIME.

1 Introduction

Temperate forests are crucial for carbon cycling and climate regulation in Central Europe, while also delivering a wide range of essential ecosystem services. However, temperate forests are also increasingly affected by climate change, which alters forest structure, productivity and species composition through rising temperatures, drought and altered disturbance regimes (Patacca et al., 2023). Therefore, monitoring the structural and functional properties of these forests is essential for understanding their current state and anticipating future changes (Coops, 2015; Ferretti et al., 2024). Key forest characteristics include (1) the above-ground biomass, needed for estimating a forest's carbon stock, (2) the leaf area index (LAI), relevant for assessing a forest's physiological processes such as tree growth or evapotranspiration, (3) the gross primary production (GPP), denoting how much carbon a forest takes up per unit of time, and (4) the net primary production (NPP), measuring the difference between GPP and autotrophic respiration in a forest per unit time.

In remote sensing, such forest characteristics are often assessed via vegetation indices computed based on forests' reflectance (Xiao et al., 2019). Many vegetation indices for analyzing forest cover and productivity consider reflectance from the red and infrared range, as leaves strongly absorb red light but reflect infrared light, making the difference between red and infrared reflectance a good indicator for leaf coverage (Zeng et al., 2022). With the start of the MODIS satellite, global data products for GPP and NPP were developed by combining these data with meteorological data and forest models (Running et al., 2004). These approaches were later extended by integrating calibration data from Eddy-Covariance measurements, which provide highly detailed information on GPP and net ecosystem exchange (NEE) at field sites (Joiner et al., 2018; Badgley et al., 2019; Xiao et al., 2019).

Although they were used as inputs in these complex analysis pipelines, vegetation indices were not originally designed to directly infer forest characteristics from reflectance data. Instead, they served as more general proxies, sometimes with well-understood limitations. For example, it remained challenging to estimate biomass from spectral data in dense forests, because the forest height, a significant driver of biomass, is difficult to infer from a forest's reflectance. Moreover, leaves in the lower canopy of dense forests only have a minor effect on a forest's reflectance, resulting in saturation effects when determining the LAI. Therefore, other methods such as LiDAR or RADAR as well as machine learning approaches deriving forest attributes from high resolution images have gained importance but are not available across the entire globe (Dubayah et al., 2022; Santoro et al., 2021; Lang et al., 2023; Mugabowindekwe et al., 2023). GPP and NPP, in turn, are influenced by dynamic physiological and environmental factors as well as canopy structure, making these variables even more challenging to derive from reflectance data than purely structural attributes. NPP is especially difficult to estimate, as it is also affected by the respiration of trees. This has led to critique of the MODIS NPP product (Turner et al., 2006; Park et al., 2021).

A pathway for revising vegetation indices is to consider information encoded in a broader spectrum of light and to develop more intricate mathematical linking functions to estimate forest characteristics directly from these data. For example, Brown (2000) included reflectance from the shortwave infrared (SWIR) range to improve estimates of the LAI. Several studies used simple mechanistic models to derive functional forms for vegetation indices to achieve accuracy in predicting GPP and other forest characteristics (Jin and Eklundh, 2014; Badgley et al., 2019; Camps-Valls et al., 2021).

The novel availability of hyperspectral data, such as from the German EnMAP (Environmental Mapping and Analysis Program) satellite launched in 2022 or the planned CHIME (Copernicus Hyperspectral Imaging Mission for the Environment) mission (Guanter et al., 2015; Qian, 2021; Storch et al., 2023), and machine learning methods allow extending such approaches. The new data and methods make it feasible to develop vegetation indices that harness more subtle relationships between forests and their reflectance spectra – including those arising from a complex interplay of ecological processes and constraints and therefore difficult to derive analytically. To infer these relationships, systematic data-driven approaches are needed.

Such an analysis needs to account for the effects of biome, site and environmental conditions, forest structure, plant traits, phenology, viewing geometry and sensor noise (Dong et al., 2024). These factors can influence which combination of wavelengths and linking function are best suited for estimating a certain forest characteristic. For example, variations in leaf chlorophyll content propagate to variations in blue and red reflectances, and may thus blur relationships between blue/red reflectance and LAI. Therefore, a particularly robust index may be needed if, e.g., the forests' phenological state is inhomogeneous. Thus, a systematic assessment of spectral robustness is needed to identify vegetation indices that are accurate and reliable under heterogeneous field and observation conditions.

The aim of this study is hence twofold. First, we provide a methodological framework for systematically evaluating data-driven two-band vegetation indices, including a novel class, which we call “parameter-free” indices. That way, we assess the potential of data-driven vegetation indices for estimating different forest characteristics. Second, we aim to determine the most informative and robust wavelength combinations for predicting biomass, LAI, yearly GPP and yearly NPP in the presence of environmental and physiological variability, and we analyze how these wavelength combinations change dependent on the density and structural complexity of the considered forests.

This analysis requires a large-scale dataset of forest and reflectance data, covering the whole set of possible forests within the range of expected environmental conditions and plant traits. Because such datasets are often not available for NPP and other forest characteristics of interest, we bridged the data gap using forest models, which already have a long and successful tradition in ecological research (Fischer et al., 2016; Rödig et al., 2018; Shugart et al., 2018; Maréchaux et al., 2021; Bugmann and Seidl, 2022). Simulation models have also been used in several remote sensing studies (see e.g. Fang, 2003; Lemaire et al., 2008; Wang et al., 2022; Dong et al., 2024). We simulated and analyze a total of 900 000 Central European forest stands along with their reflectance profiles subject to different sources of variation and noise. The resulting dataset contains reflectance data for more than 250 different wavelength channels along with corresponding forest properties and is published along with this study for further analysis.

2 Methods

Figure 1 Flow chart depicting the development of data-driven vegetation indices for biomass, leaf area index (LAI), gross primary production (GPP), and net primary production (NPP). The entropy of tree diameters (DBH entropy) is computed as measure of structural diversity to analyze the relationship between forest structure and optimal vegetation indices.

We analyzed the relationships between forests' solar reflectance spectra and their biomass, LAI, yearly GPP, and yearly NPP in scenarios with different heterogeneity of environmental and physiological parameters. For each scenario, we generated 100 000 forest states via a Markow Chain Monte Carlo (MCMC) approach applied to a forest model, and simulated the corresponding reflectance spectra using the radiative transfer model mSCOPE (Yang et al., 2017). Then we fitted linear and parameter-free models to investigate the links between the forest characteristics and the reflectance values at different wavelengths and compared their predictive capabilities to those of existing vegetation indices. To investigate how different sources of variation affect the investigated relationships, we applied variation terms to the parameters of the radiative transfer model in each scenario and analyzed their effects on the results. The process is depicted as a flow chart in Fig. 1. Below we provide details on each of these steps.

2.1 Generating forest stands

We used the process-based forest model Formind (Bohn and Huth, 2017; Fischer et al., 2024) along with a “ForestFactory” MCMC approach to generate a sample of forest stands representing the variety of forest compositions potentially found in the field, building on the work by Henniger et al. (2023b). Formind considers trees on the individual level and computes their gross and net productivity based on environmental conditions (e.g. temperature and precipitation), interactions between trees (e.g. light and water competition), as well as species-specific traits. As the model represents productivity and carbon fluxes explicitly and features both individual trees and their vertical leaf distribution, Formind is particularly suited to assess forest productivity and structure (Bohn and Huth, 2017; Fischer et al., 2024).

Before applying the model, we updated its parameterization for temperate forests in central Europe to better capture the carbon dynamics of heterogeneous forests. We considered forest inventory data collected at the “Hohes Holz” research site in central northern Germany (mean ann. temperature 10.8 °C; mean ann. precipitation 660 L m⁻²) and changed the temperate model parameterization (Bohn et al., 2014) in the following aspects: (1) we updated the allometric relationships according to data from Jucker et al. (2022); (2) we adjusted parameters for light climate and respiratory losses to better model the carbon dynamics in structurally heterogeneous forests; (3) we chose parameters on soil properties that matched the soil at the research site (values taken from Maidment, 1993); (4) we added the mechanistic defoliation and mortality mechanism introduced by Fischer et al. (2024); and (5) we used a new approach to consider weather data on a daily time scale whereas forest dynamics such as growth and mortality remain modelled on the yearly time scale. Details on the parameterization can be found in Sect. S5 of the Supplement. We validated the parameterization by comparing the simulated GPP, NPP, and respiration with independent estimates based on Eddy-Flux measurements (Pohl et al., 2023; see Sect. S5.3). A more in-depth validation study of an older parameterization of Formind was conducted by Holtmann et al. (2021) at the Hohes Holz site and by Rödig et al. (2017) at the Wetzstein site (Thuringia, Germany).

Typically, forest models such as Formind are initialized with existing forest inventory data or simulate forest successions starting from bare ground. However, since these simulations may not cover the entirety of forest states found in managed forests, we used a different approach here. Trees cannot survive if their GPP is insufficient to cover their respiratory needs. Hence, assuming that the GPP of any live tree found in a real forest must exceed its respiration (Bohn and Huth, 2017), we randomly assembled a sample of states satisfying this property, starting from a uniform prior of forest states.

We considered forest patches of 20 m×20 m size, for each randomly selecting a one-year sequence based on the climate at the Hohes Holz site (years 2000–2017; Muñoz Sabater, 2019). Starting from bare ground, we now randomly added or removed trees of random sizes and species. Then we computed the yearly GPP and autotrophic respiration of the forest stand given the climate sequence and evaluated whether the stand was “feasible”, i.e., (1) each tree's GPP exceeded its respiration, and (2) the tree crowns did not exceed the available space. If the stand was feasible, we repeated the tree addition/removal procedure; otherwise we returned to the last feasible state and repeated the attempt to add and remove trees until a feasible state was found.

To increase the range of qualitatively different forest stands, we used a hierarchical sampling approach when adding new trees, randomly constraining the potential tree sizes and species in some stands (Henniger et al., 2023b). Furthermore, we facilitated computational efficiency by dynamically adjusting the “step size”, i.e., the number of individual tree additions or removals between two forest feasibility checks: we started at 5 additions/removals and increased or decreased this number by factor 2 if the last five feasibility checks were positive or negative, respectively. After 200 tree additions or removals, we terminated the process. Further details regarding this approach are provided in Sect. S2.

We generated 100 000 forest stands for environmental conditions similar to those found at the Hohes Holz site. The forest stands differ, among others, in biomass, tree number, and species composition. For each stand, we computed the biomass, LAI, yearly GPP, and yearly NPP (mean values: biomass: 321 t ODM ha⁻¹, i.e., tons organic dry matter per hactar, LAI: 3, GPP: 28.8 t ODM ha⁻¹ yr⁻¹, NPP: 3.4 t ODM ha⁻¹ yr⁻¹). Furthermore, we determined the basal-area-weighted DBH entropy for the tree size distribution (Fischer et al., 2024), which is a measure for the structural heterogeneity of forests (mean: −4.54). Here, a lower value indicates lower heterogeneity and/or stronger dominance of a single individual. We provide the formula of basal-area-weighted DBH entropy in Sect. S4. The distribution of the considered forest characteristics is visualized in Sect. S3. Lastly, we computed the forests' solar reflectance spectra using the radiative transfer model mSCOPE, as described below.

2.2 Simulating solar reflectance spectra

We simulated the solar reflectance spectra of the generated forest patches for wavelengths between 400 and 2400 nm. To that end, we used an extended version of the radiative transfer model mSCOPE, designed for application in forests with vertically heterogeneous leaf traits (Yang et al., 2017; Henniger et al., 2023a). Our extended model version adds on to the original in also considering vertically heterogeneous leaf densities. Furthermore, we adjusted the probability of observing sunflecks on the ground, as the structure and covering of the ground are difficult to model exactly.

We parameterized the radiative transfer model similar to Henniger et al. (2023a) but adjusted individual parameters based on data from the TRY plant trait database (Kattge et al., 2011). A list of the parameters can be found in Sect. S1. We set the view zenith angle to 0° and the sun zenith angle to 31.5°, corresponding to maximal zenith angle observed at the Hohes Holz site. Furthermore, we set the relative azimuth angle between sun and view direction to a value of 140° (Henniger et al., 2023a), but note that this value is only relevant if the view zenith is changed from its default 0°. We set the probability to observe direct sunlight (sunflecks) at the ground to 0, which reduced the simulated reflectance values in the visible range. This improved the agreement between model results and field observations at the study site.

We simulated forests and corresponding hyperspectral data for 9 scenarios with different dimensions of variation each. First, we considered the scenario without any parameter variations. Second, we varied the LAI input to the radiative transfer model in each simulated forest patch by a random factor $\mathit{\epsilon }\in \left[\frac{\mathrm{2}}{\mathrm{3}},\frac{\mathrm{3}}{\mathrm{2}}\right]$ from a log-uniform distribution: $\ln \left(\mathit{\epsilon }\right)\sim \mathcal{U}\left(\ln \left(\mathrm{2}/\mathrm{3}\right),\ln \left(\mathrm{3}/\mathrm{2}\right)\right)$. Here, 𝒰(a, b) is the uniform distribution in [a,b]. Using a log-uniform multiplicative perturbation ensures that the parameters remain in the positive domain and in their respective order of magnitude.

Third, we perturbed all species' leaf trait parameters (Chlorophyll a+b content, leaf mass per unit area, equivalent water thickness, senescence material, carotenoid content) independently by random factors $\mathit{\epsilon }\in \left[\frac{\mathrm{2}}{\mathrm{3}},\frac{\mathrm{3}}{\mathrm{2}}\right]$ (log-uniform distribution). Fourth, we altered the leaf structure parameter in each patch by a random factor $\mathit{\epsilon }\in \left[\frac{\mathrm{2}}{\mathrm{3}},\frac{\mathrm{3}}{\mathrm{2}}\right]$ (log-uniform distribution), but constrained the resulting values to the admissible range [1,3] where necessary. Fifth, we changed the leaf inclination parameters by random values $\mathit{ϵ}\sim \mathcal{U}\left(-\mathrm{0.2},\mathrm{0.2}\right)$ and constrained them to the admissible range $\left[-\mathrm{1},\mathrm{1}\right]$ if required. Furthermore, we normalized these parameters by their joint Euclidean norm if the norm exceeded 1.

Sixth, we drew the soil wetness parameter in each simulated patch randomly from 𝒰(0, 1). Seventh, we changed the sun and view zenith angles as well as the relative azimuth by independent random values from $\mathcal{U}\left(-\mathrm{15},\mathrm{15}\mathit{°}\right)$. Eighth, we imposed independent random perturbations to all simulated reflectance values. The perturbations followed normal distributions with mean 0 and standard deviation ${\mathit{\sigma }}_i=\mathrm{0.05}{R}_i+\mathrm{0.005}$, where i indicates the wavelength and R_i the corresponding reflectance value. Finally, ninth, we combined all the perturbations listed above.

2.3 Analyzing vegetation indices

For each considered variation scenario, we generated 100 000 forest states with corresponding reflectance spectra and filtered out all forest states with a biomass below 50 t ODM ha⁻¹, as the hyperspectral model was primarily designed for areas with full forest cover. This concerned less than 2 % of the forest patches. We split the datasets into 75 % training and 25 % validation data. Then we selected all possible pairs of two wavelengths and fitted a model estimating either biomass, LAI, GPP, or NPP based on the corresponding reflectance values of these wavelengths in the training data.

Here, we considered two types of models: (1) an affine-linear model of the form $f=a_{\mathrm{0}}+a_{\mathrm{1}}{R}_{w_{\mathrm{1}}}+a_{\mathrm{2}}{R}_{w_{\mathrm{2}}}$, where f is the forest attribute, w₁ and w₂ are the wavelengths, and $R_{w_{\mathrm{1}}}$ and $R_{w_{\mathrm{2}}}$ are the corresponding reflectance values, and (2) a parameter-free model based on regular binning of the data. The shape of parameter-free models is not given by a simple formula but directly derived from the data, making these models particularly useful for data-driven vegetation indices. To obtain the parameter-free model, we binned the data into 25 regular intervals along each wavelength axis, resulting in 625 potential bins for the tuples $\left(R_{w_{\mathrm{1}}},R_{w_{\mathrm{2}}}\right)$. Then we determined the mean value of the considered forest property in each bin and assigned this value to the middle point of this bin. That is, for a bin $\left[a_{w_{\mathrm{1}}},b_{w_{\mathrm{1}}}\right]\times \left[a_{w_{\mathrm{2}}},b_{w_{\mathrm{2}}}\right]$, we set

$$\begin{array}{}\text{(1)}& f\left({\displaystyle \frac{a_{w_{\mathrm{1}}}+b_{w_{\mathrm{1}}}}{\mathrm{2}}},{\displaystyle \frac{a_{w_{\mathrm{2}}}+b_{w_{\mathrm{2}}}}{\mathrm{2}}}\right):={\displaystyle \frac{\mathrm{1}}{\left|\mathcal{B}\right|}}\sum _{i\in \mathcal{B}}{f}_i,\end{array}$$

where ℬ is the set of all data points in the bin and |B| its cardinality. We then applied a bi-linear interpolation between these points to obtain model predictions for all other reflectance values. For reflectances outside the convex hull of the training data, we applied a nearest neighbour extrapolation based on the closest bin centre.

For both model types, we computed the Pearson R² of each fitted model based on the validation data. Furthermore, we created heatmaps for the R² values as function of the wavelengths. Afterwards, we determined for each forest characteristic the maximal obtained R² and the percentage P_0.9 of models achieving an R² of at least 90 % of the maximum. We computed P_0.9 by determining the surface area in the heatmap where models achieved a sufficiently large R² and dividing it by the total surface area of all considered bands (Table 1). Specifically, we defined a function f_0.9(w₁,w₂) with $f_{\mathrm{0.9}}\left(w_{\mathrm{1}},w_{\mathrm{2}}\right)=\mathrm{1}$ for wavelength pairs with sufficient R² and $f_{\mathrm{0.9}}\left(w_{\mathrm{1}},w_{\mathrm{2}}\right)=\mathrm{0}$ otherwise, Then, $P_{\mathrm{0.9}}=\int \int f_{\mathrm{0.9}}\left(w_{\mathrm{1}},w_{\mathrm{2}}\right)\mathrm{d}{w}_{\mathrm{1}}\mathrm{d}{w}_{\mathrm{2}}$, which we approximated via a midpoint Riemann sum (numerical integration) of f_0.9 with those wavelength pairs as sampling points for which we had computed reflectance values.

To measure how important each individual band (Table 1) is for obtaining a large R², we determined the fraction of all models that use a wavelength from the respective band and achieve a high R² (90 % of max). We plotted this quantity (colour coded by band) along with the respective maximal R² for each considered forest characteristic and variation scenario.

Table 1 Considered wavelength ranges. The water absorption bands are excluded, as no reliable reflectance data can be obtained for these wavelengths in practice. The last column provides the wavelengths used for representing the band when analyzing classical vegetation indices (Table 2). The representing wavelengths were chosen from within the bands that existing satellite missions (MODIS, Landsat) record to represent the considered ranges.

To understand the impact of forest structure on these results, we repeated the analysis in the scenario without variations after grouping the data by forest structure. Specifically we split the forest set into sparse and dense forests (biomass below and above 200 t ODM ha⁻¹, respectively) and structurally heterogeneous forests vs. forests dominated by few individuals (DBH entropy S_DBH less or greater than −2.5, respectively; cf. Fischer et al., 2024). That way, we obtained the maximal R² values and the most significant wavelength bands for forests with qualitatively different structural properties. To streamline our study, we focused our analysis on the NPP.

We visualized the functional form of the best models with heatmaps (scenario without variations). We set the obtained R² values into context with existing vegetation indices by determining the correlations of classical vegetation indices with biomass, LAI, GPP, and NPP based on our simulated data without noise. Here, we considered all vegetation indices listed by Zeng et al. (2022) that use reflectances from clearly defined wavelengths or bands and do not depend on a site-specific parameterization. By squaring the correlation coefficients, we obtained the R² values that could potentially be achieved with linear models predicting the forest properties based on the vegetation indices.

3 Results

The R² values of the analyzed classical vegetation indices and the best linear and parameter-free models are displayed in Table 2. For the scenario without variations, the new data-driven models achieve much higher R² values than the best classical vegetation indices: 0.75 (data-driven) vs. 0.34 (NDVIre) for the biomass, 0.97 (data-driven) vs. 0.6 (EVI) for the LAI, 0.75 (data-driven) vs. 0.55 (EVI) for the GPP, and 0.6 (data-driven) vs. 0.34 (CAI) for the NPP. The maximal R² achieved with linear models was almost equal to the value achieved with parameter-free models and in one case (GPP) even slightly higher.

Table 2 Squared correlation coefficients between different forest characteristics and vegetation indices listed in Zeng et al. (2022). The bold font highlights the respective highest values in each column, excluding the values obtained by models developed in this paper, for which the best R² values are given in the last two rows. The values were computed based on forests generated without considering environmental or physiological variations. The formulas and index names were taken from Zeng et al. (2022) with minor corrections. The specific wavelengths used for the named bands are provided in Table 1. ^* Note on NDWI: version by Gao (1996), not McFeeters (1996). n/a: not applicable.

While the maximal achieved R² values were similar for linear and parameter-free models, the range P_0.9 of wavelengths for which high R² values (90 % of max.) were achieved differed and was much smaller for linear models, especially for biomass (drop by 66 %) and NPP (drop by 94 %; see Fig. 2). In the scenario without environmental or physiological variation, the biomass and LAI could be estimated best based on NIR reflectance combined with a reflectance from narrow bands from the SWIR range: one band close to the water absorption band between NIR and SWIR1, one at the centre of SWIR1, one close to the water absorption band between SWIR1 and SWIR2, and one encompassing the second half of SWIR2. These wavelength ranges were also well suited for estimating GPP, for which, however, an additional range combining green and red with NIR wavelengths yielded high R² values as well.

NPP was best estimated based on wavelengths in the SWIR range combined with any other range. Specifically, combinations of the centre SWIR1 range with any visible or NIR wavelength, combinations of two significantly different SWIR1 wavelengths, and combinations of SWIR2 with visible wavelengths or SWIR1 wavelengths close to the first water absorption band permitted high R² values.

Figure 2 Wavelength combinations used in the best (a) linear and (b) parameter-free models for different properties of forests without environmental or physiological variation. Each colour corresponds to a different forest characteristic. The shaded areas show the combinations of wavelengths for which models with at least 90 % the R² value of the best 2-wavelength model could be constructed. The graphs on the x and y axis show the mean reflectance profiles of all considered forest stands.

Analyzing the best models for the considered forest characteristics yielded several similarities. The best models have in common that the estimated forest properties increase with the reflectance of one wavelength, while they decrease with the reflectance of the other (Fig. 3). The reflectance values we used for fitting the models were correlated and hence accumulated along increasing lines $l:R_{w_{\mathrm{2}}}=a_{\mathrm{0}}+a_{\mathrm{2}}{R}_{w_{\mathrm{1}}}$ in the two-wavelength space (see dots in Fig. 3; $R_{w_{\mathrm{1}}}$ and $R_{w_{\mathrm{2}}}$ are the reflectances of the considered wavelengths w₁ and w₂). The gradient of the forest properties was perpendicular to this line, indicating that (1) many different reflectance pairs could lead to the same forest property estimate if they lie on a line parallel to the correlation line l and (2) the models were sensitive to changes in individual reflectance values (other reflectance held constant).

Figure 3 The best parameter-free models for (a) biomass, (b) LAI, (c) GPP, and (d) NPP, respectively, for forests without environmental or physiological variations. The axes correspond to reflectance values (R_w is the reflectance of wavelength w in nanometers). The background colour depicts the model's prediction, i.e., the respective forest characteristic corresponding to a reflectance pair. The points correspond to simulated forest patches in the validation dataset: their positions correspond to their reflectances, their colours to their biomass, LAI, GPP, and NPP.

The effect of environmental and physiological variability on the potential R² was moderate in most of the considered scenarios (dots in Fig. 4). In five of the seven considered noise scenarios with only one dimension of variation, the maximal obtained R² was almost insensitive to noise and did not decrease by more than 0.07. For biomass, LAI, and GPP, this applied for variability in leaf traits, leaf structure, leaf angle, soil wetness, and sun and view angles. White noise had a slightly higher effect (R² reduction by up to 0.15), and added variability in the LAI was the individual source of variation with the highest impact (R² reduction up to 0.2 for biomass and GPP and 0.3 for LAI). For the NPP, LAI variability was less significant for the maximal R². Instead, the leaf parameters played the most important individual role (R² decrease by 0.16). For all forest characteristics, the combined noise scenario yielded the lowest R², which was about half as big as without variations, respectively.

The ranges of wavelengths that yield a near-optimal R² (bars in Fig. 4) were more strongly affected by variations than the optimal potential R² values. Here, variability in the leaf traits and the leaf structure had the strongest effect, in particular for biomass and NPP, where the fraction P_0.9 of evaluated models with an R² exceeding 90 % of the optimum dropped by more than 82 %. For the NPP, variability in the leaf inclination parameters had a similarly strong impact. The value P_0.9 decreased least in the presence of LAI and sun/view angle variations. In the joint noise scenario, P_0.9 was strongly reduced for all forest characteristics.

In most considered scenarios, added variations did not affect the combinations of wavelength bands in which models with a high R² could be found (colours in Fig. 4). However, in the presence of leaf trait variability, no models using visible and red edge bands achieved near-optimal R² values. For the NPP, variability in the leaf structure parameters had a similar effect. Furthermore, models with near-infrared light lost their predictive capabilities for NPP if more factors than LAI and sun/view angles were variable.

Figure 4 Impact of different sources of variation on the predictability of forest characteristics. Each bar and point corresponds to a source of variation. The circles depict the maximal possible R² values. The bars show the fraction of considered wavelength pairs for which a model with an R² of at least 90 % of the respective maximum can be constructed. Here, the colours depict the bands from which the corresponding wavelengths are taken. It is visible that while the maximal R² is only moderately sensitive to variability in environmental factors and leaf properties, the range of wavelengths for which a high R² can be attained decreases significantly in the presence of noise.

Filtering the considered forest stands and developing specialized models for forests with high/low biomass and/or structural diversity generally improved the R² values for NPP (black circles in Fig. 5). NPP estimates achieved particularly high R² values (>0.7) in forests with small biomass (<200 t ODM ha⁻¹), whereas the R² values were lower in denser forests (R²<0.63). Binning by structural diversity yielded more intricate results: if the forests were also filtered by biomass, focusing on forests with lower structural heterogeneity ($S_{\mathrm{DBH}}<-\mathrm{2.5}$) led to larger R² values. Specifically, NPP could be estimated with an R² of 0.81 in forests with low biomass and low structural heterogeneity, whereas the R² dropped to 0.49 in forests structurally diverse forests with high biomass. If no biomass filter was applied, however, NPP could be better estimated in structurally diverse forests.

Structural diversity had a larger effect than the biomass on P_0.9, the range of wavelengths based on which NPP estimates with a high R² (≥90 % of max.) could be achieved (coloured circles in Fig. 5). More models achieved a near-maximal R² in structurally homogeneous forests. The biomass range only made a strong difference for structurally homogeneous forests, where P_0.9 was significantly larger in forests with a low biomass.

Filtering by biomass and structural diversity had a significant effect on the bands that were best suited to estimate NPP (colours in Fig. 5). Most prominently, the percentage of high-achieving models using NIR light was much larger in forests known to have a small biomass (namely 51 %–70 %) or large structural heterogeneity (20 %–99 %). In contrast, NIR light was much less used in the best models for unfiltered forests (14 % of the high-achieving models) or forests with neither large structural heterogeneity nor small biomass (0 %).

Figure 5 Predictability of the NPP for subsets of the forest stands. In the second and third column, the forest stands are filtered by DBH entropy, where low values indicate dominance by a single tree or multiple similarly-sized trees, whereas high values indicate that the basal area is evenly distributed over trees of different size classes. In the first column, no filtering by DBH entropy is applied. Similarly, the second and third row correspond to forests with high and low above-ground biomass, respectively, whereas no filtering was applied in the first row. The diameters of the hollow black circles correspond to the respective maximal R² values that could be achieved for the dataset. The shaded areas are proportional to the fraction of considered wavelength pairs for which a high R² (> 90 % of maximum) could be achieved. The colours depict the bands from which the corresponding wavelengths were taken. It is visible that for dense forests, most models consider a wavelength from the near infrared band. In contrast, for sparser forests, the SWIR1 band was used more frequently. In general, it was easier to estimate the NPP of dense forests and forests dominated by trees from few size classes.

4 Discussion

In this study, we conducted an extensive in-silicio study of 900 000 simulated temperate forest stands to (1) assess the general potential of data-driven vegetation indices to estimate forest characteristics and (2) analyze which wavelengths pairs yield the best estimates depending on the variability of traits and conditions within the considered datasets. To that end, we systematically evaluated a total of more than 7 million data-driven models for biomass, LAI, yearly GPP, and yearly NPP based on different wavelength pairs. We found that even if the forests were subject to significant variation of environmental and physiological properties, data-driven models yielded relatively large R² values, significantly exceeding the R² values obtained from classical vegetation indices. This suggests that data-driven vegetation indices, even if using only two wavelengths, are a promising tool for deriving forest properties from remote sensing data.

4.1 Main findings

Neither a large number of different wavelengths nor a complicated functional form were necessary to estimate forest attributes with relatively high accuracy. In fact, linear models superseded parameter-free models without pre-imposed functional form in some wavelength regions. However, focusing on a specific functional form limited the range of wavelengths suitable for analyzing forest attributes and increased the challenge of identifying the optimal wavelengths. This agrees with earlier findings by Gong et al. (2003), who evaluated different wavelength combinations and functional forms to estimate LAI based on reflectance data. The better suitability of linear models in some wavelength regions is due to the limited resolution of the binning method we applied in the parameter-free case, where we computed average values of forest characteristics for individual reflectance intervals. Increasing the resolution (i.e., using smaller intervals) or using more general parametric models will therefore lead to even better R² values than we presented, but may require larger datasets.

Despite the relatively high R² values we obtained, ranging from 0.6 for the NPP to 0.97 for the LAI, our results show that identifying the right wavelengths for estimating specific forest properties is key. In the absence of environmental and physiological variability, the classical vegetation indices rarely utilized wavelength pairs optimal for estimating any of the considered forest characteristics. Exceptions were the vegetation indices using red and NIR/red edge light (e.g. NDVI, EVI, NDVIre, NIRv), which were, however, only in the optimal range for estimating GPP and neither of the other forest characteristics. This is surprising, since LAI and GPP are typically assumed to be strongly correlated (Gitelson et al., 2014), and several vegetation indices were developed by considering the reflectance profile of leaves (Zeng et al., 2022).

In the presence of environmental or physiological variability, the range of suitable wavelengths decreased significantly, and some bands lost their suitability. For example, visible light became sub-optimal when leaf traits were uncertain, and NIR light became unsuited for estimating NPP in this case. This suggests that the choice of wavelengths for estimating forest characteristics should take into account the type and magnitude of variations expected in the field data. Identifying and quantifying these sources of variation will hence be an important preparational step to deriving data-driven vegetation indices for forest characteristics.

Aside the above-mentioned cases, the best band combinations for estimating forest characteristics remained remarkably stable throughout the considered scenarios. In particular, combinations of longer wavelengths (NIR, SWIR1, and SWIR2) were generally well suited. This is in line with previous studies (Gong et al., 2003; Lemaire et al., 2008; Psomas et al., 2011; Houborg and McCabe, 2018; Almeida et al., 2019), yet contrasts with the most common choices of wavelengths for vegetation indices (only 4 out of 23 vegetation indices listed in the review of Zeng et al. (2022) combined two wavelengths from these bands). However, these vegetation indices were typically not designed to quantitatively estimate the forest characteristics considered in this study. Instead, many are applied to assess the presence of vegetation, e.g. classify areas into forest landscapes and other landscape types, whereas in our study, we considered forests only, i.e., presumed that the landscape had been filtered to only include forests before. Nonetheless, a classification into forested and unforested land could be easily conducted in an independent preparational step, making forest-specific vegetation indices applicable in combined multi-step procedures in practice.

We found that combinations of NIR and SWIR wavelengths were particularly useful for estimating forest characteristics. Reflectance in the SWIR range is strongly related to leaf water content, lignin, proteins, nitrogen and cellulose (Curran, 1989; Fu et al., 2021; Zeng et al., 2022), and the combination of LAI and nitrogen content has been found well suited for estimating forest productivity (Reich, 2012; Zhang et al., 2023). The significance of NIR/SWIR combinations may be further understood by considering the sensitivity of forest reflectance to LAI, which is low in the NIR range but high in the SWIR range (Verrelst et al., 2015; sensitivity analysis of radiative transfer models). Furthermore, NIR and SWIR reflectances are much less sensitive to leaf constituents such as chlorophyll or carotenoid than reflectances from the visible spectrum (Mousivand et al., 2014; Verrelst et al., 2015; Prikaziuk and van der Tol, 2019). This may also explain why visible light became unsuitable for estimating GPP when the leaf properties were uncertain. We note, however, that a discussion of potential technical limitations in measuring NIR/SWIR light or required atmospheric corrections is beyond the scope of this study.

Since variations in leaf constituents may be induced by high species richness, reflectances from the visible spectrum may be unsuitable for assessing particularly diverse forests. The significance of visible light decreased similarly in forests with high biomass, potentially due to the low penetration depth of visible light (Hovi and Rautiainen, 2020) leading to a quick saturation of reflectance-LAI relationships (Mutanga et al., 2023). These findings suggest that visible light is sub-optimal to assess the considered forest characteristics in structurally complex forests.

Noteworthily, an increase in the heterogeneity of tree sizes had a contrasting effect: it increased the significance of visible light for estimating NPP. The tree size diversity, measured by the DBH entropy in this study, is low if forest patches are dominated by individual large trees, overshadowing understorey trees. If multiple differently sized trees contribute equally to a forest's basal area, it is more likely that their canopies are directly visible from above, making it easier and potentially optimal to use visible light to assess the forest state. This relationship may not hold in strongly stratified (e.g. tropical) forests, though.

4.2 Methodological contributions

A major aim of this study was to establish a simulation-based methodological framework for systematically evaluating two-band vegetation indices. Though building on previous work on generating ecologically feasible forests with simulation models (Bohn and Huth, 2017; Henniger et al., 2023b) and combining forest models with radiative transfer models (Henniger et al., 2023a), we introduced several methodological advancements. Our Markov Chain Monte Carlo approach for sampling forests allows a rigorous statistical interpretation of the distribution of the resulting forests. Furthermore, the algorithm's simplicity reduces the risk of model artifacts potentially favouring specific forest structures, makes it easier to incorporate climatic effects on forest states and to exchange the underlying forest model with other process-based forest models. The updates of the forest model Formind and its parameterization for temperate forests improved the model's efficiency and accuracy in heterogeneous forests.

The adjusted version of mSCOPE disentangles the vertical space and LAI dimensions, thereby permitting the simulation of forests with heterogeneous leaf densities in different layers of the canopy. The new soil sunfleck parameter reduces the impact of uncertain soil properties on reflectance and may lead to more accurate results in scenarios with large solar zenith angles.

The parameter-free vegetation indices proposed in this study circumvent the need to limit the predictive power of vegetation indices by pre-imposing specific functional forms. This makes this class of vegetation indices particularly useful for data-driven vegetation indices, as the optimal functional form can be directly derived from the data. At the same time, their methodological simplicity makes data-driven indices well-suited for considering large data sets and analyzing many vegetation indices systematically. However, if less data are available or more than two or three wavelengths are regarded jointly, other methods, such as machine learning approaches, should be considered.

4.3 Limitations

We used a model-based approach for a systematic analysis of the vegetation indices under controlled conditions. This is a typical method in theoretical studies (see e.g. Fang, 2003; Lemaire et al., 2008; Wang et al., 2022; Dong et al., 2024) and facilitates a general understanding of the studied relationships, as a variety of different scenarios can be considered easily. Nonetheless, the results depend on the assumptions underlying the applied models and are hence subject to associated limitations.

The applied forest model Formind considered climatic drivers such as droughts and their effects, but the relationships between these stressors with the optical leaf properties (Watt et al., 2021; Zhou et al., 2021) were not modelled explicitly. Similarly, the effects of shadows occurring on the surface of heterogeneous tree canopies (Hilker et al., 2010; Zeng et al., 2022), which could affect the estimability of forest characteristics by inducing a relationship between forests' height heterogeneity and their reflectance, are not accounted for in mSCOPE. We addressed these simplifications by modelling and investigating the effect of different sources of variation on the results, and our findings are in good agreement with earlier analyses of wavelength combinations for estimating LAI (Gong et al., 2003).

We simulated the analyzed forest stands considering environmental conditions at the Hohes Holz site, where the utilized forest model had been extensively validated in multiple respects (Rödig et al., 2017; Holtmann et al., 2021). The “ForestFactory” approach to generate forests yields stands with a broad variety of structures (Henniger et al., 2023b). This makes the simulated dataset representative for a large set of forests varying in site conditions and management. We accounted for weather-induced variability by additionally randomizing the utilized weather data. Nonetheless, the model results may not be representative for forests with significantly different site conditions or climate (e.g. montane forests). Additional studies are required to further assess the impacts of a wider range of environmental site characteristics on vegetation index performance (e.g. soil type, rooting depth, slope and aspect, mean temperature and precipitation).

4.4 Outlook and future directions

Our systematic in-silicio analysis of forests and their reflectance provides the basis for deriving data-driven vegetation indices from field data. Field data for structural forest properties such as biomass and LAI could be derived from LIDAR measurements and orthophotos. These data could then be combined with reflectance data from multi- and hyperspectral satellite missions. Data-driven vegetation indices for biomass, LAI, and yearly GPP could be derived based on data from the Landsat 8/Landsat 9 satellites (most promising band combination: band 5: NIR and band 7: SWIR 2), the Sentinel 2 satellite (band 8: NIR/band 12: SWIR 2), or the MODIS instrument (band 5: NIR/band 7: SWIR). Here, the MODIS data may be particularly well suited due to the sensor's relatively small bandwidth in the SWIR range.

Deriving a data-driven vegetation index for NPP may be more challenging, as corresponding field data are difficult to acquire. Furthermore, the wavelengths we found to be best suited for estimating NPP are often not considered by multispectral satellite missions, because these wavelengths are close to the atmospheric water absorption bands. Here, data from hyperspectral satellites such as EnMAP and CHIME could prove particularly useful, as they record a significantly larger number of wavelength channels.

Since we found the sensitivity towards wavelength choice to be minimized by parameter-free vegetation indices, we suggest taking a corresponding approach when developing data-driven vegetation indices. Once data-driven indices are established for some areas where field data are available, the results could be used to recalibrate and fine-tune our model-driven analysis pipeline, e.g., by providing insights into the variability of environmental conditions and plant traits. The refined model could then be used to generalize the results from the field studies or to develop indices specifically tailored towards specific forest types.

These indices could be derived by simulating forests with structural properties mimicking the those of the respective region of interest. Tailor-made vegetation indices could yield more precise results than presented in this study, where we considered a broad mixture ranging from even-aged monocultures to uneven-aged heterogeneous forests. Moreover, since data-driven vegetation indices require that the training data are representative for the area of interest, using carefully constrained data for fitting may even be a prerequisite for accurate estimates.

Besides generalizing results from spatially constrained field studies, the ability to jointly simulate forests and their reflectance can support the development of advanced vegetation indices that incorporate more wavelengths simultaneously. To that end, machine learning methods may be applied, which require sufficiently large training datasets. This increases the importance of simulation approaches. Our process-based modelling approach can furthermore be extended to consider variables beyond those analyzed in this study. For example, it might be possible to develop specialized vegetation indices for forest properties such as species diversity, structural diversity, net ecosystem exchange, forest health or disturbance.

5 Conclusions

We used a hybrid modelling approach to generate large datasets of forest stands and corresponding hyperspectral data in the 400–2400 nm range. Based on these datasets, we systematically evaluated the potential of vegetation indices to estimate above-ground biomass, LAI, yearly GPP, and yearly NPP. We found that estimates from data-driven indices have the potential to be significantly more accurate than predictions derived from “classical” vegetation indices.

We assessed which wavelength combinations were best suited for estimating the considered forest characteristics and observed that combinations of NIR and SWIR light yielded good results in general, with biomass, LAI, and GPP often being well inferable via the same wavelength combinations. The optimal choice of wavelengths depended on the structure of the considered forests, with visible light gaining in importance in less dense and structurally heterogeneous forests.

We proposed and evaluated a new class of vegetation indices, namely parameter-free vegetation indices. We found that the functional form of the vegetation indices did not significantly affect the maximal possible achievable accuracy, but instead constrained the range of wavelengths where this accuracy could be attained. We obtained a similar result with respect to environmental and physiological variability: evaluating different potential sources of uncertainty in physiological and environmental parameters, we observed that while uncertainty did not strongly reduce the achievable accuracy, it decreased the range of wavelengths where accurate vegetation indices could be constructed. These results, along with the simulation approach introduced in this study and the generated data, may facilitate the development of new data-driven vegetation indices, optimized for estimating individual forest characteristics of interest and tailored to the structure of the considered forests.