Importance of the legacy effect for assessing spatiotemporal correspondence between interannual tree-ring width and remote sensing products in the Sierra Nevada

Carbon uptake and tree growth are important factors for assessing productivity and long-term carbon storage. Measurements of radial stem growth are mainly performed at the individual tree scale and can be used to infer ecosystem net primary productivity (NPP). However, these measurements are spatially limited, and remote sensing provides a promising tool to track vegetation function and productivity across spatial scales, making it a viable technique for assessing variation in interannual tree growth and carbon storage. In this study we examined the correspondence between in-situ annual tree-ring width across four dominant evergreen species in the Sierra Nevada and a wide range of remotely sensed products linked to carbon uptake including NPP, gross primary productivity (GPP), net photosynthesis (PSN), normalized difference vegetation index (NDVI), near infrared reflectance of vegetation index (NIR V ), and chlorophyll/carotenoid index (CCI) from MODIS (Moderate resolution Imaging Spectroradiometer) as well as downscaled solar-induced fluorescence (SIF) products, using a 14 year dataset (2000 – 2014) across 62 forest sites. We show that variation of annual tree-ring width was best captured by the annual sum of MODIS GPP, with a legacy effect (5-month backwards shift). Across all forest sites, MODIS GPP with a 5-month legacy effect showed moderate correspondence with tree-ring width ( r = 0.60). Within each individual site, however, the interannual correspondence between MODIS GPP with a legacy effect and tree growth was stronger (median r = 0.70 vs 0.14 without a legacy effect). The importance of legacy effects in explaining tree growth variation within sites indicates that tree growth each year is influenced by carbon uptake during the latter part of the previous growing season. Additional local environmental factors also explained annual variation in tree-ring width, including (in descending order of importance) local tree density, latitude, slope, DBH, elevation and aspect.


