3.1. Assessing model performance using accuracy metrics

Table 4 shows the root mean squared errors and coefficients of determination for all three different ML models and the two different input variable configurations. In general, the models performed better when trained on the full input set, as seen from the higher $ {R}^2 $ and lower RMSE. However, the $ {R}^2 $ is high, at almost 0.8, for reduced input set models as well. The difference between $ {R}^2 $ and RMSE for the different input set configurations is around $ \simeq $ 0.05–0.06 and $ \simeq $ 0.30–0.32, respectively. This difference is small enough to consider even the reduced input set models to model NEE well, but also high enough to indicate that the additional variables do contain information about NEE for slightly more accurate predictions.

For the reduced input set, all three models perform well, including RF, which is only slightly worse than the more advanced models XGBoost and LGBM. As RF is a slightly less complex model compared to XGBoost and LGBM, containing less tunable hyperparameters, this result is to be expected. For the full input set, RF is again the worst model, which is to be expected due to the previously stated reason; however, here the difference is noticeably larger, highlighting how the more complex models find more information from the full input set. Overall, it can be deduced that all models capture the relationship between NEE and the input variables well, and that these data are well suited for modeling NEE efficiently, even when using the reduced input set containing only four variables.

In a similar study, Abbasian et al. (Reference Abbasian, Solgi, Hosseini and Kia2022) used ML on flux-tower data to predict daily NEE in various European ecosystems. For evergreen needle-leaf forests, their models achieved $ {R}^2 $ values between 0.80 and 0.88. Our scores fall into the same range; however, we used slightly different input variables and instead modeled 30-min fluxes compared to daily fluxes. Another study by Liu et al. (Reference Liu, Zuo, Wang, Yuan, Zhu, Zhang, Zhang, Sun, Guo, Guo, Song, Song and Xu2021) showed similar ranges of $ {R}^2 $ of 0.89–0.90 for modeling NEE using ML in evergreen needle-leaf forests for daily fluxes. An important notion is that the dataset in the present study is dominated by Hyytiälä data, as seen from Table 2, and according to our previous study (Ezhova et al. Reference Ezhova, Laanti, Lintunen, Kolari, Nieminen, Mammarella, Heljanko and Kulmala2025a), Hyytiälä was modeled better than Värriö. This suggests that Värriö and Järvselja may be more challenging for the ML models, possibly due to factors such as limited data or higher variability, which might explain the reduced overall performance.

3.2. SHAP summary: visualization of feature effects

We now consider SHAP values, allowing us to analyze how models see the relationship between input variables and the target variable, NEE. With baseline $ {\phi}_0=-3.269194\mu \mathrm{mol}\hskip0.1em {\mathrm{m}}^{-2}\hskip0.1em {\mathrm{s}}^{-1} $ , SHAP values measure how each feature shifts the prediction from this mean: negative values indicate stronger uptake (more negative NEE), while positive values indicate stronger respiration (more positive NEE).

Here, we consider the SHAP values, using the beeswarm summary plot seen in Figure 1. This summary plot presents a compressed overview of how the input variables impact the model’s output. In order to make the plot easier to understand, all ordering is based on the results of the XGBoost model, the best-performing model.

Figure 1.

Beeswarm summary plots of SHAP values showing the influence of input variables on NEE predictions across models and input configurations. PAR, DiffPAR, NDVI, and soil temperature consistently show the strongest effects. Despite variations in model type and input set, overall SHAP patterns remain similar, indicating robust relationships between input variables and NEE.

Figure 1. Long description

A grid of six beeswarm plots organized into three rows and two columns.

Axes and Legend:

* The vertical y-axis lists features: P A R, Diff P A R, N D V I, Soil Temp, L A I, Air Temp, V P D, Rel Hum, Soil Wat Cont, and Fric Vel.

* The horizontal x-axis represents the S H A P value for N E E in micromoles per second per square meter, ranging from negative 10 to positive 5.

