Forests for forests: combining vegetation indices with solar-induced chlorophyll fluorescence in random forest models improves gross primary productivity prediction in the boreal forest

Remote sensing is a powerful tool for understanding and scaling measurements of plant carbon uptake via photosynthesis, gross primary productivity (GPP), across space and time. The success of remote sensing measurements can be attributed to their ability to capture valuable information on plant structure (physical) and function (physiological), both of which impact GPP. However, no single remote sensing measure provides a universal constraint on GPP and the relationships between remote sensing measurements and GPP are often site specific, thereby limiting broader usefulness and neglecting important nuances in these signals. Improvements must be made in how we connect remotely sensed measurements to GPP, particularly in boreal ecosystems which have been traditionally challenging to study with remote sensing. In this paper we improve GPP prediction by using random forest models as a quantitative framework that incorporates physical and physiological information provided by solar-induced fluorescence (SIF) and vegetation indices (VIs). We analyze 2.5 years of tower-based remote sensing data (SIF and VIs) across two field locations at the northern and southern ends of the North American boreal forest. We find (a) remotely sensed products contain information relevant for understanding GPP dynamics, (b) random forest models capture quantitative SIF, GPP, and light availability relationships, and (c) combining SIF and VIs in a random forest model outperforms traditional parameterizations of GPP based on SIF alone. Our new method for predicting GPP based on SIF and VIs improves our ability to quantify terrestrial carbon exchange in boreal ecosystems and has the potential for applications in other biomes.