Introduction
Carbon uptake and the allocation of resources to growth have important roles in the carbon cycle and are strongly influenced by climate variability (Babst et al., 2013;Messori et al., 2019).Quantifying net primary productivity (NPP) is a challenge as it incorporates both carbon uptake (net photosynthesis) and biomass growth including foliage, wood and roots (Litton et al., 2007;Luyssaert et al., 2010).Furthermore, climate change exacerbates uncertainties regarding how ecosystem carbon uptake and sequestration may respond in the future (Le Quéré et al., 2009).On one hand, warming temperatures may increase productivity and growth due to a longer growing season (Jeong et al., 2011) and improved photosynthetic efficiency (Myneni et al., 1997;Nemani et al., 2003).On the other, increased frequency and severity of drought may negatively influence carbon uptake and tree growth (Buermann et al., 2018;Meehl and Tebaldi, 2004).Ultimately, climate change will have a substantial impact on tree productivity and growth, making it increasingly important to develop tools to monitor these changes across spatial and temporal scales (Babst et al., 2018).
Tree-ring width corresponds strongly with tree growth, which is an important contributor to long term carbon storage of forests (Körner, 2017;Litton et al., 2007;Luyssaert et al., 2010) and is strongly influenced by climate (Babst et al., 2013;Dolanc et al., 2013;Fritts, 1966;Walker and Johnstone, 2014).Assessing tree-ring width is useful for understanding historical climate-growth responses traditionally associated with temperature (Babst et al., 2018).However, recent trends show a reduction of sensitivity of tree ring widths to temperature, and instead annual growth may be limited by other factors such as the phenology of photosynthesis, water availability, nutrient availability, carbohydrate pools, soil frost, pathogens and herbivory (D'Arrigo et al., 2008).In addition, there may be carryover, or legacy, effects between years that can influence tree growth (Babst et al., 2014a).Thus when considering tree-ring width at the annual timescale, there may be a temporal separation, or legacy effect, between climate and tree growth (Fritts, 1966).When a legacy effect is accounted for, the relationship between climate and tree growth can improve (Gough et al., 2008;Lagergren et al., 2019).This is in part associated with storage of prior year resources that can be remobilized for growth in the following year (Royce and Barbour, 2001).
The annual variation of tree-ring width is important in the context of NPP, as wood biomass represents a large proportion of resource allocation from the annual carbon budget (Litton et al., 2007;Luyssaert et al., 2010).However, acquiring tree-ring width data requires in situ collection of tree cores or dendrometry from individual trees.Therefore, there is a demand for using large-scale techniques to monitor the variation of annual tree-ring width (Babst et al., 2018).At the site scale, the eddy covariance technique enables the monitoring of ecosystem carbon fluxes, providing a powerful tool to monitor carbon exchange between the land and atmosphere (Baldocchi, 2003).Recent findings have linked tree-ring width estimated from tree cores and dendrometers with eddy covariance derived NPP showing promise for large scale monitoring of NPP and carbon storage (Babst et al., 2014b;Campioli et al., 2011;Granier et al., 2008;Ohtsuka et al., 2009;Zweifel et al., 2010).This successful linkage of carbon fluxes with tree growth suggests that large scale techniques that may track variation of carbon uptake also reflect variations of annual tree-ring width and thus NPP.
Recent advances have highlighted the potential of satellite-based remote sensing for tracking variations of annual tree-ring width (Coulthard et al., 2017;Kaufmann et al., 2004;Vicente-Serrano et al., 2016).However, it should be noted that there are temporal and spatial constraints of the remote sensing products, that may limit comparisons between remote sensing products and tree-ring width (Babst et al., 2018), with the later having already established long-term tree ring databases like the International Tree-Ring Data Bank (ITRDB).Satellites such as Landsat 5, 7 and 8 contain combined datasets dating back to 1984 with a 30 m spatial resolution, but is limited by poor temporal resolution (16 day revisit) that leads to potential data gaps due to cloud contamination (Wulder et al., 2008).In contrast, MODIS (Moderate Resolution Imaging Spectroradiometer), which was launched in 2000 aboard the Aqua and Terra satellites, provides high temporal resolution (1-2 day revisit) with a 250+ m pixel size spatial resolution.While the spatial resolution of MODIS may be limiting, the high temporal resolution may better capture intra-annual dynamics while also being less susceptible to data gaps from cloud contamination.
MODIS has been used to estimate global annual NPP, derived from 8day MODIS estimates of gross primary productivity (GPP) and net photosynthesis (PSN).These products use a combination of remotely sensed vegetation indices to represent canopy 'greenness' and meteorological data to constrain plant functional variability (Running et al., 2004;Zhao et al., 2005).The MODIS NPP, GPP and PSN products therefore can be used to estimate ecosystem carbon uptake, which may correspond with variations in annual tree-ring width.A recent study by Levesque et al. (2019), demonstrated that satellite-derived NPP (including from MODIS) correlated with tree-ring isotopes, a proxy for reconstructing historical NPP.Given this result, satellite-based remote sensing may be suitable to track variations in tree growth at the ecosystem scale where trees experience similar conditions for parallel growth patterns.
A traditional remotely sensed vegetation index for tracking forest growth, is the normalized difference vegetation index (NDVI) due to its sensitivity to canopy chlorophyll content (i.e.'greenness') and use as a proxy of the fraction of absorbed photosynthetically active radiation and leaf area index (Carlson and Ripley, 1997;Myneni and Williams, 1994).Cumulative NDVI has been observed to correlate with interannual variations of forest growth (Kaufmann et al., 2008;Vicente-Serrano et al., 2016).However, this relationship between NDVI and forest growth may vary by vegetation type and temporal scale (Brehaut and Danby, 2018;Vicente-Serrano et al., 2020).This is especially important for evergreen conifer forests, since vegetation greenness and NDVI (a) have the potential to saturate during the growing season and (b) often remain relatively high even during low productivity periods, decoupling NDVI from photosynthetic activity and carbon accumulation (Sims et al., 2006;Walther et al., 2016).Therefore, we hypothesize that remotely sensed products more closely associated with carbon uptake may better respond to evergreen conifer forest growth than NDVI.
Novel remote sensing products may have potential to overcome the limitations of NDVI for tracking evergreen conifer forest productivity.The near-infrared reflectance of vegetation index (NIR V ), has been recently proposed to reduce saturation effects of NDVI and sensitivity to background conditions (soil and snow) by multiplying NDVI with nearinfrared reflectance, showing close correspondence to ecosystem GPP at large spatiotemporal scales (Badgley et al., 2017).The performance of NIR V requires further validation in evergreen conifer forests, although recent studies have demonstrated its promise (Badgley et al., 2019;Wong et al., 2020).A physiology-based vegetation index was recently described as the chlorophyll/carotenoid index (CCI), which is sensitive to foliar pigment dynamics that are associated with photosynthetic activity (Gamon et al., 2016).In evergreen conifer forests, CCI has been shown to track seasonal dynamics of photosynthesis (Wong et al., 2020(Wong et al., , 2019) ) with further support from its analog, the photochemical reflectance index (PRI) (Garbulsky et al., 2011;Zhang et al., 2016).
Another physiological approach to track the seasonal variation of photosynthesis from satellites is solar-induced fluorescence (SIF) (Frankenberg et al., 2011;Joiner et al., 2011).SIF represents chlorophyll fluorescence emission, which is associated with changes in photosynthetic activity (Porcar-Castell et al., 2014).Strong relationships between SIF and carbon uptake has been observed in evergreen conifer forests (Magney et al., 2019;Walther et al., 2016).However, a current limitation of satellite-based SIF is the poor spatial resolution (40+ km pixel size), which would lead to severe spatial scale mismatch with sitespecific dynamics and tree-ring width comparisons.Recent advances have sought to enhance the spatial resolution of SIF by introducing downscaled SIF products with finer spatial resolutions (~5.5 km pixel size) derived from fusing SIF with additional products from different remote sensing platforms (Duveiller et al., 2020;Wen et al., 2020;Zhang et al., 2018).Therefore, we aim to explore the potential of downscaled SIF products for assessing annual variation of carbon uptake in relation to tree-ring width.
In this study, we evaluated the relationships between annual treering width with remotely sensed MODIS products (NPP, GPP, PSN, NDVI, NIR V and CCI) and downscaled SIF products (CSIF, GOME-2 and SCIAMACHY) across four dominant evergreen conifer species in the Sierra Nevada, spanning 62 sites across a wide range of topographic and environmental conditions.Our primary objectives were to: 1) explore the relationships between tree-ring width and satellite-based products associated with carbon uptake; 2) evaluate a potential legacy effect on these relationships; and 3) characterize the influence of site-specific environmental factors on these relationships across sites.