* A color scale at the bottom indicates Feature Value, ranging from dark purple for Low to bright yellow for High.

Grid Layout:

* Left Column (Panels A, C, E): Full input set containing all ten features.

* Right Column (Panels B, D, F): Reduced input set containing only the top four features: P A R, Diff P A R, N D V I, and Soil Temp.

* Top Row (Panels A and B): X G Boost model.

* Middle Row (Panels C and D): L G B M model.

* Bottom Row (Panels E and F): R F model.

Data Trends:

* Across all panels, P A R and Diff P A R show the widest horizontal spread, indicating the strongest influence on N E E. High feature values (yellow) for these variables generally correlate with negative S H A P values.

* N D V I and Soil Temp show moderate spread, with high values typically associated with positive S H A P values.

* The distribution patterns for the top four features remain highly consistent regardless of whether the model uses the full or reduced input set or which of the three model architectures is employed.

For the full input set (Figures 1a, 1c, 1e), all ML models agree that the four variables that have the largest scale of effect on NEE are PAR, DiffPAR, NDVI, and soil temperature as indicated by the highest mean absolute SHAP values (global importance). These variables can be considered to be the ones with the most explanatory power, with PAR being the clear first in each case. This ranking is expected: these variables control light quantity and quality for photosynthesis, temperature-driven respiration, and canopy greenness. Because the study sites are temperature-limited and relatively humid, VPD and water availability contribute less at the annual scale (Ezhova et al., Reference Ezhova, Aarne, Arola, Lipponen, Lintunen, Kokkola, Ylivinkka, Yli-Juuti, Petäjä, Kerminen, Virtanen and Kulmala2025b).

PAR is theoretically the most important variable to explain photosynthesis, especially during the peak growing season (Palmroth and Hari, Reference Palmroth and Hari2001; Moffat et al., Reference Moffat, Beckstein, Churkina, Mund and Heimann2010). These results are consistent with previous results (Ezhova et al., Reference Ezhova, Laanti, Lintunen, Kolari, Nieminen, Mammarella, Heljanko and Kulmala2025a) from the same forest sites, where PAR, DiffPAR, and SoilTemp were the main drivers. The effect of DiffPAR on photosynthesis is also well known (Gu et al., Reference Gu, Baldocchi, Verma, Black, Vesala, Falge and Dowty2002), with canopies with higher LAI benefiting more from diffuse radiation (Durand et al., Reference Durand, Murchie, Lindfors, Urban, Aphalo and Robson2021).

While the light-related variables are considered to be the most important to explain NEE, the third variable in importance is NDVI, most likely important due to it being a good proxy for the fraction of absorbed PAR (Balzarolo et al., Reference Balzarolo, Vicca, Nguy-Robertson, Bonal, Elbers, Fu, Grünwald, Horemans, Papale, Peñuelas, Suyker and Veroustraete2016), and reflecting seasonality and the average age of the leaves. NDVI is also distinct enough in the different sites to act as a possible site-identifier proxy. Soil temperature, on the other hand, can model respiration, as respiration is highly sensitive to temperature, while also capturing seasonality. While the model clearly uses air temperature to some degree as well, it has been noted in some previous studies that soil temperature may be a better parameter than air temperature for modeling ecosystem respiration (Kolari et al., Reference Kolari, Kulmala, Pumpanen, Launiainen, Ilvesniemi, Hari and Nikinmaa2009; Lasslop et al., Reference Lasslop, Migliavacca, Bohrer, Reichstein, Bahn, Ibrom, Jacobs, Kolari, Papale, Vesala, Wohlfahrt and Cescatti2012). Several of the remaining input variables are highly correlated with the four most important ones—for example, soil temperature is closely related to air temperature—which may cause the model to either prioritize one or split the explanatory power between them (Moffat et al., Reference Moffat, Beckstein, Churkina, Mund and Heimann2010). However, they still contribute to the model, and since SHAP can split contributions among correlated inputs (Molnar, Reference Molnar2025), the importance of these variables may be divided.