Introduction
Uncertainty in future climate projections is largely driven by terrestrial ecosystem feedbacks on the carbon cycle (Friedlingstein et al 2014). A major contributor to terrestrial ecosystem feedbacks is plant carbon uptake via photosynthesis. Plant carbon uptake can be estimated locally at the tower/site level, as gross primary productivity (GPP), but remote sensing is necessary to scale and understand carbon exchange across space and time (Anav et al 2015). This is especially relevant in arctic-boreal ecosystems which play a major, but highly uncertain, role in the global carbon cycle (Bonan 2008, Thurner where PAR is the photosynthetically active radiation, f PAR is the fraction of photosynthetically active radiation absorbed by chlorophyll, and LUE P is the light-use-efficiency of photosynthesis. Proposed remote sensing metrics for approximating GPP typically track either the physical (f PAR) or the physiological (LUE P ) components of this equation. Solar-induced chlorophyll fluorescence is an especially powerful remote sensing metric for understanding GPP over a variety of ecosystems (Sun et al , 2018 and across a variety of scales (Yang et al 2017, Mohammed et al 2019, Porcar-Castell et al 2021, Pierrat et al 2022 due to its connection to both the physical (PAR × f PAR) and the physiological (LUE P ) components of GPP. Solarinduced fluorescence (SIF) can be described using the light-use-efficiency model: where LUE F is the light-use-efficiency of fluorescence, and f esc is the fraction of emitted SIF photons which escape the canopy. SIF and GPP share the physical drivers PAR × f PAR, which largely explains the relationship between SIF and GPP in cropping systems (Dechant et al 2020). Equations (1) and (2) can be combined to eliminate the absorbed photosynthetically active radiation (APAR = PAR × f PAR): The physiological component of SIF, LUE F covaries with LUE P under moderate light conditions (Porcar-Castell et al 2014 and averaged over broad spatio-temporal scales , Pierrat et al 2022. This is the basis for SIF as a proxy for GPP when measured from space across a broad range of ecosystems (Sun et al , 2018. However, these conditions are not necessarily always met (Marrs et al 2020), especially in the boreal forest where complex canopy structure and divergence between the light and carbon fixation reactions of photosynthesis can complicate these signals (Maguire et al 2020, Pierrat et al 2022. Therefore, a more nuanced relationship between SIF and GPP that takes into consideration divergence between variations in LUE F , LUE P , and f esc will further improve the utility of SIF as a proxy for GPP. Additional reflectance-based remotely sensed metrics, i.e. vegetation indices (VIs), can provide more information on both the structural and physiological processes impacting GPP and can therefore improve our ability to track and understand GPP with remote sensing.
Physically, SIF and GPP are impacted by illumination conditions and canopy structure (i.e. the combined effects of leaf-area index (LAI), vertical distribution of LAI, leaf orientation, clumping, etc) which mediate both f PAR and f esc . At a constant PAR and LUE P , GPP amplifies under cloudy sky conditions because a higher diffuse fraction allows light to penetrate deeper into the canopy, thus increasing f PAR (Gu et al 2002, Alton et al 2007, Durand et al 2021. Although the amplification of GPP under diffuse skies has been well documented, it is often not considered when approximating GPP with remote sensing proxies. Canopy structure is often approximated using greenness based VIs which are sensitive to chlorophyll content in an instrument field of view. Thus, they are generally a good approximation for f PAR. The normalized difference vegetation index (NDVI) effectively tracks vegetation productivity in ecosystems where chlorophyll content and carbon uptake are closely correlated (i.e. larger variations in f PAR than LUE P ) (Tucker 1979, Yang et al 2017, Wang and Friedl 2019. This is not the case in the boreal forest where changes in carbon uptake do not correlate with changes in chlorophyll content (Sims et al 2006, Garbulsky et al 2010, Gamon et al 2013. Additionally, NDVI is extremely sensitive to the presence of snow cover (Pierrat et al 2021a(Pierrat et al , 2022. The nearinfrared reflectance from vegetation (NIRv) expands on NDVI by multiplying NDVI by the total scene NIR reflectance, thereby amplifying the vegetated signal. NIRv has shown stronger correlations with f PAR and GPP than NDVI alone (Badgley et al 2017. NIRv can also be used as a proxy for f esc (Zeng et al 2019). Both NDVI and NIRv provide useful information on the physical/structural influences on SIF and GPP, however, it is unclear how sensitive these indices are to changes in f PAR and f esc and if there exists a universal relationship across ecosystems.
Physiologically, the light-use-efficiencies of both fluorescence and photosynthesis are mediated by non-photochemical quenching-a heat dissipation mechanism that plants utilize to avoid damage from excess sunlight (Raczka et al 2019, Walter-McNeill et al 2021. The extent to which non-photochemical quenching is necessary to protect plant tissue depends on a host of environmental controls that determine a plant's ability to photosynthesize (Demmig-Adams and Adams 2006). Of particular relevance in this respect is plant's sensitivity to temperature. Freezing over winter or heat waves over summer create stress conditions which subsequently impact LUE P . Monitoring the extent of non-photochemical quenching thus provides insight into LUE P (Adams et al 2004), and likely LUE F , although those dependencies are not yet well quantified.
Non-photochemical quenching manifests in two forms of photoprotection in boreal ecosystems, both of which can be detected with remote sensing measures (Gamon et al 1997, 2016, Demmig-Adams et al 2008. The first is rapidly reversible and therefore important under short-term high light stress. It can be detected using the photochemical reflectance index (PRI). Therefore, PRI tracks changes in LUE P over the diurnal cycle (Gamon et al 1997, Wong and Gamon 2015a, 2015b. The second form is sustained and thus important over longer time periods such as winter in regions where freezing temperatures limit the ability for plants to photosynthesize (Verhoeven 2014). This form can be tracked using the chlorophyllcarotenoid index (CCI) (Gamon et al 2016). CCI has successfully tracked variations in GPP over the seasonal cycle because it is sensitive to long-term changes in LUE P .
Including information on heat dissipation dynamics through remotely sensed products (PRI and CCI) improves modeling of photosynthetic phenology (Wong et al 2022), especially when used in conjunction with SIF (Wang et al 2020, Hikosaka andTsujimoto 2021). In addition, correcting the SIF signal for canopy and structural effects using structurally sensitive metrics also improves its relationship with GPP , Lu et al 2020. It has been well documented and mechanistically explained why VIs are sensitive to changes in non-photochemical quenching or canopy structure but there are no universal quantitative relationships among them. To develop such a relationship, the impacts from canopy structure, view/sun angle effects, snowcover, and potential differences in instrumentation (field of view, instrument sensitivity) must all be accounted for. Until mechanistic models are able to effectively account for all the aforementioned effects, the capacity for these indices to inform SIF and GPP is limited to qualitative or site-specific empirical relationships. In order to make more effective use of remote sensing as a proxy for GPP, it is necessary to have a quantitative framework that can relate SIF, VIs, and GPP.
Advances in machine learning have provided exciting opportunities for data analysis and predictive modeling. In particular, random forest models are non-parametric in nature and are therefore well suited for approximating non-linear, multi-parameter relationships in complex systems (Breiman 2001). Because random forest models do not assume functional dependencies between input variables and predicted output and rather 'learn' relationships based on input data, they are particularly useful for systems where statistical parameterized models either do not exist or are not well constrained. Random forest models are also highly interpretable compared with other machine learning techniques due to predictor importance estimates, which make them more useful for explaining and understanding observed relationships. Finally, random forest models have already been used to understand nuance in the SIF-GPP relationship (Jiao et al 2019, Bai et al 2022, Pierrat et al 2022. We propose that using random forest models as a tool to understand and predict GPP based on a combination of remote sensing metrics will present a significant improvement over traditional parameterized models (Monteith 1972, Dechant et al 2020 or other machine learning approaches. The central question of this study is therefore: can we improve our ability to predict GPP from SIF by using random forest models as a quantitative framework that can incorporate additional physical and physiological information provided by VIs? To answer this question, we present 2.5 years of tower-based remote sensing data across two boreal forest locations, qualitatively evaluate the seasonal and diurnal trends among them, and present an interpretation on the physical and physiological information contained in them (section 3.1). We justify the use of random forest models as a tool for understanding SIF, VI, and GPP dynamics by showing that they can accurately reproduce SIF-GPP-PAR relationships (section 3.2). Finally, we explore the utility of random forest models and remotely sensed products for improving GPP estimation by showing how random forest models driven by SIF, VIs, and temperature improve the prediction of GPP and shed light on important physical and physiological processes (sections 3.3 and 3.4).

Site descriptions: Southern Old Black Spruce (SOBS) and National Ecological Observatory Network (NEON) Delta Junction (DEJU)
We collected data at the SOBS site (FLUXNET site code CA-Obs) and the National Ecological Observatory Network (NEON), DEJU which represent the northern and southern limits of the North American boreal forest and associated environmental conditions (figure 1, table 1). SOBS is located near the southern limit of the boreal forest in Saskatchewan, Canada (53.98 • N, 105.12 • W) (Jarvis et al 1997). It is a mixed forest stand with stem density predominantly (90%) black spruce (Picea mariana), and scattered

Data collection: tower-based remote sensing, GPP, and environmental variables
We collected tower-based remotely sensed measurements (far-red SIF, NDVI, NIRv, CCI, PRI) using PhotoSpec (see Grossmann et al 2018) for detailed instrument description) at both the SOBS and DEJU field sites. Measurements ran from August 2019 to December 2021 at both field locations (figures 1 and 2). At both sites, Photospec was installed atop the scaffolding tower facing due north. It has a narrow field of view (0.7 • ), 2D scanning capabilities, and simultaneously measures SIF and VIs at the same point in the canopy (Grossmann et al 2018). Individual measurements take approximately 20 s. We took canopy representative scans at both field locations within a 30 min window and averaged measurements together to compare with the temporal resolution of GPP and environmental variables. SIF was retrieved in the far-red (745-758 nm) wavelength range using a Fraunhofer-line based fitting algorithm (Grossmann et al 2018). The Fraunhofer-line based approach makes SIF retrievals insensitive to atmospheric scattering and therefore robust even under cloudy sky conditions (Frankenberg et al 2011, Mohammed et al 2019, Chang et al 2020. We filtered data for low quality retrievals and retrievals with unstable sky conditions (Pierrat et al 2021a(Pierrat et al , 2022 for both sites. The VIs were calculated as follows, with ρ nm:nm = the average reflectance across a wavelength range in nm: In boreal ecosystems, the onset of photosynthesis often occurs prior to complete snowmelt (Starr andOberbauer 2003, Parazoo et al 2018). Therefore, any Table 1. Summary of methods for measurement and data processing at the two field locations.

Measurement
SOBS DEJU SIF, VIs, PAR, and Df PhotoSpec mounted atop the tower and processed following (Grossmann et al 2018, Pierrat et al 2022.
PhotoSpec mounted atop the tower and processed following (Grossmann et al 2018, Pierrat et al 2022. GPP

Eddy-covariance
Taken using a 3D sonic anemometer (CSAT3, Campbell Scientific, Logan, UT) in combination with a closed-path infrared gas (CO2/H2O) analyzer (LI-7200 analyzer, Li-Cor, Lincoln, NE) operated in absolute mode. We performed quality assurance on the data using the standard Fluxnet-Canada method following , Barr et al 2004.
Obtained from (National Ecological Observatory Network (NEON) 2022a) using a Campbell Scientific CSAT-3 3D Sonic Anemometer and LI-COR-LI7200 gas analyzer. We performed quality assurance on carbon fluxes based on turbulent and storage fluxes separately, using a bivariate statistical procedure for each, to overcome quality flag restrictions in the 'expanded' NEON eddy-covariance bundle. The 3% outliers (3% of rarest events from the tails of each distribution) were excluded from joint probability distributions for all available data for (a) turbulent flux and PAR, and separately for ( proposed approach for predicting GPP based on remote metrics must be effective even in the presence of snow. Because of this, we did not filter for snow but visually identified snow dates using PhenoCam imagery . Filtering for snow using an NDVI threshold >0.5 (Magney et al 2019, Cheng et al 2020) did not change the results of this study (figures S1-S4).
To include impacts of illumination conditions on our analysis, we determined a measure of direct vs diffuse radiation (Df). Df reflects the deviation of PAR at a given solar zenith angle from the expected PAR during a clear sky reference day so that Df = 1 is clear sky conditions (Pierrat et al 2021b). Df values were calculated for every PhotoSpec measurement (∼20 s resolution) and averaged together in 30 min windows to compare with GPP and environmental measurements (table 1). Values where Df > 0.8 are considered clear sky conditions. Data collection and processing for eddy covariance and meteorological data for both SOBS and DEJU are summarized in table 1.

Data analysis: random forest and parameterized models
We trained and tested a variety of random forest models (table 2). All random forest models were produced using Matlab's TreeBagger function (MATLAB 2019) which is based on the random forest algorithm from (Breiman 2001). We used the full data collection window (August 2019-December 2021, figure 2) to train all models. All models were created using 100 regression trees and sampled with replacement on an in bag fraction of 0.7 (i.e. 70% of the data were randomly chosen to train the models saving 30% to test the models). Out-of-bag (OOB) predictor importance estimates were determined using the permuted predictor delta error following the standard Classification And Regression Tree (CART) algorithm (Breiman 2001) using the remaining 30% of data. Model performance was evaluated by calculating the Pearson's correlation coefficient (Dickinson Gibbons et al 1985) between measurements and predictions on the full dataset (R 2 ) and the reserved test dataset (OOB R 2 score). Specific models are identified in section 3 with the naming conventions in table 2. We compared random forest models with common parameterized light-use-efficiency models for the relationships among SIF, GPP, and PAR (equations (1)-(3)). All curve fitting and goodness of fit statistics (R 2 values) were done using Matlab's fit function (MATLAB 2019). Specifics of fitted equations are provided in section 3 and figure captions where relevant.

Trends among SIF, VIs, and GPP
We present the relationships and trends among SIF, VIs, and environmental parameters across seasonal and diurnal scales for both field locations in figures 2 and 3. Compared with SOBS, DEJU generally exhibits a shorter growing season, a smaller summer maximum GPP, and lower values for all remotely sensed metrics (figures 2 and 3).
SIF tracks the seasonal cycle and daily-weekly variability in GPP (figure 2). Monthly diurnal profiles of SIF show good agreement with monthly diurnal profiles of GPP across both sites (figure 3 SIF increases in spring prior to changes in GPP (figure 2) and shows a small diurnal profile over winter at both sites, while GPP does not (figure 3). This early spring increase has been reported in other evergreen locations (Magney et al 2019, Pierrat et al 2022, Yang et al 2022 and is attributed to persistent photosystem II activity (Porcar-Castell 2011, Bowling et al 2018. This leads to a winter light response of SIF which means SIF increases coincident with PAR in winter-/early spring in a way that does not reflect changes in GPP. This does not preclude the use of SIF as a proxy for GPP, but it must be accounted for to prevent overestimation of GPP in winter (Pierrat et al 2021a(Pierrat et al , 2022. Structurally sensitive indices (NDVI and NIRv) show little variability over summer and are sensitive to snow in winter across both sites (figure 2). Over summer at SOBS, NDVI and NIRv show greater changes than at DEJU, which reflects the fact that SOBS is a mixed-species forest with scattered deciduous larch trees. SOBS therefore experiences greater changes in canopy structure over the seasonal cycle than evergreen DEJU. NDVI peaks in mid-morning and mid-afternoon over the summer months at SOBS and has no clear pattern at DEJU (figure 3). NIRv is slightly higher in the morning at DEJU and has no clear pattern at SOBS (figure 3). These patterns suggest that NDVI and NIRv are useful for accounting for shifts in viewing geometry and illumination effects across the diurnal cycle.
Physiologically sensitive indices (CCI and PRI) generally track the seasonal cycle of GPP but show sensitivity to snow cover in winter. CCI is more sensitive to snow at DEJU due to a higher fraction of snow cover on the canopy. Both CCI and PRI increase in spring following GPP at both sites (figure 2). This is consistent with a two-phased spring transition beginning with the onset of photosynthesis, followed by a reduction in sustained photoprotection (Pierrat et al 2021a). CCI at DEJU is lower than at SOBS which may indicate differences in chlorophyll:carotenoid ratios between the sites. Diurnally, PRI shows a small decrease around midday during summer months at both sites (figure 3). This highlights the sensitivity of PRI to rapidly reversible photoprotection which will occur under high-light, high-stress conditions. CCI shows a slight diurnal pattern at SOBS and does not show a consistent diurnal pattern at DEJU. The lack of a diurnal pattern at DEJU supports the idea that CCI is sensitive to sustained cold-season photoprotection which will not change over the course of the diurnal cycle. PRI and CCI thus provide information on plant heat dissipation processes at both seasonal and diurnal timescales.
Environmental conditions at the two sites (PAR, Df, and air temperature) are consistent with expected environmental patterns. PAR increases prior to an increase in GPP and begins to decrease in fall prior to a decline in GPP at both sites. Df does not show a defined seasonal cycle, nor does it show a prominent diurnal pattern at either site. Thus, diffuse sky conditions are largely independent of season or time of day. Air temperature tracks the seasonal cycle of GPP, and air temperatures above 0 • C are a good indicator of growing season length (Parazoo et al 2018, Pierrat et al 2021a, Stettz et al 2022. Air temperature peaks in the afternoon following the peak in SIF and GPP.

SIF and GPP light response curves: light-use-efficiency and random forest models
We tested the ability of random forest models to reproduce quantitative SIF-GPP dynamics by comparing light response curves following the parameterized light-use-efficiency model (figure 4 row (a)) with light response curves produced by environmentally driven random forest models (table 2, ENV-SIF and ENV-GPP, figure 4 row (b)), at a half-hourly resolution at the SOBS site. The same analysis tested at the DEJU site shows consistent results (figures S5 and S6).
The half-hourly data and parameterized lightuse-efficiency model (GPP = GPPmax×PAR c×PAR , (Michaelis andMenten 1913, Monteith 1972) and (SIF =c × PAR, equation (2) show a curved light response for GPP, consistent with the light saturation of GPP, and a near linear light response for SIF ( figure 4 row (a)). Both GPP and SIF show a seasonally dynamic light  (Michaelis andMenten 1913, Monteith 1972) and SIF =c × PAR (equation (2)). Row (b) shows light response curves of half-hourly GPP and SIF produced from two random forest models (ENV-GPP and ENV-SIF). response with minima over winter and maxima over summer. The light response of SIF approaches but does not go to zero over the winter months due to persistent winter photosystem II activity (Porcar-Castell 2011, Bowling et al 2018, Yang et al 2022. Both ENV-SIF and ENV-GPP models had strong performance with high R 2 (R 2 = 0.93 and 0.94 respectively) and OOB R 2 scores (OOB R 2 = 0.86 and 0.89 respectively) ( figure S7). This highlights the ability of random forest models to capture the environmental dependencies of both SIF and GPP. The two trained models were run for each month using the range of PAR values for that month and the monthly average for the rest of the environmental predictor variables to reproduce the light response curves of SIF and GPP ( figure 4 row (b)). The monthly light response curves produced by the random forest models have the same patterns as both the data and the parameterized light-use-efficiency models. The light response of GPP is curved and the light response of SIF is largely linear. SIF shows a slight curvature at high PAR values which we attribute to changes in f esc at high PAR (Pierrat et al 2022). Both GPP and SIF exhibit a monthly variable light response and the light response of SIF does not drop to zero over winter, consistent with the light-use-efficiency model.
These results highlight the efficacy of random forest models for reproducing quantitative relationships and environmental dependencies for GPP and SIF without prescribing a parameterized model onto the data (Chen et al 2021). Because of this, and the additional information that can be provided by additional remote sensing metrics (section 3.1), we justify the use of random forest models for GPP prediction.

Random forest models compared with parameterized models for predicting GPP
We tested the ability of random forest models to improve GPP prediction by testing traditional parameterized models (a linear SIF-GPP relationship, a non-linear SIF-GPP relationship that takes into account the light saturation of GPP (Monteith 1972) and a monthly variable non-linear SIF-GPP relationship (Damm et al 2015, Pierrat et al 2022 against our new random forest approach combining multiple remote sensing indices for predicting GPP (figures 5 and 6). The two random forest models used to test this (table 2, RS-SOBS and RS-DEJU) were driven by remote sensing and meteorological variables that reflect the physical and physiological controls on photosynthesis relevant at a half-hourly temporal resolution.
The linear fit between SIF and GPP (figure 5 row (a)) performs moderately-poorly with an R 2 = 0.58 at SOBS and R 2 = 0.43 at DEJU ( figure 6 row (a)). The non-linear fit based on (Damm et al 2015) shows little to no improvement from the linear fit with R 2 = 0.60 at SOBS and R 2 = 0.43 at DEJU (figure 5 row (b)). Both models with a fixed SIF-GPP relationship overestimate GPP in winter due to a persistent winter light response of SIF. The monthly variable non-linear SIF-GPP relationship (figure 5 row (c)) shows a marked improvement in the predictability of GPP at both the SOBS and DEJU sites with R 2 = 0.76 and R 2 = 0.66 respectively ( figure 6 row (c)). In addition to improved R 2 values, the monthly variable nonlinear SIF-GPP relationship helps account for the persistent winter light response and the winter overestimation of GPP is no longer observed. The random forest models driven by remotely sensed products (figure 5 row (d)) improved the predictability of GPP compared to all parameterized models at both sites with R 2 = 0.90 and 0.86 and OOB R 2 = 0.81 and 0.73 at SOBS and DEJU respectively (figure 5 row (d)). The overestimation of GPP in winter is not observed using the random forest approach and residuals between predicted and measured GPP are more homoscedastic across seasons. , (Damm et al 2015). Row (c) shows the same non-linear fit but fitted monthly to create a monthly variable SIF-GPP relationship (Pierrat et al 2022). Row (d) shows the input variables predictor importance estimates for random forest models RS-SOBS and RS-DEJU.
The predictive power of input variables is evaluated with the predictor importance estimates (figure 5 row (d)). SIF was the most important predictor at SOBS and second most important predictor at DEJU, highlighting the power of SIF for GPP prediction. DEJU showed a higher dependence on air temperature than SOBS which could point to stronger temperature limitations at this site. Df was an important predictor across both sites which highlights the dependence of GPP on diffuse vs. direct radiative conditions (Durand et al 2021, Pierrat et al 2021b. NIRv was the third most important predictor at SOBS which could reflect the fact that SOBS is a mixed-species forest interspersed with deciduous trees and will thus have more dramatic changes in canopy structure over the course of the season than DEJU. CCI was moderately important at both SOBS and DEJU which highlights the relevance of CCI for capturing changes in sustained photoprotection over winter. NDVI was not a particularly important predictor for either SOBS or DEJU which could reflect the fact that NIRv is more effective at capturing changes in canopy structure (Badgley et al 2017 than NDVI. PRI was also not particularly important for either SOBS or DEJU which suggests that rapidly reversible non-photochemical quenching dynamics may be more effectively captured in the SIF signal when a prescribed relationship between SIF and GPP is not used. Including PAR as an input variable does not improve model results but does decrease the predictor importance of SIF (figures S8 and S9). Substituting PAR for SIF as a predictor variable also does not change model performance (figures S10 and S11). This suggests that VIs effectively capture the structural (f PAR) and physiological (LUE P ) factors relevant for predicting GPP.

Random forest models for predicting GPP across boreal forest sites
Satellite remote sensing enables the approximation of GPP over a broader spatial range than towerbased measurements, making it essential for understanding regional GPP dynamics. However, relationships between remotely sensed products and GPP are often site specific and thus require a new model or set of parameters for each ecosystem or plant functional type. We tested the potential of our random forest modeling approach across sites for a 'universal' model for GPP based on remotely sensed products. We trained a model (RS-total) to predict GPP from data that are readily accessible or can be inferred from satellite measurements (table 2, RS-total). We used data at a daily resolution averaged between 10:00 and 14:00 under all (figure 7) and only sunny (figure S12) sky conditions to replicate satellite observations across both the SOBS and DEJU sites together. SOBS and DEJU are at the latitudinal extremes of the boreal ecosystem, thus, testing the feasibility of interpolation across the boreal region.
Our results show excellent performance on the predictability of GPP based on remotely sensed metrics (R 2 = 0.94 and OOB R 2 = 0.89) (figure 7). Residuals between predicted and measured GPP are highly homoscedastic and GPP is not overestimated in winter. Predictor importance estimates show that SIF is the most valuable predictor for GPP, but air temperature, NDVI, CCI, and NIRv are all relevant predictors. To test the universality of this approach, we tested the same model but included a site flag as a predictor (either SOBS or DEJU) (figure S13). This model showed no improvement in performance (R 2 = 0.94 and OOB R 2 = 0.89) and had site flag as the least relevant predictor. This may be because importance estimates are biased towards predictors with many classes or different values (Loh and Shih 1997) and we only had two flags for the two sites. Alternatively, this, as well as the success of the model without the site flag, points to our random forest modeling approach being independent of site, and therefore potentially generalizable across the boreal biome at the satellite level , Li et al 2018. This approach also works whether or not the data have been filtered for snow contamination (figure S4) and shows near identical results when only clear sky days are used (figure S12) which are both advantageous in the boreal. The combined random forest modeling approach also outperformed all other parameterized models at a daily midday resolution for the two sites separately (figures S14 and S15) and worked well for the combined sites at a half-hourly resolution (figure S16). Substituting PAR for SIF also showed good model performance and similar predictor importance estimates (figure S17). This suggests that VIs are effective at capturing the physical and physiological effects relevant for predicting GPP. Our results support the use of this approach for improving the predictability of GPP from remote sensing observations because it can account for physical and physiological mechanisms impacting remotely sensed signals.

Conclusions
In this paper, we present a quantitative framework for predicting GPP using random forest models driven by a set of remotely sensed products and air temperature. Multiple years of tower-based remote sensing and GPP data across two field locations at the northern and southern ends of the North American boreal forest show that VIs and SIF contain valuable information on the physical and physiological drivers of GPP. Additionally, random forest models driven by environmental variables are able to reproduce light response curves of SIF and GPP and are thus able to capture quantitative relationships of plant physiology. These results justify the use of random forest models to predict GPP based on a set of random forest parameters. Random forest models outperform traditional parameterized models based on SIF alone for predicting GPP because they are able to incorporate physical and physiological information provided by additional remote sensing metrics without prescribing a parameterized model. This approach is not site specific and therefore has the potential to be scaled across the boreal domain using satellite measurements. Finally, this approach has potential for use in other ecosystems where remote sensing data is available. Random forest models improve the utility of SIF and vegetation index data for scaling and understanding GPP across space and time and thus present an exciting opportunity to better understand vegetation's role in the global carbon cycle.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: 10.5281/zenodo.7231157.
The meteorological, soil and eddy-covariance measurements at SOBS were made with support from the Global Institute for Water Security, University of Saskatchewan.
Data from DEJU were obtained from The National Ecological Observatory Network. The National Ecological Observatory Network is a program sponsored by the National Science Foundation and operated under cooperative agreement by Battelle. This material is based in part upon work supported by the National Science Foundation through the NEON Program. Some data used in this research were obtained through the NEON Assignable Assets program.

Acknowledgments
We would like to thank the reviewers and editor for their thoughtful comments which improved the quality of our manuscript. We would also like to acknowledge the support provided by NEON staff for their help in the collection of data used in this study.

Funding
This work was supported by NASA's Earth Science Division IDS (Awards 80NSSC17K0108 at UCLA, 80NSSC17K0110 at JPL) and ABoVE programs (Award 80NSSC19M0130). A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This material is also based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1650604 and DGE-2034835. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation. DRB and TSM were funded during this time by the Macrosystems Biology and NEON-Enabled Science program at NSF (Award 1926090). We acknowledge funding by the Canadian Space Agency and Natural Sciences and Engineering Research Council of Canada (NSERC).