Study sites
Sampling sites were identified by driving and hiking along roads and trails along an elevational gradient at each of the four focal latitudes (Fig. 1) and identifying relatively homogeneous stands in which Douglas-fir (P.menziesii) was dominant or co-dominant.We sought to identify sites ranging from the upper and lower elevational limits of the distribution of Douglas-fir.Where possible, we sampled a site with a generally north-facing aspect and a site with a generally south-facing aspect at each elevation.Within each sampled site, we identified a random point near the center of the relatively homogeneous area to identify the site location.For each site, topographic information including slope, aspect and elevation were recorded.We avoided areas with clear evidence of moderate-to high-severity wildfire or prescribed fire during the growth period analyzed.

Tree rings/radial growth
Tree growth data was obtained from 794 individual trees across 63 sites in the Sierra Nevada (Fig. 1).Tree diameter at breast height (DBH) measured with metal DBH tapes and tree cores (one per tree) collected with a 4.3 mm-diameter increment borer were sampled in mid-to late summer in 2013 and 2014.From each site, tree measurements were obtained from 3 to 27 individual trees within a roughly 2000 to 4000 m 2 area (Fig. S1).We limited analysis to sites with ≥5 trees sampled resulting in 62 sites in the analysis with a mean of 12.3 trees per site (5.3 standard deviation).
We sampled all co-dominant conifer species.The primary species were P. menziesii, followed by P. ponderosa, A. concolor, and P. lambertiana, which were spatially interspersed within each site.For the least abundant co-dominant species, we sampled the canopy dominant or co-dominant (height at least 67% of the height of neighboring trees) individuals with DBH >20 cm that were most accessible from (generally, closest to) the site center point.We then sampled the closest individual of each other co-dominant species to each of the sampled trees of the least abundant co-dominant species.This was to ensure interspersion of species sampled across the site.Only trees with DBH >20 cm were sampled, and 95% of sampled trees had DBH >30 cm.
There was no minimum stand density for sampling and most stands were not closed-canopy (consistent with most mixed-conifer stands in the Sierra Nevada); we captured an approximation of stand density by computing Voronoi polygon area around each sampled tree (see Section 2.4).Some stands, especially those on south-facing slopes near the lower elevational limit along each transect, contained shrubby understory and/or exposed bedrock.Despite the large spatial range of tree-core measurements at each site, we assume that interannual variability in tree growth is representative of the forest captured within a satellite pixel.
At each site, we collected increment cores from canopy dominant or co-dominant individuals stratified across all species.Tree cores were collected from 1.3 m height on the slope side of the treethe same location used for measuring DBH.After collection and drying, we mounted cores into grooved aluminum blocks, sanded them, and scanned them at 1200 dots per inch (DPI).We measured tree rings in scanned images using CooRecorder v9.01 (Cybis Dendrochronology AB, Sweden) and then performed visual crossdating using CooRecorder and its companion software CDendro following standard crossdating methods (Speer, 2010).We excluded cores and sections of cores that could not be confidently crossdated (e.g., due to growth aberrations such as knots or multiple ambiguous ring boundary decisions that produced similar alignment with the reference chronology).Out of 794 cored trees, 29 were dropped entirely due to uncertain crossdating, yielding 765 cores that contributed to the analysis prior to excluding the site with <5 cored trees (now n = 62).Ring width data were not present for every year in the analysis period (2000 to 2013) for every core (e.g., outer rings broke off on some cores; inner rings were dropped on some due to rot, aberrations such as branches, or uncertain crossdating); the total number of rings analyzed each year was less than 762; it ranged between 533 and 725 (Table S1).Ring width records date from 1956 to 2014.Most trees had annual rings from prior to 1956 but we arbitrarily truncated measurement and crossdating at this year.

Climate data
We obtained annual weather and historical climate normal (1981-2010 mean) data by extracting and summarizing values from the 4 km resolution PRISM dataset (PRISM Climate Group, 2014) for precipitation and the 800 m resolution TopoWx dataset (Oyler et al., 2015) for temperature.For each site, we extracted minimum monthly temperature, maximum monthly temperature, and total monthly precipitation values from the raster data files using bilinear interpolation.We then computed total annual precipitation for each water year (October to September) by summing monthly values, and we computed mean annual temperature for each water year by averaging across monthly values (including averaging minimum and maximum temperature).We computed historical climate normals for temperature and precipitation by calculating the mean value across water years ending in 1981 through 2010.To provide an indicator of drought severity via accumulated water excess or deficit, we used the 4 km pixel Palmer drought severity index (PDSI) from the WestWide Drought Tracker (Abatzoglou et al., 2017) where values range from − 4 to 4 to represent severity where negative is drought, zero is normal and positive is excess water conditions.Daylength was obtained for each site location from the Daymet V3 dataset (Thornton et al., 2016) retrieved using Google Earth Engine.