When considering the contributions of the individual feature values, the SHAP values demonstrate that they capture these contributions consistently, regardless of the model or input configuration used. The plots show that higher radiation drives more negative NEE (higher carbon uptake), while lower radiation produces higher NEE. This pattern aligns with both diurnal and seasonal cycles, so that the highest-radiation points (midday in peak growing season) contribute most strongly to NEE. NDVI behaves differently from radiation variables because it does not change over the course of a single day but instead reflects seasonality. Instead, NDVI’s influence is subtler: data points associated with more negative NEE (greater uptake) tend to have slightly higher NDVI than those linked to positive NEE, even though NDVI remains essentially constant between morning, midday, and evening. The fourth parameter is soil temperature, for which SHAP clearly shows that the dependence is nonlinear, as we see points with high soil temperature on both sides of the SHAP values. Leaf area index (LAI) is the variable that differs most across the sites, where low LAI values correspond exclusively to Värriö, which also has the lowest overall NEE from the studied sites due to it being in a much northern region, making it most likely the variable the models use to identify the site, similarly to NDVI. VPD, relative air humidity, soil water content, and friction velocity all contribute to the predictions on a small scale, and as previously mentioned, they are highly dependent on other variables, leading to them having a smaller scale of impact.

All three different ML models exhibit similar overall patterns in SHAP values, but LGBM shows more distinctive behavior, while XGBoost and random forest are more alike. The shapes and the colors of the distributions of the SHAP values are similar; however, the scales vary. PAR has a larger scale in LGBM (Figure 1c), providing SHAP values up to −9, while random forest and XGBoost have only values up to around −7.5 (Figure 1a, e). Reducing the input variables makes the scales of the four variable larger, as is expected due to the need to compensate for the effect of the removed variables, making the PAR SHAP values as high as −12 in LGBM and XGBoost (Figure 1b, d). It would seem that when using the reduced variable set, the positive effects of any features omitted from that set are absorbed by PAR and NDVI, while the negative effects are captured by soil temperature (Figure 1b, d, f). However, the distributions remain similar in shape and color, suggesting that they are still modeled largely in the same way, whether additional variables are included, suggesting a robust relationship between the main input variables and NEE. Due to the large dimensionality of the SHAP values, only general trends can be observed from Figure 1, highlighting the need for a different kind of visualization to observe more detailed relationships, provided in the next section.

3.3. Dependence of NEE on different variables visualized with SHAP values

In this section, we consider SHAP values based on individual features and plot them using scatter plots and box plots to analyze the interactions between different variables. While previous studies exist that use SHAP values to analyze ML models of NEE (Gao et al., Reference Gao, Yao, Gao, Yuan, He, Wang, Li, Ma, Russell, Pressley and Zou2024; Chen et al., Reference Chen, Wu, Chen and Wan2025), they mostly focus on an overview of the variable effects using plots similar to those seen in Figure 1. While this type of figure effectively captures the general relationship between input variables and the target variable through SHAP values, it does not reveal how features interact with one another or the underlying nature of these dependencies.

We can analyze the impact of variables of interest on NEE in more detail. All SHAP values in this and the following section come from the XGBoost model trained on the reduced input variable set (PAR, DiffPAR, NDVI, and Soil Temperature). To reveal how variables interact, the scatter plots are colored according to another selected variable. This approach highlights how the SHAP values of one predictor change in relation to variations in another, making the dependencies between drivers more visible.