Estimating local density from Voronoi polygons
We estimated the area of Voronoi polygons (Kenkel, 1991) surrounding each cored tree as an index of growing space, local density, or competition potentially experienced by each tree (Goodwin et al., 2020).Voronoi polygons are created by drawing lines between the focal tree and each of its neighbors, such that each line bisects, and is perpendicular to, the line between the focal tree and its neighbor.The polygon at the interior of these intersecting lines represents the Voronoi polygon (Fig. S2).We only considered trees with a DBH at least half that of the focal tree as potentially contributing to the Voronoi polygon.Trees contributing to the polygon were mapped in the field (by recording their compass azimuth and distance relative to the focal tree) to calculate the Voronoi polygon area.

MODIS data
The MODIS NPP, GPP, PSN, NDVI and NIR V datasets used in this study were obtained from each of the 62 sites using Google Earth Engine to retrieve MODIS data from 2000 to 2014.Annual NPP was acquired from the MOD17A3HGF (V6) collection with a 500 m pixel size.GPP and PSN were acquired from the MOD17A2H (V6) collection consisting of a cumulative 8-day composite with 500 m pixel size.NDVI and NIR V were estimated from the MOD09GQ (V6) collection consisting of daily surface spectral reflectance with a 250 m pixel size and screened for cloud cover using quality flags from the 500 m MOD09GA (V6) collection (Vermote and Wolfe, 2015).We also used the MOD10A1 (V6) normalized difference snow index (NDSI) product obtained via Google Earth Engine to further screen the reflectance data (Hall and Riggs, 2016).For all MODIS vegetation products, dates with NDSI values greater than zero were filtered out (Fig. S3).NDVI and NIR V were calculated using bands 1 (620-670 nm) and 2 (841-876 nm) as: The chlorophyll/carotenoid index (CCI) was calculated from the daily MCD19A1 (V6) MODIS Terra and Aqua combined Multi-Angle Implementation of Atmospheric Correction Land Surface Bidirectional Reflectance Factor gridded Level 2 product (Lyapustin and Wang, 2018).
The 1 km pixel resolution was filtered using quality assessment flags for high quality, cloud and snow free days.Wintertime CCI was sensitive to snow (Fig. S4), therefore we filtered the data using the NDSI where values were greater than zero.CCI was calculated using bands (620-670 nm) and 11 (526-536 nm) as: The MODIS GPP, PSN, NDVI, NIR V and CCI products were averaged monthly and summed annually.To explore the legacy effect, we applied backward monthly shifts ranging from 0 to 12 months for determining the annual sum.For site-specific snow cover, NDSI was used to determine day of year for snow on and snow off, and snow length (number of days with potential snow cover between snow on and snow off).Snow on was determined as the first day of snow when NDSI was >0.Snow off was determined as the first day where NDSI equaled zero for consecutive days.

Solar-induced fluorescence products
Three downscaled SIF products (0.05 • or ~ 5.5 km resolution) were obtained for the 62 sites.The clear-sky instantaneous contiguous SIF (CSIF) dataset, which used SIF from Orbiting Carbon Observatory-2 (OCO-2) combined with surface reflectance from MODIS to generate a 4 day resolution downscaled SIF product from 2000 to 2017 (Zhang et al., 2018).The GOME-2 downscaled SIF dataset is based on Global Ozone Monitoring Experiment 2 (GOME-2) SIF using explanatory variables from MODIS to generate a 8 day resolution downscaled SIF product from 2007 to 2018 (Duveiller et al., 2020).The SCIAMACHY SIF v2.2 dataset is based on the fusion of SIF retrievals from Scanning Imaging Absorption spectrometer for Atmospheric CHartographY (SCIA-MACHY) and GOME-2 for a monthly resolution downscaled SIF product from 2002 to 2018 (Wen et al., 2020).Like the MODIS products, the downscaled SIF products were monthly averages and summed annually, while applying a backward monthly shift ranging from 0 to 12 months for determining the annual sum to explore the legacy effect.

Data analysis
Pearson correlation tests were conducted between the remotely sensed products and annual tree-ring width across years for each site individually to assess the relative correlation coefficients (r) and slopes across sites.Prior to determining the slopes, annual tree-ring width and the remotely sensed products were scaled 0 to 1 for comparisons of slopes.
To explore the interacting relationships between annual tree-ring width, annual remotely sensed products and site-specific climate and geographic characteristics, we used R (R Development Core Team, 2020) to perform a principal component analysis (PCA).To explore the importance of each remotely sensed product, species identity of the cored tree (i.e., P. menziesii, P. ponderosa, A. concolor, and P. lambertiana), site-specific climate and topography in relation to the variation of annual tree-ring width, we used Random Forest analysis using the 'Boruta' package to run iterations until all parameters were stabilized with a maximum iteration of 2000 times to determine variable importance via the mean decrease accuracy measure (Breiman, 2001).

Results
Fig. 2 summarizes the individual site and median of the variation of annual tree-ring width, mean annual temperature, total annual precipitation, PDSI, and the annual sums of the MODIS and downscaled SIF products (for monthly timeseries see Fig. S3).The variation of annual tree-ring width highlights years with larger annual tree-ring width in 2005, 2006, 2010, 2011and 2012 (Fig. 2a (Fig. 2a).These years align with periods of lower mean annual temperature, higher total annual precipitation and abundant soil water via PDSI (Fig. 2b, c, d).For the satellite products, NPP, GPP and PSN showed pronounced annual variability, while NDVI, NIR V , CCI and SIF were relatively stable across all years (Fig. 2e-l).
Across all site and years, there was no relationship between MODIS NPP and annual tree-ring width (Fig. 3).The other satellite products showed slightly higher correlations (0.28 to 0.48) with annual tree-ring width based on calendar year from January to December (zero-month shift) (Fig. S5).Exploring the legacy effect by applying a monthly backwards shift for the annual sum revealed an optimal shift of 5months (August to July) for GPP and PSN, which was consistent across most sites (Fig. 4a, S6) that led to a near three-fold increase in the correlation coefficients (Fig. 4b).The other products (NDVI, NIR V , CCI and SIF) did not show a clear optimal shift across sites (Fig. 4a, S6), although a shift of 6 to 7 months showed a near two-fold increase in the median correlation coefficient (Fig. 4b).Density plots revealed that the site-specific range of correlation coefficients were highest for GPP and PSN with density peaks near 0.75 for GPP and 0.70 for PSN (Fig. 4c).In contrast, the other products showed higher distribution and lower density peaks ranging located from 0.2 to 0.6 (Fig. 4c).The density range of slopes across sites showed distribution peaks for GPP and PSN, while the remaining satellite products were broadly distributed (Fig. 4d).
Applying a 5-month shift legacy effect on the satellite products, the individual site relationships with annual tree-ring width were mostly positive (Fig. 5).When considering the overall relationships across all sites and years, the GPP and annual tree-ring width relationship was nonlinear with larger site-specific slopes and variability in high annual tree-ring width sites than in low annual tree-ring width sites (Fig. 5a).In contrast, the PSN and annual tree-ring width relationships were highly scattered and had a weak overall relationship (Fig. 5b).Both annual tree-ring width and NDVI and NIR V relationships showed a weaker overall relationship compared to the annual tree-ring width and GPP relationship (Fig. 5c, d).There was a poor relationship between CCI and annual tree-ring width (Fig. 5e).The annual tree-ring width and SIF relationships were generally poor with correlation coefficients ranging from 0.24 to 0.38 (Fig. 5f-h).
The PCA biplot revealed the relative correlations between sitespecific climate and geographic characteristics and satellite products (Fig. 6a).A group of variables including all satellite products, DBH and latitude were located near tree-ring width (Fig. 6a).Slightly beyond this group, were precipitation and temperature.To further explore the importance of the explanatory parameters f on the variation of annual tree-ring width, we constructed a Random Forest analysis using all sites and years from 2001 to 2012 with 5-month backwards shifted satellite products (Fig. 6b).The Random Forest model consisting of site-specific For monthly timeseries of satellite products (e-l), see Fig. S3.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)climate, topography and satellite products explained 63% of the variation of annual tree-ring width.The most important explanatory parameters for the variation of annual tree-ring width were GPP and species identity of the sampled tree cores, followed by a grouping of NDVI, local density, NIR V , and CCI in that order.The remaining parameters had importance's of <10% with the weakest predictor being daylength.A second Random Forest model was constructed using only satellite products further highlighting predictor importance (Fig. 6c).Comparing observed and predicted annual tree-ring width from the Random Forest models shows good model performance when all variables are considered including site characteristics and satellite products (r = 0.83) compared to only satellite products (r = 0.38) (Fig. 6d).

Discussion
In this study, we explored the relationships between annual tree-ring width with annual sums of satellite products to represent carbon uptake across 62 evergreen conifer sites in the Sierra Nevada.We found that when a legacy effect is considered (5-month backward shift), the relationship generally improved between tree-ring width and all satellite products.The greatest improvement was for GPP and PSN with a near three-fold increase in the correlation coefficient (r = 0.70 across all sites (Fig. 4b), with the majority of within site variation ranging from r = 0.60 to 0.80).By accounting for this legacy effect on annual carbon uptake, we found that some remotely sensed products of carbon uptake can reveal interannual variation of tree-ring width of evergreen conifers in the Sierra Nevada, but that challenges with newer remote sensing indices such as NIRv, CCI, and downscaled SIF persist.