Figure 2 shows the PAR, DiffPAR, NDVI, and SoilTemp on the X axis with the corresponding SHAP values on the Y axis. The coloring of the plots was chosen based on which variable interactions we wanted to observe in more detail. Figure 2a, featuring PAR and PAR SHAP, shows the slightly nonlinear interaction between PAR and NEE. As PAR increases, so does its contribution toward more negative NEE. This is consistent with increasing radiation contributing toward more photosynthesis (and thus reduced NEE) as mentioned in the previous subsection. NDVI was chosen as the coloring to demonstrate how seasonality interacts with PAR. A low NDVI, corresponding to points in winter, early spring, or late autumn, restricts the effect of PAR, as lower leaf biomass means less absorption of radiation and less photosynthesis, but also less respiration. Conversely, high NDVI values correspond to the highest range of SHAP change across PAR levels, meaning that the highest increase in carbon uptake due to PAR is reached when the NDVI is high, but also the highest emissions of carbon under low PAR occur when the NDVI is high. The PAR SHAP values have the highest absolute effect range on NEE, reflecting the perceived importance of the variable.

Figure 2b visualizes DiffPAR and DiffPAR SHAP. It demonstrates that as DiffPAR increases, so does the carbon sink. Here we can also see how soil temperature interacts with the effect of DiffPAR on NEE: lower soil temperature is linked to more negative NEE when DiffPAR is low, around <200 $ \mu \mathrm{mol}\;{\mathrm{m}}^{-2}\;{\mathrm{s}}^{-1} $ , while higher soil temperature is linked to more negative NEE when DiffPAR is high, around >300 $ \mu \mathrm{mol}\;{\mathrm{m}}^{-2}\;{\mathrm{s}}^{-1} $ . This interaction reflects seasonal patterns, where DiffPAR has little effect on NEE under cold soil conditions, likely due to limited photosynthetic activity (Lintunen et al., Reference Lintunen, Paljakka, Salmon, Dewar, Riikonen and Hölttä2020). The strongest contribution of DiffPAR to carbon uptake appears under moderate soil temperatures, typically in late spring or early summer, while its effect diminishes again at the highest temperatures.

Figure 2.

SHAP dependence plots for the reduced input XGBoost model, showing how PAR, DiffPAR, NDVI, and soil temperature influence NEE predictions. Each subfigure plots a variable against its SHAP value, with color indicating a second interacting variable. The plots reveal seasonal and interaction effects, such as stronger PAR and DiffPAR contributions to carbon uptake under high NDVI or moderate soil temperatures, and that soil temperature acts more as a modulating factor than a primary driver.

Figure 2. Long description

Four scatter plots arranged in a two-by-two grid. All y-axes represent S H A P values in micromoles per second per meter squared.

Panel A, top-left. X-axis is P A R from 0 to 2000 micromoles per second per square meter. The data shows a strong negative linear trend where S H A P values decrease from positive 5 to negative 10 micromoles per second per meter squared as P A R increases. Points are color-coded by N D V I from 0.5 purple to 0.8 yellow, showing higher N D V I values correlate with more negative S H A P values at high P A R.

Panel B, top-right. X-axis is Diff P A R from 0 to 1000 micromoles per second per square meter. A negative trend is visible, with S H A P values dropping from 3 to negative 6 micromoles per second per square meter. Points are color-coded by Soil Temp from negative 10 degrees Celsius purple to 20 degrees Celsius yellow.

Panel C, bottom-left. X-axis is N D V I unitless from 0.5 to 0.8. The data remains relatively flat near a S H A P value of 2 micromoles per second per square meter until N D V I reaches 0.75, where it drops sharply toward negative 5 micromoles per second per square meter. Points are color-coded by P A R from 500 purple to 2000 yellow, in units of micromoles per second per square meter.

Panel D, bottom-right. X-axis is Soil Temp from negative 10 to 20 degrees Celsius. The data forms a V-shape with a trough near 8 degrees Celsius where S H A P values reach negative 2, then rise back toward positive 4 micromoles per second per square meter at higher temperatures. Points are color-coded by Diff P A R from 250 purple to 1000 yellow, in units of micromoles per second per square meter.