Satellite estimation of the interannual variation in tree-ring width
Annual tree-ring width varied across both sites and years, which aligned with annual variations in temperature and precipitation (Fig. 2).Using tree-ring width as a ground-based indicator for NPP (Zweifel et al., 2010), we showed that satellite-based NPP, as well as other satellite products as indicators of annual carbon uptake, were unable to resolve interannual variation of tree-ring width (Figs. 3,4,S5).Further exploration revealed a legacy effect for the satellite products and their relationships with annual tree-ring width (Fig. 4b).Legacy effects from the prior year can impact tree growth in following years (Anderegg et al., 2015;Scharnweber et al., 2020).This can lead to a mismatch in timing between annual carbon uptake and resource allocation for growth (Gough et al., 2008).Lagergren et al. (2019) report a legacy effect of a 6-9 months backward shift of carbon uptake in conifers to improve the relationship between annual carbon uptake and annual tree-ring width.For our sites, we identify that a 5-month shift was optimal for GPP and PSN, while it ranged from 6 to 7 months for the remaining products (Fig. 4b, S7).The timing of this shift indicates that carbon uptake during the prior autumn season and subsequent spring season contributes to annual tree-ring width the following year.For GPP and PSN, the 5month shift aligns annual carbon uptake closer to the water year in California, which by convention begins in October (Davis et al., 2016;Flint et al., 2013).This coupling between carbon uptake and water year with annual tree-ring width suggests that water availability, represented by precipitation and PDSI (Fig. 2c, d), is a main limitation for tree growth at our sites, which is common in mountain ecosystems (Dolanc et al., 2013;Littell et al., 2008).Additionally, this timing of carbon uptake and resource availability may reflect the storage of nonstructural carbohydrates during the autumn that can support spring growth resumption (Kozlowski, 1992;Richardson et al., 2013;Tixier et al., 2019).
The 5-month shift strengthened the GPP and PSN relationships with tree-ring width the most, with a nearly three-fold increase in the correlation coefficient (vs.two-fold for NDVI, NIR V , CCI and SIF) (Fig. 4b).The site level performance of the satellite products revealed that GPP and PSN are considerably better than NDVI, NIR V , CCI and SIF for assessing annual tree-ring width based on the range of the correlation coefficients and the higher slope (i.e., greater sensitivity) (Fig. 4c, d).Our results differed from Vicente-Serrano et al. (2020, 2016), who found good correlations between cumulative NDVI and tree-ring growth across evergreen and deciduous forests.Vicente-Serrano et al. (2020) reported relatively low mean NDVI ranging from 0.23 to 0.46, which suggests a potentially large annual dynamic range in foliar greenness.Therefore, we suspect that potential saturation, narrower dynamic range and large site-to-site variation of NDVI, also observed in NIR V , in our sites (Fig. S3c, d) are likely tracking variation in canopy structure like leaf area index (Carlson and Ripley, 1997).This will influence the relationship between tree growth and NDVI, such that NDVI is less sensitive to actual carbon uptake -especially in evergreen conifers which retain their needles year-round (Gamon et al., 1995;Sims et al., 2006;Walther et al., 2016).This can lead to variation in the tree growth and NDVI relationships observed in mountainous ecosystems (Brehaut and Danby, 2018).NIR V performed similarly to NDVI for assessing annual tree-ring width (Fig. 4).Part of this may be due to the sensitivity of NDVI and NIR V to snow cover, which resulted in wintertime data being filtered (Fig. S3c, d), but also suggests that ancillary information on environmental conditions (as in the MODIS GPP, PSN products) is necessary.The lack of seasonality of NDVI and NIR V (Fig. S3c, d) also likely limited the impact of applying a legacy effect as the median correlation coefficient remained relatively similar across all month shifts (Fig. 4b).Further, there is likely decoupling between NIR V and carbon uptake at the beginning and end of the photosynthetically active season (Wang et al., 2020).
The physiology-sensitive satellite products, CCI and SIF, also showed poor relationships with annual tree-ring width even with the consideration of a legacy effect (Figs. 4,5).This was surprising as CCI and SIF generally perform well at tracking annual variations of GPP in evergreen conifer forests (Gamon et al., 2016;Walther et al., 2016).We suspect part of this may be due to spatial scale mismatch as CCI (1 km pixel) and the downscaled SIF products (~5.5 km pixel) were much larger than GPP (500 m).In addition, we note that CCI was influenced by snow cover represented by NDSI resulting in winter months (Nov to Apr) being filtered out from the annual cumulative sums (Fig. S4).These limitations, as well as the large topographic variation within a 1-5 km satellite pixel in the Sierra Nevada, likely severely impacted the CCI and downscaled SIF products and their relationships with annual tree-ring width.
The GPP and PSN products incorporate parameters derived from satellite-based NDVI and climate data including temperature and vapor pressure deficit to account for variation in photosynthetic light-use efficiency (Running et al., 2004).This ultimately leads to a more robust GPP and PSN product for tracking annual carbon uptake compared to NDVI and NIR V on their own, indicated by the differences in their relationships with annual tree-ring width (Fig. 4).Comparing the tree-ring width relationships with GPP and PSN revealed that GPP displayed an overall nonlinear pattern across all sites and years, whereas PSN was more linear and had a lower Pearson correlation coefficient (Fig. 5).This indicates that while PSN reflects local variation in annual tree-ring width, it is unable to resolve these trends across sites.PSN incorporates modeled leaf and root maintenance respiration using lookup tables for modeled temperature and biome properties (Running and  Zhao, 2015).This approach may further reduce sensitivity of the PSN product to local site-specific differences driven by environmental conditions (Neumann et al., 2015;Turner et al., 2006Turner et al., , 2005)).In contrast, the GPP product represents gross carbon uptake, exclusive of respiration, leading to a simpler product that correlates better with annual treering width both within and across sites.This result corresponds with eddy covariance based GPP estimates for tracking variation of annual tree-ring width (Babst et al., 2014b;Lagergren et al., 2019).The satellite-based GPP results imply that interannual variations of carbon uptake correlate with variations of tree-ring width, due to the synchronous response of carbon uptake resource allocation towards wood growth (ie.aboveground NPP).

Site-specific performance
Across the 62 sites used in our analysis, there were site-specific differences in the performance of satellite products for reflecting interannual variation of tree-ring width.GPP and PSN showed most sites had optimal positive relationships with tree-ring width at 5-months (Fig. 4a), although about 2 and 4 sites, respectively, showed negative correlations indicating a poor relationship at the respective sites (Fig. S6).In contrast the other products (NDVI, NIR V , CCI, and SIF) did not have a distinct optimal shift and instead were uniform across many potential month shifts with >8 sites with negative correlations.This is notable as the range of correlation coefficients and slopes was large and did not converge (Fig. 4c, d).Part of this may be driven by the large site-specific variation and range as seen in NDVI, NIR V and CCI suggesting an effect of canopy structure or density (Fig. S3).In contrast, the site-specific variation of GPP, PSN and the SIF products were more uniform in range.
Random forest modeling found that predictions of tree-ring width were substantially improved when site-specific biotic and topographic parameters such as species identity, local density, latitude, slope, elevation and aspect were included (Fig. 6d).Satellite products alone may reflect interannual variation of relative change (ring-width index vs Z-scores, see Fig. S8).This decoupling of tree-ring growth from the fluxrelated parameters estimated by satellite products may indicate that carbon partitioning to wood growth depends on local factors such as resource availability, stand age and stress tolerance, as predicted by optimal partitioning theory (Litton et al., 2007;Poorter et al., 2015Poorter et al., , 2012;;Puglielli et al., 2020;Shipley and Meziane, 2002).While GPP may track variation of annual tree-ring width, for predictive capabilities, other site-specific information is needed.

Limitations
We note that there is a mismatch in spatial scale between tree cores taken from individual trees and satellite pixels, which subsume wide variation among individual trees (Fig. S1).To better match the satellite pixels, we averaged all trees within each 3000 m 2 site and assume these trees are representative of the local region.To better address the spatial scale mismatch, future studies should seek to increase sample size of tree cores across larger spatial areas or investigate this relationship at a smaller spatial scale with proximal sensors.The products with the largest pixels, CCI and SIF, did not perform well in our study, which is likely attributable to the spatial scale mismatch, as these products are only currently available at 1 km and ~ 5.5 km pixels, respectively.We hypothesize that finer resolution products may yield better success in tracking species-and site-specific carbon uptake dynamics for reflecting the variation of annual tree-ring width, as canopy-scale photochemical reflectance index (PRI; analogous to CCI) has shown promise for tracking daily radial growth of evergreen conifers (Eitel et al., 2020).
Further improvements to our understanding of how radial growth depends on the interaction of carbon uptake and environmental factors could advance tree growth models and ultimately improve quantification of long-term aboveground carbon storage (Babst et al., 2018;Voelker, 2011).We note that while we included climate in our analysis, there is also a mismatch between tree cores and climate information such as temperature and PDSI (4 km) and precipitation (800 m).Local meteorological stations would better represent site-specific differences in climate.Site history was not considered fully in our site selection process as fully documenting management and disturbance history would be infeasible given the many jurisdictions and recordkeeping systems involved.Our intent in site selection was to capture a representative sample of disturbance histories, including thinning and prescribed fires, and it is one of the reasons we chose to sample many smaller sites as opposed to fewer larger sites.We did however avoid areas with clear evidence of recent moderate-to high-severity fires.We recognize that the comparison between environmental conditions and tree growth will be tree-age dependent, but in principle this should not impact the reliability of the relationship between a remote sensing derived product and tree-ring growth in a given year.Nonetheless, to adequately link climate conditions with tree-ring growth, future studies need to consider how species-specific growth patterns and age-related growth decline influence the link between carbon uptake and annual tree-ring growth (Babst et al., 2013), as species identity is important in our Random Forest model (Fig. 6b).To explore interannual trends independent of species identity, we also scaled tree-ring width and the satellite products using tree width index and z-scores, respectively, which revealed similar Pearson correlation coefficients with unscaled data (Figs.S8, 5).Therefore, for evaluating interannual trends in forest growth, data scaling may be required across larger geographical ranges and across species.
Standard remote sensing carbon uptake products (eg.GPP and PSN) have inherent limitations as well, leading to generally poor performance for detecting species-and site-specific responses (Turner et al., 2006(Turner et al., , 2005)).This is due to how these products account for photosynthetic efficiency by utilizing biome-specific lookup tables (Ruimy et al., 1994;Running et al., 2004).Physiological-based products such as CCI and SIF may address this gap by providing site-specific information of light harvesting and photoprotective mechanisms as proxies of photosynthetic activity (Gamon et al., 2016;Magney et al., 2019).Therefore, as CCI and SIF products become more readily available and at improved spatial and temporal resolutions, we suggest re-evaluating cumulative sums of CCI and SIF with the legacy effect for reflecting interannual variations of tree-ring growth and aboveground NPP.