In Figure 2c NDVI’s SHAP is shown against NDVI with PAR coloring. This plot integrates the effect of station and seasonality on NEE. The high NDVI region corresponds to the peak growing season at the southern stations, when high radiation greatly increases photosynthesis, whereas during autumn/spring, high PAR corresponds to high carbon emissions. As high NDVI corresponds to plants being able to absorb radiation (Balzarolo et al., Reference Balzarolo, Vicca, Nguy-Robertson, Bonal, Elbers, Fu, Grünwald, Horemans, Papale, Peñuelas, Suyker and Veroustraete2016), this is also a logical conclusion from the models. Lower NDVI values are also more common in Värriö, which has systematically more positive NEE, which might also help to explain why lower NDVI values have higher SHAP values.

Finally, Figure 2d shows soil temperature and its SHAP values. Among the four input variables in the reduced input set, soil temperature has the lowest impact, suggesting that the ML models treat it more as a supporting factor rather than a direct driver of NEE. The seasonal pattern is visible here: as soil temperature rises to around 7°C, the carbon sink increases. However, at higher temperatures, the effect reverses, and further increases in soil temperature lead to positive SHAP contributions for NEE, likely due to stronger temperature-dependent respiration that pushes NEE toward higher values. The color shading represents diffuse PAR levels, showing that diffuse PAR has a positive effect on the carbon sink at higher soil temperatures, but a negative effect when the soil is cold.

While these plots are insightful and help validate that the models behave intuitively and reproduce the expected ecology, they do not contain much new information. To gain more insight from SHAP values, we next employ plots with a different approach.

3.4. DiffPAR SHAP values

Diffuse radiation was chosen as an input variable of interest for more detailed analysis due to open questions about the effects of diffuse radiation on forest carbon sinks, such as whether it is more mediated through clouds or aerosol (van Diepen et al., Reference van Diepen, Kaiser, Hartogensis, Graf, de Arellano and Moene2025).

SHAP values can reveal context-dependent interactions between diffuse radiation, NEE, and other variables. This motivates moving from SHAP as a model validation step to using SHAP to stratify the data. Here we use DiffPAR SHAP values to form groups and compare them with simple DiffPAR thresholding; the results suggest a clearer diffuse-light (cloud-related) enhancement of carbon uptake, motivating our climate-relevant discussion of how cloudiness and aerosol-driven scattering may influence NEE.

Figures 3 and 4 summarize how the model attributes diffuse radiation effects and how this relates to observed NEE. Figure 3 plots DiffPAR against its SHAP values with points colored by co-variates, using only observations with NDVI >0.8 to reduce the influence of seasonal variation. In Figure 3a (colored by diffuse fraction, a variable not included in the model), within intermediate DiffPAR (approximately 100–500 $ \mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ ), SHAP contributions become increasingly negative as diffuse fraction increases, consistent with a stronger carbon sink under high diffuse fraction indicating overcast skies. This can be due to differences in the spectral distribution of diffuse radiation between rather clear and overcast skies, that is, the distribution of incoming radiation across wavelengths can differ between these two cases. Moreover, Figure 3a shows a distinct pattern around DiffPAR SHAP $ \approx -2 $ : when DiffPAR SHAP values are $ \ge -2 $ , the separation between overcast and clear-sky cases becomes especially clear, whereas values $ <-2 $ are associated with patchy cloud conditions. Diffuse radiation is particularly enhanced under patchy clouds, which are known to exert the strongest positive effect on photosynthesis. For this reason, we selected $ -2 $ as the threshold in Figure 4, where we contrast high- and low-diffuse-radiation cases.

Figure 3.

Alternative view of the DiffPAR–SHAP relationship (as in subfigure 2b), restricted to observations with NDVI $ >0.8 $ to reduce the influence of seasonal variation, and colored by diffuse fraction (3a) and soil temperature (3b). Subfigure 3a confirms that points with high DiffPAR and low PAR correspond to overcast conditions (high diffuse fraction). Subfigure 3b shows no clear pattern with soil temperature, indicating that the observed diffuse radiation effect on NEE is not driven by temperature.