Conclusions
We demonstrated that satellite-based GPP with a legacy effect correlated well with variations in annual tree-ring width, both within and across geographically diverse sites in the Sierra Nevada.The legacy effect may reflect a misalignment of the timing of resource gain with resource use; i.e., the use of carbon stores accumulated in one year for growth in the subsequent year (Lagergren et al., 2019).The dynamics of carbon uptake and tree growth in evergreen species may be conserved across other evergreen species in different regions as carbon uptake is often limited by water availability, light availability and temperature (Andreu et al., 2007;Fritts, 1966).This suggests a broader potential for satellite-based detection of the variation of annual tree-ring width.Therefore, the annual carbon uptake and tree-ring width relationships along with the timing of a legacy effect should be validated in other systems, which has already shown promise in a broader context with NDVI (Vicente-Serrano et al., 2016).The relationship between annual tree-ring width and satellite products associated with carbon uptake suggest a promising large-scale tool for tracking evergreen conifer annual tree-ring width variation, an indicator of NPP, to better understand the fate of carbon as climate change will impact the carbon sequestration of forests.

Fig. 1 .
Fig. 1.The location of each evergreen conifer forest site in the Sierra Nevada, California (a) where tree cores were collected and satellite pixels selected (n = 63).Insets represent the north cluster (b) and south cluster (c).

Fig. 2 .
Fig. 2. Boxplot of annual tree-ring width (a), mean annual temperature (b), total annual precipitation (c), mean annual PDSI (d), and annual sums of GPP and median NPP (red diamond) (e), PSN (f), NDVI (g) NIRV (h), normalized CCI (i), SIF CSIF (j), SIF GOME-2 (k) and SIF SCIAMACHY (l).Boxplots represent median, 25th and 75th percentile and interquartile range.Points and colour represent each site (n = 62) and red diamond in (e) represents annual median NPP across all sites.For monthly timeseries of satellite products (e-l), see Fig. S3.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3.No relationship between annual tree-ring width and MODIS NPP.Coloured symbols represent site year data and coloured lines represent site-specific relationships.The black line represents overall relationship of all sites and years with the Pearson correlation coefficient (r).p-value: ns, p > 0.05; *, p < 0.05; **, p < 0.01; *** p < 0.001.

Fig. 4 .
Fig. 4. Density of sites for maximum correlation coefficient across different legacy effects based on backwards shifts from 0 to 12 months (a).Median Pearson correlation coefficients of all sites (n = 62) between annual tree-ring width and annual sum of GPP, NPP, NDVI, NIRV, normalized CCI, SIF CSIF, SIF GOME-2 and SIF SCIAMACHY across different legacy effects (b).Shaded regions indicate standard deviation.Dashed line highlights the optimum 5-month backward shift for GPP and PSN.Density plots of Pearson correlation coefficients (c) and slopes (d) between annual tree-ring width and 5-month backward shift normalized satellite products from all sites.

Fig. 6 .
Fig.6.Principal components analysis of annual data from 62 sites between 2001 and 2012 (a).Variable importance for predicting annual tree-ring width from 2001 to 2012 using Random Forests with 63% explained variation, characterized by percent increase in mean square error where a higher value indicates that a variable is more important to the classification (b) and a Random Forest model using only remotely sensed data (c).Bars represent the iteration median, minimum and maximum range of variable importance.MODIS products include the 5-month backwards shift.Relationship between predicted and observed tree-ring width from the Random Forest models, using only satellite products (RS, red line and symbols) or all variables (ALL, black line and symbols) (d) with the Pearson correlation coefficients (r).p-value: ns, p > 0.05; *, p < 0.05; **, p < 0.01; *** p < 0.001.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)