Figure 3. Long description

The figure consists of two vertically stacked scatter plots, labeled A and B, both sharing the same axes. The x-axis represents Diff P A R in micromoles per second per square meter, ranging from 0 to 1000. The y-axis represents Diff P A R S H A P for N E E in micromoles per second per square meter, ranging from negative 6 to 4 micromoles per second per square meter. Both plots show a strong negative linear correlation where S H A P values decrease as Diff P A R increases.

* Panel A is titled Color Diffuse Fraction. The data points are colored based on a gradient from dark purple (0.0) to bright yellow (1.0). Points with low Diff P A R (0 to 250 micromoles per second per square meter) and high S H A P values (0 to 4 micromoles per second per square meter) are predominantly yellow, indicating a high diffuse fraction. As Diff P A R increases, the points transition to orange and then purple.

* Panel B is titled Color Soil Temp. The data points are colored based on a gradient from dark purple (10 degrees Celsius) to bright yellow (25 degrees Celsius). Unlike panel A, the colors in panel B are mixed throughout the data cloud, showing no distinct spatial pattern or correlation between soil temperature and the Diff P A R S H A P relationship.

Figure 4.

Comparison of reduced input XGBoost interpretation (4a) versus traditional data grouping (4b) for assessing the effect of DiffPAR. In 4a, data points with DiffPAR SHAP below −2 show lower NEE, supporting the model’s attribution of enhanced carbon uptake to diffuse radiation. In contrast, 4b shows a weaker relationship when NEE is grouped solely by raw DiffPAR thresholds, highlighting how SHAP reveals interaction-driven effects not apparent in direct data plots.

Figure 4. Long description

Two vertically stacked boxplot graphs labeled A and B. Both share a Y-axis of N E E in micromoles per second per square meter, ranging from negative 40 to 20, and an X-axis of P A R Interval in micromoles per second per square meter, divided into nine bins from 0 to 2250.

Panel A is titled N E E by P A R interval and Diff P A R S H A P negative 2 micromoles per second per square meter threshold. It compares two groups: Diff P A R S H A P greater than or equal to negative 2 micromoles per second per square meter in yellow and Diff P A R S H A P less than negative 2 micromoles per second per square meter in purple. As P A R increases, the purple boxes consistently show lower median N E E values compared to the yellow boxes, particularly in the 500 to 1500 micromoles per second per square meter range. Sample sizes n are noted below each box.

Panel B is titled N E E by P A R interval and Diff P A R 400 micromoles per second per square meter threshold. It compares Diff P A R less than 400 micromoles per second per square meter in blue and Diff P A R greater than or equal to 400 micromoles per second per square meter in red. The separation between the blue and red boxes is less pronounced than the separation in Panel A, with significant overlap in the interquartile ranges across most P A R intervals.

Both panels show a general trend of decreasing N E E as P A R increases, with outliers represented by individual circles above and below the whiskers.

In Figure 3b (colored by soil temperature), no comparable gradient is evident, indicating that the diffuse-radiation signal is not a proxy for thermal conditions and is unlikely to be explained by collinearity with temperature. These patterns pool observations across all sites, reinforcing the spatial robustness of the effect.

Figure 4 contrasts the DiffPAR SHAP grouping of $ -2 $ with raw DiffPAR grouping of 400 $ \mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ , clarifying how diffuse radiation influences NEE through interactions rather than as a simple monotonic response. Importantly, both the SHAP-based and raw DiffPAR-based groupings contain a similar number of observations (i.e., roughly balanced $ n $ ), so the sharper separation in 4A, is not an artifact of uneven sample sizes. When defining the raw DiffPAR grouping, we selected a threshold of 400 $ \mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ . This choice was motivated by Supplementary Figure A4C, where DiffPAR values separated by the $ \mathrm{DiffPAR}\ \mathrm{SHAP}<-2 $ threshold show division where, in the $ \mathrm{DiffPAR}\ \mathrm{SHAP}<-2 $ group, 400 lies below the lower quartile, while in the $ \mathrm{DiffPAR}\ \mathrm{SHAP}\ge -2 $ group it lies above the upper quartile. This makes it a fitting boundary between conditions dominated by strong diffuse-radiation effects and those without. Moreover, a comparable threshold of approximately 200 W, m $ {}^{-2} $ of diffuse radiation (equivalent to roughly 400 $ \mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ PAR) was identified in our recent related work (Ezhova et al., Reference Ezhova, Aarne, Arola, Lipponen, Lintunen, Kokkola, Ylivinkka, Yli-Juuti, Petäjä, Kerminen, Virtanen and Kulmala2025b) as the point where diffuse radiation begins to have a beneficial effect on photosynthesis.

In Figure 4a, observations with $ \mathrm{DiffPAR}\ \mathrm{SHAP}<-2 $ , that is, cases where the model attributes a strong diffuse-radiation contribution, exhibit markedly lower (more negative) NEE on average; these points also have high DiffPAR (mean $ =571.75\mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ ), consistent with an enhanced carbon sink under diffuse-light conditions. By comparison, Figure 4b groups observations by raw DiffPAR thresholds alone and shows a much weaker separation in NEE, indicating that total DiffPAR without context is insufficient to isolate the effect. The stronger contrast in 4A implies that the model is able to disentangle the diffuse radiation effect from other variables; therefore, the model suggests that the diffuse radiation effect is larger compared to standard threshold-based analysis.

Overall, after excluding the lowest and highest PAR bins (0–250 and 2000–2250), the SHAP-stratified subset ( $ \mathrm{DiffPAR}\ \mathrm{SHAP}<-2 $ ) has a median NEE that is $ 1.545\mu \mathrm{mol},{\mathrm{m}}^{-2},{\mathrm{s}}^{-1} $ more negative than the subset defined by raw $ \mathrm{DiffPAR}\ge 400 $ , indicating that the model-attributed signal isolates a stronger diffuse-light enhancement than raw DiffPAR alone.

Neimane-Šroma et al. (Reference Neimane-Šroma, Durand, Lintunen, Aalto and Robson2024) recently studied the diffuse radiation fertilization effect in a boreal Scots pine forest (SMEAR II forest) at ecosystem, branch, and ground vegetation scales, demonstrating that at equivalent total PAR, diffuse radiation conditions yield higher ecosystem productivity by enhancing photosynthesis across scales more than direct radiation conditions. Our SHAP-based analysis is novel in independently identifying this diffuse radiation fertilization as a distinct pattern in the model’s predictions, pointing to a previously unreported phenomenon from a modeling perspective. This explainable ML approach, therefore, complements and extends the physiological insights of Neimane-Šroma et al. (Reference Neimane-Šroma, Durand, Lintunen, Aalto and Robson2024), bridging their multi-scale gas-exchange findings with ecosystem-scale model behavior.

A plausible explanation is that under overcast (highly diffuse) skies, light penetrates deeper into the canopy and alleviates direct-radiation stress on sun-exposed foliage, giving the forest a net photosynthetic boost not offset by lower total PAR.

Since SHAP reveals the model’s internal logic rather than proving causality, targeted branch- and forest floor-level gas-exchange measurements during high-diffuse, low-PAR periods (as in Neimane-Šroma et al., Reference Neimane-Šroma, Durand, Lintunen, Aalto and Robson2024) are needed to verify the mechanism. If confirmed, this result implies that projected increases in cloudiness or aerosol-driven light scattering could make boreal forests stronger short-term CO₂ sinks than current models assume, challenging the classical assumption that reduced PAR necessarily limits carbon uptake.